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,
141 BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
152 void reinit(
const GridType& grid);
155 void reinit(ManagerType& leafs,
Real dx);
172 void measure(
Real& area,
Real& volume,
bool useWorldUnits =
true);
179 void measure(
Real& area,
Real& volume,
Real& avgMeanCurvature,
bool useWorldUnits =
true);
186 const TreeType* mTree;
188 InterruptT* mInterrupter;
194 bool checkInterrupter();
196 typedef typename TreeType::LeafNodeType LeafT;
197 typedef typename LeafT::ValueOnCIter VoxelCIterT;
198 typedef typename ManagerType::LeafRange LeafRange;
199 typedef typename LeafRange::Iterator LeafIterT;
203 Measure2(
LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
205 if (parent->mGrainSize>0) {
206 tbb::parallel_for(parent->mLeafs->
leafRange(parent->mGrainSize), *
this);
211 Measure2(
const Measure2& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
212 void operator()(
const LeafRange& range)
const;
213 LevelSetMeasure* mParent;
214 typename GridT::ConstAccessor mAcc;
218 Measure3(LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
220 if (parent->mGrainSize>0) {
221 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
223 (*this)(parent->mLeafs->leafRange());
226 Measure3(
const Measure3& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
227 void operator()(
const LeafRange& range)
const;
228 LevelSetMeasure* mParent;
229 typename GridT::ConstAccessor mAcc;
231 inline double reduce(
double* first,
double scale)
233 double* last = first + mLeafs->leafCount();
234 tbb::parallel_sort(first, last);
236 while(first != last) sum += *first++;
243 template<
typename Gr
idT,
typename InterruptT>
246 : mTree(&(grid.tree()))
248 , mInterrupter(interrupt)
249 , mDx(grid.voxelSize()[0])
253 if (!grid.hasUniformVoxels()) {
255 "The transform must have uniform scale for the LevelSetMeasure to function");
259 "LevelSetMeasure only supports level sets;"
260 " try setting the grid class to \"level set\"");
265 template<
typename Gr
idT,
typename InterruptT>
268 ManagerType& leafs,
Real dx, InterruptT* interrupt)
269 : mTree(&(leafs.tree()))
271 , mInterrupter(interrupt)
278 template<
typename Gr
idT,
typename InterruptT>
282 if (!grid.hasUniformVoxels()) {
284 "The transform must have uniform scale for the LevelSetMeasure to function");
288 "LevelSetMeasure only supports level sets;"
289 " try setting the grid class to \"level set\"");
291 mTree = &(grid.tree());
293 mDx = grid.voxelSize()[0];
297 template<
typename Gr
idT,
typename InterruptT>
302 mTree = &(leafs.
tree());
309 template<
typename Gr
idT,
typename InterruptT>
313 if (mInterrupter) mInterrupter->start(
"Measuring level set");
316 const bool newLeafs = mLeafs == NULL;
317 if (newLeafs) mLeafs =
new ManagerType(*mTree);
318 const size_t leafCount = mLeafs->leafCount();
319 if (leafCount == 0) {
323 mArray =
new double[2*leafCount];
327 const double dx = useWorldUnits ? mDx : 1.0;
329 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
337 if (mInterrupter) mInterrupter->end();
341 template<
typename Gr
idT,
typename InterruptT>
344 Real& avgMeanCurvature,
347 if (mInterrupter) mInterrupter->start(
"Measuring level set");
349 const bool newLeafs = mLeafs == NULL;
350 if (newLeafs) mLeafs =
new ManagerType(*mTree);
351 const size_t leafCount = mLeafs->leafCount();
352 if (leafCount == 0) {
353 area = volume = avgMeanCurvature = 0;
356 mArray =
new double[3*leafCount];
360 const double dx = useWorldUnits ? mDx : 1.0;
362 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
363 avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
371 if (mInterrupter) mInterrupter->end();
378 template<
typename Gr
idT,
typename InterruptT>
383 tbb::task::self().cancel_group_execution();
389 template<
typename Gr
idT,
typename InterruptT>
391 LevelSetMeasure<GridT, InterruptT>::
392 Measure2::operator()(
const LeafRange& range)
const
396 mParent->checkInterrupter();
397 const Real invDx = 1.0/mParent->mDx;
398 const DiracDelta<Real> DD(1.5);
399 const size_t leafCount = mParent->mLeafs->leafCount();
400 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
401 Real sumA = 0, sumV = 0;
402 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
403 const Real dd = DD(invDx * (*voxelIter));
405 const Coord p = voxelIter.getCoord();
406 const Vec3T g = invDx*Grad::result(mAcc, p);
407 sumA += dd * g.dot(g);
408 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
411 double* v = mParent->mArray + leafIter.pos();
418 template<
typename Gr
idT,
typename InterruptT>
420 LevelSetMeasure<GridT, InterruptT>::
421 Measure3::operator()(
const LeafRange& range)
const
423 typedef math::Vec3<ValueType> Vec3T;
424 typedef math::ISGradient<math::CD_2ND> Grad;
425 typedef math::ISMeanCurvature<math::CD_SECOND, math::CD_2ND> Curv;
426 mParent->checkInterrupter();
427 const Real invDx = 1.0/mParent->mDx;
428 const DiracDelta<Real> DD(1.5);
429 ValueType alpha, beta;
430 const size_t leafCount = mParent->mLeafs->leafCount();
431 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
432 Real sumA = 0, sumV = 0, sumC = 0;
433 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
434 const Real dd = DD(invDx * (*voxelIter));
436 const Coord p = voxelIter.getCoord();
437 const Vec3T g = invDx*Grad::result(mAcc, p);
438 const Real dA = dd * g.dot(g);
440 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
441 Curv::result(mAcc, p, alpha, beta);
442 sumC += dA * alpha/(2*
math::Pow2(beta))*invDx;
445 double* v = mParent->mArray + leafIter.pos();
456 template<
class Gr
idT>
457 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
462 m.
measure(area, volume, useWorldSpace);
466 template<
class Gr
idT>
467 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
471 "level set area is supported only for scalar, floating-point grids");
474 template<
class Gr
idT>
478 return doLevelSetArea<GridT>(grid, useWorldSpace);
483 template<
class Gr
idT>
484 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
489 m.
measure(area, volume, useWorldSpace);
493 template<
class Gr
idT>
494 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>,
Real>::type
498 "level set volume is supported only for scalar, floating-point grids");
501 template<
class Gr
idT>
505 return doLevelSetVolume<GridT>(grid, useWorldSpace);
510 template<
class Gr
idT>
511 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
515 m.
measure(area, volume, useWorldSpace);
518 template<
class Gr
idT>
519 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
523 "level set measure is supported only for scalar, floating-point grids");
526 template<
class Gr
idT>
530 doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
535 template<
class Gr
idT>
536 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
541 m.
measure(area, volume, avgCurvature, useWorldSpace);
544 template<
class Gr
idT>
545 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
549 "level set measure is supported only for scalar, floating-point grids");
552 template<
class Gr
idT>
556 doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
563 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:293
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:350
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:114
#define OPENVDB_VERSION_NAME
Definition: version.h:43
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:307
Definition: Exceptions.h:86
Definition: Exceptions.h:39
Type Pow2(Type x)
Return .
Definition: Math.h:514
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
double Real
Definition: Types.h:64
Type Pow3(Type x)
Return .
Definition: Math.h:518
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:594
Definition: Exceptions.h:87