8 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
20 #include <boost/math/constants/constants.hpp>
21 #include <tbb/parallel_for.h>
22 #include <tbb/parallel_sort.h>
23 #include <tbb/parallel_invoke.h>
24 #include <type_traits>
39 template<
class Gr
idType>
41 levelSetArea(
const GridType& grid,
bool useWorldSpace =
true);
51 template<
class Gr
idType>
61 template<
class Gr
idType>
72 template<
class Gr
idType>
79 template<
typename RealT>
84 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {}
88 const RealT mC, mD, mE;
99 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
108 static_assert(std::is_floating_point<ValueType>::value,
109 "level set measure is supported only for scalar, floating-point grids");
135 Real area(
bool useWorldUnits =
true);
140 Real volume(
bool useWorldUnits =
true);
145 Real totMeanCurvature(
bool useWorldUnits =
true);
150 Real totGaussianCurvature(
bool useWorldUnits =
true);
155 Real avgMeanCurvature(
bool useWorldUnits =
true) {
return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
160 Real avgGaussianCurvature(
bool useWorldUnits =
true) {
return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
164 int eulerCharacteristic();
169 int genus() {
return 1 - this->eulerCharacteristic()/2;}
173 using LeafT =
typename TreeType::LeafNodeType;
174 using VoxelCIterT =
typename LeafT::ValueOnCIter;
175 using LeafRange =
typename ManagerType::LeafRange;
176 using LeafIterT =
typename LeafRange::Iterator;
177 using ManagerPtr = std::unique_ptr<ManagerType>;
178 using BufferPtr = std::unique_ptr<double[]>;
184 const GridType *mGrid;
187 InterruptT *mInterrupter;
188 double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
190 bool mUpdateArea, mUpdateCurvature;
193 bool checkInterrupter();
197 MeasureArea(
LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
199 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring area and volume of level set");
200 if (parent->mGrainSize>0) {
201 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
203 (*this)(parent->mLeafs->leafRange());
205 tbb::parallel_invoke([&](){parent->mArea = parent->reduce(0);},
206 [&](){parent->mVolume = parent->reduce(1)/3.0;});
207 parent->mUpdateArea =
false;
208 if (parent->mInterrupter) parent->mInterrupter->end();
210 MeasureArea(
const MeasureArea& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
211 void operator()(
const LeafRange& range)
const;
212 LevelSetMeasure* mParent;
213 mutable math::GradStencil<GridT, false> mStencil;
216 struct MeasureCurvatures
218 MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
220 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring curvatures of level set");
221 if (parent->mGrainSize>0) {
222 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
224 (*this)(parent->mLeafs->leafRange());
226 tbb::parallel_invoke([&](){parent->mTotMeanCurvature = parent->reduce(0);},
227 [&](){parent->mTotGausCurvature = parent->reduce(1);});
228 parent->mUpdateCurvature =
false;
229 if (parent->mInterrupter) parent->mInterrupter->end();
231 MeasureCurvatures(
const MeasureCurvatures& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
232 void operator()(
const LeafRange& range)
const;
233 LevelSetMeasure* mParent;
234 mutable math::CurvatureStencil<GridT, false> mStencil;
237 double reduce(
int offset)
239 double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
240 tbb::parallel_sort(first, last);
242 while(first != last) sum += *first++;
249 template<
typename Gr
idT,
typename InterruptT>
252 : mInterrupter(interrupt)
258 template<
typename Gr
idT,
typename InterruptT>
262 if (!grid.hasUniformVoxels()) {
264 "The transform must have uniform scale for the LevelSetMeasure to function");
268 "LevelSetMeasure only supports level sets;"
269 " try setting the grid class to \"level set\"");
273 "LevelSetMeasure does not support empty grids;");
276 mDx = grid.voxelSize()[0];
277 mLeafs = std::make_unique<ManagerType>(mGrid->tree());
278 mBuffer = std::make_unique<double[]>(2*mLeafs->leafCount());
279 mUpdateArea = mUpdateCurvature =
true;
282 template<
typename Gr
idT,
typename InterruptT>
286 if (mUpdateArea) {MeasureArea m(
this);};
292 template<
typename Gr
idT,
typename InterruptT>
296 if (mUpdateArea) {MeasureArea m(
this);};
297 double volume = mVolume;
298 if (useWorldUnits) volume *=
math::Pow3(mDx) ;
302 template<
typename Gr
idT,
typename InterruptT>
306 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
307 return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
310 template<
typename Gr
idT,
typename InterruptT>
314 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
315 return mTotGausCurvature;
318 template<
typename Gr
idT,
typename InterruptT>
322 const Real x = this->totGaussianCurvature(
true) / (2.0*boost::math::constants::pi<Real>());
328 template<
typename Gr
idT,
typename InterruptT>
333 tbb::task::self().cancel_group_execution();
339 template<
typename Gr
idT,
typename InterruptT>
341 LevelSetMeasure<GridT, InterruptT>::
342 MeasureArea::operator()(
const LeafRange& range)
const
346 mParent->checkInterrupter();
347 const Real invDx = 1.0/mParent->mDx;
348 const DiracDelta<Real> DD(1.5);
349 const size_t leafCount = mParent->mLeafs->leafCount();
350 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
351 Real sumA = 0, sumV = 0;
352 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
353 const Real dd = DD(invDx * (*voxelIter));
355 mStencil.moveTo(voxelIter);
356 const Coord& p = mStencil.getCenterCoord();
357 const Vec3T g = mStencil.gradient();
358 sumA += dd*g.length();
359 sumV += dd*(g[0]*
Real(p[0]) + g[1]*
Real(p[1]) + g[2]*
Real(p[2]));
362 double* ptr = mParent->mBuffer.get() + leafIter.pos();
369 template<
typename Gr
idT,
typename InterruptT>
371 LevelSetMeasure<GridT, InterruptT>::
372 MeasureCurvatures::operator()(
const LeafRange& range)
const
374 using Vec3T = math::Vec3<ValueType>;
376 mParent->checkInterrupter();
377 const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
378 const DiracDelta<Real> DD(1.5);
379 ValueType mean, gauss;
380 const size_t leafCount = mParent->mLeafs->leafCount();
381 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
382 Real sumM = 0, sumG = 0;
383 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
384 const Real dd = DD(invDx * (*voxelIter));
386 mStencil.moveTo(voxelIter);
387 const Vec3T g = mStencil.gradient();
388 const Real dA = dd*g.length();
389 mStencil.curvatures(mean, gauss);
391 sumG += dA*gauss*dx2;
394 double* ptr = mParent->mBuffer.get() + leafIter.pos();
406 template<
class Gr
idT>
408 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
409 doLevelSetArea(
const GridT& grid,
bool useWorldUnits)
411 LevelSetMeasure<GridT> m(grid);
412 return m.area(useWorldUnits);
415 template<
class Gr
idT>
417 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
418 doLevelSetArea(
const GridT&,
bool)
421 "level set area is supported only for scalar, floating-point grids");
427 template<
class Gr
idT>
431 return doLevelSetArea<GridT>(grid, useWorldUnits);
439 template<
class Gr
idT>
441 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
442 doLevelSetVolume(
const GridT& grid,
bool useWorldUnits)
444 LevelSetMeasure<GridT> m(grid);
445 return m.volume(useWorldUnits);
448 template<
class Gr
idT>
450 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
451 doLevelSetVolume(
const GridT&,
bool)
454 "level set volume is supported only for scalar, floating-point grids");
460 template<
class Gr
idT>
464 return doLevelSetVolume<GridT>(grid, useWorldUnits);
472 template<
class Gr
idT>
474 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
475 doLevelSetEulerCharacteristic(
const GridT& grid)
477 LevelSetMeasure<GridT> m(grid);
478 return m.eulerCharacteristic();
481 template<
class Gr
idT>
483 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
484 doLevelSetEulerCharacteristic(
const GridT&)
487 "level set euler characteristic is supported only for scalar, floating-point grids");
494 template<
class Gr
idT>
498 return doLevelSetEulerCharacteristic(grid);
506 template<
class Gr
idT>
508 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
509 doLevelSetEuler(
const GridT& grid)
511 LevelSetMeasure<GridT> m(grid);
512 return m.eulerCharacteristics();
516 template<
class Gr
idT>
518 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
519 doLevelSetGenus(
const GridT& grid)
521 LevelSetMeasure<GridT> m(grid);
525 template<
class Gr
idT>
527 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
528 doLevelSetGenus(
const GridT&)
531 "level set genus is supported only for scalar, floating-point grids");
537 template<
class Gr
idT>
541 return doLevelSetGenus(grid);
547 template<
class Gr
idT>
551 bool useWorldUnits =
true)
554 area = m.
area(useWorldUnits);
555 volume = m.
volume(useWorldUnits);
563 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED