35 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
38 #include <tbb/parallel_for.h>
39 #include <tbb/parallel_sort.h>
40 #include <boost/bind.hpp>
41 #include <boost/function.hpp>
42 #include <boost/type_traits/is_floating_point.hpp>
43 #include <boost/utility/enable_if.hpp>
44 #include <boost/math/constants/constants.hpp>
45 #include <openvdb/math/Math.h>
46 #include <openvdb/Types.h>
47 #include <openvdb/Grid.h>
48 #include <openvdb/tree/LeafManager.h>
49 #include <openvdb/tree/ValueAccessor.h>
50 #include <openvdb/math/FiniteDifference.h>
51 #include <openvdb/math/Operators.h>
52 #include <openvdb/util/NullInterrupter.h>
67 template<
class Gr
idType>
69 levelSetArea(
const GridType& grid,
bool useWorldSpace =
true);
79 template<
class Gr
idType>
93 template<
class Gr
idType>
108 template<
class Gr
idType>
111 bool useWorldSpace =
true);
114 template<
typename RealT>
118 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {}
121 const RealT mC, mD, mE;
132 template<
typename GridT,
142 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
153 void reinit(
const GridType& grid);
156 void reinit(ManagerType& leafs,
Real dx);
173 void measure(
Real& area,
Real& volume,
bool useWorldUnits =
true);
180 void measure(
Real& area,
Real& volume,
Real& avgMeanCurvature,
bool useWorldUnits =
true);
187 if (mTask) mTask(const_cast<LevelSetMeasure*>(
this), range);
192 typedef typename GridT::ConstAccessor AccT;
193 typedef typename TreeType::LeafNodeType LeafT;
194 typedef typename LeafT::ValueOnCIter VoxelCIterT;
195 typedef typename ManagerType::BufferType BufferT;
196 typedef typename RangeType::Iterator LeafIterT;
200 InterruptT* mInterrupter;
203 typename boost::function<void (LevelSetMeasure*, const RangeType&)> mTask;
207 bool checkInterrupter();
210 void measure2(
const RangeType& );
213 void measure3(
const RangeType& );
215 inline double reduce(
double* first,
double scale)
217 double* last = first + mLeafs->leafCount();
218 tbb::parallel_sort(first, last);
220 while(first != last) sum += *first++;
227 template<
typename Gr
idT,
typename InterruptT>
232 , mInterrupter(interrupt)
233 , mDx(grid.voxelSize()[0])
238 if (!grid.hasUniformVoxels()) {
240 "The transform must have uniform scale for the LevelSetMeasure to function");
244 "LevelSetMeasure only supports level sets;"
245 " try setting the grid class to \"level set\"");
250 template<
typename Gr
idT,
typename InterruptT>
253 ManagerType& leafs,
Real dx, InterruptT* interrupt)
256 , mInterrupter(interrupt)
264 template<
typename Gr
idT,
typename InterruptT>
268 if (!grid.hasUniformVoxels()) {
270 "The transform must have uniform scale for the LevelSetMeasure to function");
274 "LevelSetMeasure only supports level sets;"
275 " try setting the grid class to \"level set\"");
278 mAcc = grid.getConstAccessor();
279 mDx = grid.voxelSize()[0];
283 template<
typename Gr
idT,
typename InterruptT>
288 mAcc = AccT(leafs.
tree());
295 template<
typename Gr
idT,
typename InterruptT>
299 if (mInterrupter) mInterrupter->start(
"Measuring level set");
300 mTask = boost::bind(&LevelSetMeasure::measure2, _1, _2);
302 const bool newLeafs = mLeafs == NULL;
303 if (newLeafs) mLeafs =
new ManagerType(mAcc.tree());
304 const size_t leafCount = mLeafs->leafCount();
305 if (leafCount == 0) {
309 mArray =
new double[2*leafCount];
312 tbb::parallel_for(mLeafs->leafRange(mGrainSize), *
this);
314 (*this)(mLeafs->leafRange());
317 const double dx = useWorldUnits ? mDx : 1.0;
319 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
327 if (mInterrupter) mInterrupter->end();
331 template<
typename Gr
idT,
typename InterruptT>
336 if (mInterrupter) mInterrupter->start(
"Measuring level set");
337 mTask = boost::bind(&LevelSetMeasure::measure3, _1, _2);
339 const bool newLeafs = mLeafs == NULL;
340 if (newLeafs) mLeafs =
new ManagerType(mAcc.tree());
341 const size_t leafCount = mLeafs->leafCount();
342 if (leafCount == 0) {
343 area = volume = avgMeanCurvature = 0;
346 mArray =
new double[3*leafCount];
349 tbb::parallel_for(mLeafs->leafRange(mGrainSize), *
this);
351 (*this)(mLeafs->leafRange());
354 const double dx = useWorldUnits ? mDx : 1.0;
356 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
357 avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
365 if (mInterrupter) mInterrupter->end();
372 template<
typename Gr
idT,
typename InterruptT>
377 tbb::task::self().cancel_group_execution();
383 template<
typename Gr
idT,
typename InterruptT>
385 LevelSetMeasure<GridT, InterruptT>::measure2(
const RangeType& range)
389 this->checkInterrupter();
390 const Real invDx = 1.0/mDx;
391 const DiracDelta<Real> DD(1.5);
392 const size_t leafCount = mLeafs->leafCount();
393 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
394 Real sumA = 0, sumV = 0;
395 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
396 const Real dd = DD(invDx * (*voxelIter));
398 const Coord p = voxelIter.getCoord();
399 const Vec3T g = invDx*Grad::result(mAcc, p);
400 sumA += dd * g.dot(g);
401 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
404 double* v = mArray + leafIter.pos();
412 template<
typename Gr
idT,
typename InterruptT>
414 LevelSetMeasure<GridT, InterruptT>::measure3(
const RangeType& range)
416 typedef math::Vec3<ValueType> Vec3T;
417 typedef math::ISGradient<math::CD_2ND> Grad;
418 typedef math::ISMeanCurvature<math::CD_SECOND, math::CD_2ND> Curv;
419 this->checkInterrupter();
420 const Real invDx = 1.0/mDx;
421 const DiracDelta<Real> DD(1.5);
422 ValueType alpha, beta;
423 const size_t leafCount = mLeafs->leafCount();
424 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
425 Real sumA = 0, sumV = 0, sumC = 0;
426 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
427 const Real dd = DD(invDx * (*voxelIter));
429 const Coord p = voxelIter.getCoord();
430 const Vec3T g = invDx*Grad::result(mAcc, p);
431 const Real dA = dd * g.dot(g);
433 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
434 Curv::result(mAcc, p, alpha, beta);
435 sumC += dA * alpha/(2*
math::Pow2(beta))*invDx;
438 double* v = mArray + leafIter.pos();
449 template<
class Gr
idT>
450 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
455 m.
measure(area, volume, useWorldSpace);
459 template<
class Gr
idT>
460 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
464 "level set area is supported only for scalar, floating-point grids");
467 template<
class Gr
idT>
471 return doLevelSetArea<GridT>(grid, useWorldSpace);
476 template<
class Gr
idT>
477 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
482 m.
measure(area, volume, useWorldSpace);
486 template<
class Gr
idT>
487 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
491 "level set volume is supported only for scalar, floating-point grids");
494 template<
class Gr
idT>
498 return doLevelSetVolume<GridT>(grid, useWorldSpace);
503 template<
class Gr
idT>
504 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
508 m.
measure(area, volume, useWorldSpace);
511 template<
class Gr
idT>
512 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
516 "level set measure is supported only for scalar, floating-point grids");
519 template<
class Gr
idT>
523 doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
528 template<
class Gr
idT>
529 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
534 m.
measure(area, volume, avgCurvature, useWorldSpace);
537 template<
class Gr
idT>
538 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
542 "level set measure is supported only for scalar, floating-point grids");
545 template<
class Gr
idT>
549 doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
556 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
Definition: LeafManager.h:126
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:284
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Type Pow2(Type x)
Return .
Definition: Math.h:498
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:109
double Real
Definition: Types.h:65
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Definition: Exceptions.h:86
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Definition: Exceptions.h:87
Definition: Exceptions.h:39
Type Pow3(Type x)
Return .
Definition: Math.h:502
Definition: Exceptions.h:88
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:594
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:302
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71