35 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED 36 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED 38 #include <openvdb/math/Math.h> 39 #include <openvdb/Types.h> 40 #include <openvdb/Grid.h> 41 #include <openvdb/tree/LeafManager.h> 42 #include <openvdb/tree/ValueAccessor.h> 43 #include <openvdb/math/FiniteDifference.h> 44 #include <openvdb/math/Operators.h> 45 #include <openvdb/util/NullInterrupter.h> 46 #include <boost/math/constants/constants.hpp> 47 #include <tbb/parallel_for.h> 48 #include <tbb/parallel_sort.h> 49 #include <type_traits> 65 template<
class Gr
idType>
67 levelSetArea(
const GridType& grid,
bool useWorldSpace =
true);
77 template<
class Gr
idType>
91 template<
class Gr
idType>
106 template<
class Gr
idType>
109 bool useWorldSpace =
true);
112 template<
typename RealT>
116 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {}
119 const RealT mC, mD, mE;
130 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
138 static_assert(std::is_floating_point<ValueType>::value,
139 "level set measure is supported only for scalar, floating-point grids");
170 void measure(
Real& area,
Real& volume,
bool useWorldUnits =
true);
177 void measure(
Real& area,
Real& volume,
Real& avgMeanCurvature,
bool useWorldUnits =
true);
186 InterruptT* mInterrupter;
192 bool checkInterrupter();
194 using LeafT =
typename TreeType::LeafNodeType;
195 using VoxelCIterT =
typename LeafT::ValueOnCIter;
196 using LeafRange =
typename ManagerType::LeafRange;
197 using LeafIterT =
typename LeafRange::Iterator;
201 Measure2(
LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
203 if (parent->mGrainSize>0) {
204 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
206 (*this)(parent->mLeafs->leafRange());
209 Measure2(
const Measure2& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
210 void operator()(
const LeafRange& range)
const;
212 typename GridT::ConstAccessor mAcc;
216 Measure3(
LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
218 if (parent->mGrainSize>0) {
219 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
221 (*this)(parent->mLeafs->leafRange());
224 Measure3(
const Measure3& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
225 void operator()(
const LeafRange& range)
const;
227 typename GridT::ConstAccessor mAcc;
229 inline double reduce(
double* first,
double scale)
231 double* last = first + mLeafs->leafCount();
232 tbb::parallel_sort(first, last);
234 while(first != last) sum += *first++;
241 template<
typename Gr
idT,
typename InterruptT>
244 : mTree(&(grid.tree()))
246 , mInterrupter(interrupt)
247 , mDx(grid.voxelSize()[0])
251 if (!grid.hasUniformVoxels()) {
253 "The transform must have uniform scale for the LevelSetMeasure to function");
257 "LevelSetMeasure only supports level sets;" 258 " try setting the grid class to \"level set\"");
263 template<
typename Gr
idT,
typename InterruptT>
267 : mTree(&(leafs.tree()))
269 , mInterrupter(interrupt)
276 template<
typename Gr
idT,
typename InterruptT>
280 if (!grid.hasUniformVoxels()) {
282 "The transform must have uniform scale for the LevelSetMeasure to function");
286 "LevelSetMeasure only supports level sets;" 287 " try setting the grid class to \"level set\"");
289 mTree = &(grid.tree());
291 mDx = grid.voxelSize()[0];
295 template<
typename Gr
idT,
typename InterruptT>
300 mTree = &(leafs.tree());
307 template<
typename Gr
idT,
typename InterruptT>
311 if (mInterrupter) mInterrupter->start(
"Measuring level set");
314 const bool newLeafs = mLeafs ==
nullptr;
316 const size_t leafCount = mLeafs->leafCount();
317 if (leafCount == 0) {
321 mArray =
new double[2*leafCount];
325 const double dx = useWorldUnits ? mDx : 1.0;
327 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
335 if (mInterrupter) mInterrupter->end();
339 template<
typename Gr
idT,
typename InterruptT>
342 Real& avgMeanCurvature,
345 if (mInterrupter) mInterrupter->start(
"Measuring level set");
347 const bool newLeafs = mLeafs ==
nullptr;
349 const size_t leafCount = mLeafs->leafCount();
350 if (leafCount == 0) {
351 area = volume = avgMeanCurvature = 0;
354 mArray =
new double[3*leafCount];
358 const double dx = useWorldUnits ? mDx : 1.0;
360 volume = this->reduce(mArray + leafCount,
math::Pow3(dx) / 3.0);
361 avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
369 if (mInterrupter) mInterrupter->end();
376 template<
typename Gr
idT,
typename InterruptT>
381 tbb::task::self().cancel_group_execution();
387 template<
typename Gr
idT,
typename InterruptT>
394 mParent->checkInterrupter();
395 const Real invDx = 1.0/mParent->mDx;
397 const size_t leafCount = mParent->mLeafs->leafCount();
398 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
399 Real sumA = 0, sumV = 0;
400 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
401 const Real dd = DD(invDx * (*voxelIter));
403 const Coord p = voxelIter.getCoord();
404 const Vec3T g = invDx*Grad::result(mAcc, p);
405 sumA += dd * g.dot(g);
406 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
409 double* v = mParent->mArray + leafIter.pos();
416 template<
typename Gr
idT,
typename InterruptT>
424 mParent->checkInterrupter();
425 const Real invDx = 1.0/mParent->mDx;
428 const size_t leafCount = mParent->mLeafs->leafCount();
429 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
430 Real sumA = 0, sumV = 0, sumC = 0;
431 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
432 const Real dd = DD(invDx * (*voxelIter));
434 const Coord p = voxelIter.getCoord();
435 const Vec3T g = invDx*Grad::result(mAcc, p);
436 const Real dA = dd * g.dot(g);
438 sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
439 Curv::result(mAcc, p, alpha, beta);
440 sumC += dA * alpha/(2*
math::Pow2(beta))*invDx;
443 double* v = mParent->mArray + leafIter.pos();
459 template<
class Gr
idT>
461 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
462 doLevelSetArea(
const GridT& grid,
bool useWorldSpace)
466 m.
measure(area, volume, useWorldSpace);
470 template<
class Gr
idT>
472 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
473 doLevelSetArea(
const GridT&,
bool)
476 "level set area is supported only for scalar, floating-point grids");
483 template<
class Gr
idT>
487 return doLevelSetArea<GridT>(grid, useWorldSpace);
497 template<
class Gr
idT>
499 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
500 doLevelSetVolume(
const GridT& grid,
bool useWorldSpace)
504 m.
measure(area, volume, useWorldSpace);
508 template<
class Gr
idT>
510 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
511 doLevelSetVolume(
const GridT&,
bool)
514 "level set volume is supported only for scalar, floating-point grids");
521 template<
class Gr
idT>
525 return doLevelSetVolume<GridT>(grid, useWorldSpace);
535 template<
class Gr
idT>
537 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value>::type
538 doLevelSetMeasure(
const GridT& grid,
Real& area,
Real& volume,
bool useWorldSpace)
541 m.
measure(area, volume, useWorldSpace);
544 template<
class Gr
idT>
546 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value>::type
547 doLevelSetMeasure(
const GridT&,
Real&,
Real&,
bool)
550 "level set measure is supported only for scalar, floating-point grids");
557 template<
class Gr
idT>
561 doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
571 template<
class Gr
idT>
573 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value>::type
574 doLevelSetMeasure(
const GridT& grid,
Real& area,
Real& volume,
Real& avgCurvature,
578 m.
measure(area, volume, avgCurvature, useWorldSpace);
581 template<
class Gr
idT>
583 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value>::type
584 doLevelSetMeasure(
const GridT&,
Real&,
Real&,
Real&,
bool)
587 "level set measure is supported only for scalar, floating-point grids");
594 template<
class Gr
idT>
598 doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
605 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
double Real
Definition: Types.h:67
Compute the mean curvature in index space.
Definition: Operators.h:556
Gradient operators defined in index space of various orders.
Definition: Operators.h:123
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:51
Definition: Exceptions.h:91
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:647
Definition: Exceptions.h:40
Definition: Exceptions.h:90
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:110
Type Pow3(Type x)
Return x3.
Definition: Math.h:506
Type Pow2(Type x)
Return x2.
Definition: Math.h:502
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188