8 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
20 #include <tbb/parallel_for.h>
21 #include <tbb/parallel_sort.h>
22 #include <tbb/parallel_invoke.h>
23 #include <type_traits>
38 template<
class Gr
idType>
50 template<
class Gr
idType>
60 template<
class Gr
idType>
71 template<
class Gr
idType>
78 template<
typename RealT>
83 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*math::
pi<RealT>()*mC), mE(eps) {}
87 const RealT mC, mD, mE;
98 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
107 static_assert(std::is_floating_point<ValueType>::value,
108 "level set measure is supported only for scalar, floating-point grids");
134 Real area(
bool useWorldUnits =
true);
139 Real volume(
bool useWorldUnits =
true);
144 Real totMeanCurvature(
bool useWorldUnits =
true);
149 Real totGaussianCurvature(
bool useWorldUnits =
true);
154 Real avgMeanCurvature(
bool useWorldUnits =
true) {
return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
159 Real avgGaussianCurvature(
bool useWorldUnits =
true) {
return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
163 int eulerCharacteristic();
168 int genus() {
return 1 - this->eulerCharacteristic()/2;}
172 using LeafT =
typename TreeType::LeafNodeType;
173 using VoxelCIterT =
typename LeafT::ValueOnCIter;
174 using LeafRange =
typename ManagerType::LeafRange;
175 using LeafIterT =
typename LeafRange::Iterator;
176 using ManagerPtr = std::unique_ptr<ManagerType>;
177 using BufferPtr = std::unique_ptr<double[]>;
183 const GridType *mGrid;
186 InterruptT *mInterrupter;
187 double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
189 bool mUpdateArea, mUpdateCurvature;
192 bool checkInterrupter();
196 MeasureArea(
LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
198 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring area and volume of level set");
199 if (parent->mGrainSize>0) {
200 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
202 (*this)(parent->mLeafs->leafRange());
204 tbb::parallel_invoke([&](){parent->mArea = parent->reduce(0);},
205 [&](){parent->mVolume = parent->reduce(1)/3.0;});
206 parent->mUpdateArea =
false;
207 if (parent->mInterrupter) parent->mInterrupter->end();
209 MeasureArea(
const MeasureArea& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
210 void operator()(
const LeafRange& range)
const;
211 LevelSetMeasure* mParent;
212 mutable math::GradStencil<GridT, false> mStencil;
215 struct MeasureCurvatures
217 MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
219 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring curvatures of level set");
220 if (parent->mGrainSize>0) {
221 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
223 (*this)(parent->mLeafs->leafRange());
225 tbb::parallel_invoke([&](){parent->mTotMeanCurvature = parent->reduce(0);},
226 [&](){parent->mTotGausCurvature = parent->reduce(1);});
227 parent->mUpdateCurvature =
false;
228 if (parent->mInterrupter) parent->mInterrupter->end();
230 MeasureCurvatures(
const MeasureCurvatures& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
231 void operator()(
const LeafRange& range)
const;
232 LevelSetMeasure* mParent;
233 mutable math::CurvatureStencil<GridT, false> mStencil;
236 double reduce(
int offset)
238 double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
239 tbb::parallel_sort(first, last);
241 while(first != last) sum += *first++;
248 template<
typename Gr
idT,
typename InterruptT>
251 : mInterrupter(interrupt)
257 template<
typename Gr
idT,
typename InterruptT>
261 if (!grid.hasUniformVoxels()) {
263 "The transform must have uniform scale for the LevelSetMeasure to function");
267 "LevelSetMeasure only supports level sets;"
268 " try setting the grid class to \"level set\"");
272 "LevelSetMeasure does not support empty grids;");
275 mDx = grid.voxelSize()[0];
276 mLeafs = std::make_unique<ManagerType>(mGrid->tree());
277 mBuffer = std::make_unique<double[]>(2*mLeafs->leafCount());
278 mUpdateArea = mUpdateCurvature =
true;
281 template<
typename Gr
idT,
typename InterruptT>
285 if (mUpdateArea) {MeasureArea m(
this);};
291 template<
typename Gr
idT,
typename InterruptT>
295 if (mUpdateArea) {MeasureArea m(
this);};
296 double volume = mVolume;
297 if (useWorldUnits) volume *=
math::Pow3(mDx) ;
301 template<
typename Gr
idT,
typename InterruptT>
305 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
306 return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
309 template<
typename Gr
idT,
typename InterruptT>
313 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
314 return mTotGausCurvature;
317 template<
typename Gr
idT,
typename InterruptT>
321 const Real x = this->totGaussianCurvature(
true) / (2.0*math::pi<Real>());
327 template<
typename Gr
idT,
typename InterruptT>
332 tbb::task::self().cancel_group_execution();
338 template<
typename Gr
idT,
typename InterruptT>
340 LevelSetMeasure<GridT, InterruptT>::
341 MeasureArea::operator()(
const LeafRange& range)
const
345 mParent->checkInterrupter();
346 const Real invDx = 1.0/mParent->mDx;
347 const DiracDelta<Real> DD(1.5);
348 const size_t leafCount = mParent->mLeafs->leafCount();
349 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
350 Real sumA = 0, sumV = 0;
351 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
352 const Real dd = DD(invDx * (*voxelIter));
354 mStencil.moveTo(voxelIter);
355 const Coord& p = mStencil.getCenterCoord();
356 const Vec3T g = mStencil.gradient();
357 sumA += dd*g.length();
358 sumV += dd*(g[0]*
Real(p[0]) + g[1]*
Real(p[1]) + g[2]*
Real(p[2]));
361 double* ptr = mParent->mBuffer.get() + leafIter.pos();
368 template<
typename Gr
idT,
typename InterruptT>
370 LevelSetMeasure<GridT, InterruptT>::
371 MeasureCurvatures::operator()(
const LeafRange& range)
const
373 using Vec3T = math::Vec3<ValueType>;
375 mParent->checkInterrupter();
376 const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
377 const DiracDelta<Real> DD(1.5);
378 ValueType mean, gauss;
379 const size_t leafCount = mParent->mLeafs->leafCount();
380 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
381 Real sumM = 0, sumG = 0;
382 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
383 const Real dd = DD(invDx * (*voxelIter));
385 mStencil.moveTo(voxelIter);
386 const Vec3T g = mStencil.gradient();
387 const Real dA = dd*g.length();
388 mStencil.curvatures(mean, gauss);
390 sumG += dA*gauss*dx2;
393 double* ptr = mParent->mBuffer.get() + leafIter.pos();
405 template<
class Gr
idT>
407 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
408 doLevelSetArea(
const GridT& grid,
bool useWorldUnits)
410 LevelSetMeasure<GridT> m(grid);
411 return m.area(useWorldUnits);
414 template<
class Gr
idT>
416 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
417 doLevelSetArea(
const GridT&,
bool)
420 "level set area is supported only for scalar, floating-point grids");
426 template<
class Gr
idT>
430 return doLevelSetArea<GridT>(grid, useWorldUnits);
438 template<
class Gr
idT>
440 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
441 doLevelSetVolume(
const GridT& grid,
bool useWorldUnits)
443 LevelSetMeasure<GridT> m(grid);
444 return m.volume(useWorldUnits);
447 template<
class Gr
idT>
449 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
450 doLevelSetVolume(
const GridT&,
bool)
453 "level set volume is supported only for scalar, floating-point grids");
459 template<
class Gr
idT>
463 return doLevelSetVolume<GridT>(grid, useWorldUnits);
471 template<
class Gr
idT>
473 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
474 doLevelSetEulerCharacteristic(
const GridT& grid)
476 LevelSetMeasure<GridT> m(grid);
477 return m.eulerCharacteristic();
480 template<
class Gr
idT>
482 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
483 doLevelSetEulerCharacteristic(
const GridT&)
486 "level set euler characteristic is supported only for scalar, floating-point grids");
493 template<
class Gr
idT>
497 return doLevelSetEulerCharacteristic(grid);
505 template<
class Gr
idT>
507 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
508 doLevelSetEuler(
const GridT& grid)
510 LevelSetMeasure<GridT> m(grid);
511 return m.eulerCharacteristics();
515 template<
class Gr
idT>
517 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
518 doLevelSetGenus(
const GridT& grid)
520 LevelSetMeasure<GridT> m(grid);
524 template<
class Gr
idT>
526 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
527 doLevelSetGenus(
const GridT&)
530 "level set genus is supported only for scalar, floating-point grids");
536 template<
class Gr
idT>
540 return doLevelSetGenus(grid);
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...
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
Definition: openvdb/Exceptions.h:63
Definition: openvdb/Exceptions.h:64
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:26
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:85
Type Pow3(Type x)
Return x3.
Definition: Math.h:555
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:822
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
constexpr T pi()
Pi constant taken from Boost to match old behaviour.
Definition: Math.h:118
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
double Real
Definition: openvdb/Types.h:38
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:315
Definition: openvdb/Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:101
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:153