8#ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
9#define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
20#include <openvdb/thread/Threading.h>
23#include <tbb/parallel_for.h>
24#include <tbb/parallel_sort.h>
25#include <tbb/parallel_invoke.h>
42template<
class Gr
idType>
54template<
class Gr
idType>
64template<
class Gr
idType>
75template<
class Gr
idType>
82template<
typename RealT>
87 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*math::pi<RealT>()*mC), mE(eps) {}
89 inline RealT
operator()(RealT phi)
const {
return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
91 const RealT mC, mD, mE;
102template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
111 static_assert(std::is_floating_point<ValueType>::value,
112 "level set measure is supported only for scalar, floating-point grids");
138 Real area(
bool useWorldUnits =
true);
143 Real volume(
bool useWorldUnits =
true);
148 Real totMeanCurvature(
bool useWorldUnits =
true);
153 Real totGaussianCurvature(
bool useWorldUnits =
true);
158 Real avgMeanCurvature(
bool useWorldUnits =
true) {
return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
163 Real avgGaussianCurvature(
bool useWorldUnits =
true) {
return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
167 int eulerCharacteristic();
172 int genus() {
return 1 - this->eulerCharacteristic()/2;}
176 using LeafT =
typename TreeType::LeafNodeType;
177 using VoxelCIterT =
typename LeafT::ValueOnCIter;
178 using LeafRange =
typename ManagerType::LeafRange;
179 using LeafIterT =
typename LeafRange::Iterator;
180 using ManagerPtr = std::unique_ptr<ManagerType>;
181 using BufferPtr = std::unique_ptr<double[]>;
187 const GridType *mGrid;
190 InterruptT *mInterrupter;
191 double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
193 bool mUpdateArea, mUpdateCurvature;
196 bool checkInterrupter();
200 MeasureArea(
LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
202 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring area and volume of level set");
203 if (parent->mGrainSize>0) {
204 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
206 (*this)(parent->mLeafs->leafRange());
208 tbb::parallel_invoke([&](){parent->mArea = parent->reduce(0);},
209 [&](){parent->mVolume = parent->reduce(1)/3.0;});
210 parent->mUpdateArea =
false;
211 if (parent->mInterrupter) parent->mInterrupter->end();
213 MeasureArea(
const MeasureArea& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
214 void operator()(
const LeafRange& range)
const;
215 LevelSetMeasure* mParent;
216 mutable math::GradStencil<GridT, false> mStencil;
219 struct MeasureCurvatures
221 MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
223 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring curvatures of level set");
224 if (parent->mGrainSize>0) {
225 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
227 (*this)(parent->mLeafs->leafRange());
229 tbb::parallel_invoke([&](){parent->mTotMeanCurvature = parent->reduce(0);},
230 [&](){parent->mTotGausCurvature = parent->reduce(1);});
231 parent->mUpdateCurvature =
false;
232 if (parent->mInterrupter) parent->mInterrupter->end();
234 MeasureCurvatures(
const MeasureCurvatures& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
235 void operator()(
const LeafRange& range)
const;
236 LevelSetMeasure* mParent;
237 mutable math::CurvatureStencil<GridT, false> mStencil;
242 double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
243 tbb::parallel_sort(first, last);
245 while(first != last) sum += *first++;
252template<
typename Gr
idT,
typename InterruptT>
255 : mInterrupter(interrupt)
261template<
typename Gr
idT,
typename InterruptT>
265 if (!grid.hasUniformVoxels()) {
267 "The transform must have uniform scale for the LevelSetMeasure to function");
271 "LevelSetMeasure only supports level sets;"
272 " try setting the grid class to \"level set\"");
276 "LevelSetMeasure does not support empty grids;");
279 mDx = grid.voxelSize()[0];
280 mLeafs = std::make_unique<ManagerType>(mGrid->tree());
281 mBuffer = std::make_unique<double[]>(2*mLeafs->leafCount());
282 mUpdateArea = mUpdateCurvature =
true;
285template<
typename Gr
idT,
typename InterruptT>
289 if (mUpdateArea) {MeasureArea m(
this);};
295template<
typename Gr
idT,
typename InterruptT>
299 if (mUpdateArea) {MeasureArea m(
this);};
300 double volume = mVolume;
301 if (useWorldUnits) volume *=
math::Pow3(mDx) ;
305template<
typename Gr
idT,
typename InterruptT>
309 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
310 return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
313template<
typename Gr
idT,
typename InterruptT>
317 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
318 return mTotGausCurvature;
321template<
typename Gr
idT,
typename InterruptT>
325 const Real x = this->totGaussianCurvature(
true) / (2.0*math::pi<Real>());
331template<
typename Gr
idT,
typename InterruptT>
336 thread::cancelGroupExecution();
342template<
typename Gr
idT,
typename InterruptT>
344LevelSetMeasure<GridT, InterruptT>::
345MeasureArea::operator()(
const LeafRange& range)
const
349 mParent->checkInterrupter();
350 const Real invDx = 1.0/mParent->mDx;
351 const DiracDelta<Real> DD(1.5);
352 const size_t leafCount = mParent->mLeafs->leafCount();
353 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
354 Real sumA = 0, sumV = 0;
355 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
356 const Real dd = DD(invDx * (*voxelIter));
358 mStencil.moveTo(voxelIter);
359 const Coord& p = mStencil.getCenterCoord();
360 const Vec3T g = mStencil.gradient();
361 sumA += dd*g.length();
362 sumV += dd*(g[0]*
Real(p[0]) + g[1]*
Real(p[1]) + g[2]*
Real(p[2]));
365 double* ptr = mParent->mBuffer.get() + leafIter.pos();
372template<
typename Gr
idT,
typename InterruptT>
374LevelSetMeasure<GridT, InterruptT>::
375MeasureCurvatures::operator()(
const LeafRange& range)
const
377 using Vec3T = math::Vec3<ValueType>;
379 mParent->checkInterrupter();
380 const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
381 const DiracDelta<Real> DD(1.5);
382 ValueType mean, gauss;
383 const size_t leafCount = mParent->mLeafs->leafCount();
384 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
385 Real sumM = 0, sumG = 0;
386 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
387 const Real dd = DD(invDx * (*voxelIter));
389 mStencil.moveTo(voxelIter);
390 const Vec3T g = mStencil.gradient();
391 const Real dA = dd*g.length();
392 mStencil.curvatures(mean, gauss);
394 sumG += dA*gauss*dx2;
397 double* ptr = mParent->mBuffer.get() + leafIter.pos();
411typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
412doLevelSetArea(
const GridT& grid,
bool useWorldUnits)
414 LevelSetMeasure<GridT> m(grid);
415 return m.area(useWorldUnits);
420typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
421doLevelSetArea(
const GridT&,
bool)
424 "level set area is supported only for scalar, floating-point grids");
434 return doLevelSetArea<GridT>(grid, useWorldUnits);
444typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
445doLevelSetVolume(
const GridT& grid,
bool useWorldUnits)
447 LevelSetMeasure<GridT> m(grid);
448 return m.volume(useWorldUnits);
453typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
454doLevelSetVolume(
const GridT&,
bool)
457 "level set volume is supported only for scalar, floating-point grids");
467 return doLevelSetVolume<GridT>(grid, useWorldUnits);
477typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
478doLevelSetEulerCharacteristic(
const GridT& grid)
480 LevelSetMeasure<GridT> m(grid);
481 return m.eulerCharacteristic();
486typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
487doLevelSetEulerCharacteristic(
const GridT&)
490 "level set euler characteristic is supported only for scalar, floating-point grids");
501 return doLevelSetEulerCharacteristic(grid);
511typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
512doLevelSetEuler(
const GridT& grid)
514 LevelSetMeasure<GridT> m(grid);
515 return m.eulerCharacteristics();
521typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
522doLevelSetGenus(
const GridT& grid)
524 LevelSetMeasure<GridT> m(grid);
530typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
531doLevelSetGenus(
const GridT&)
534 "level set genus is supported only for scalar, floating-point grids");
544 return doLevelSetGenus(grid);
553#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
555#ifdef OPENVDB_INSTANTIATE_LEVELSETMEASURE
559#define _FUNCTION(TreeT) \
560 Real levelSetArea(const Grid<TreeT>&, bool)
564#define _FUNCTION(TreeT) \
565 Real levelSetVolume(const Grid<TreeT>&, bool)
569#define _FUNCTION(TreeT) \
570 int levelSetEulerCharacteristic(const Grid<TreeT>&)
574#define _FUNCTION(TreeT) \
575 int levelSetGenus(const Grid<TreeT>&)
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Definition: Exceptions.h:63
Definition: Exceptions.h:64
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
T reduce(RangeT range, const T &identity, const FuncT &func, const JoinT &join)
Definition: Reduce.h:42
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:819
Type Pow3(Type x)
Return x3.
Definition: Math.h:552
Type Pow2(Type x)
Return x2.
Definition: Math.h:548
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
double Real
Definition: Types.h:60
@ GRID_LEVEL_SET
Definition: Types.h:416
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:212
#define OPENVDB_INSTANTIATE_CLASS
Definition: version.h.in:153
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:157