8 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
20 #include "openvdb/thread/Threading.h"
22 #include <tbb/parallel_for.h>
23 #include <tbb/parallel_sort.h>
24 #include <tbb/parallel_invoke.h>
26 #include <type_traits>
41 template<
class Gr
idType>
53 template<
class Gr
idType>
63 template<
class Gr
idType>
74 template<
class Gr
idType>
81 template<
typename RealT>
86 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*math::
pi<RealT>()*mC), mE(eps) {}
90 const RealT mC, mD, mE;
101 template<
typename Gr
idT,
typename InterruptT = util::NullInterrupter>
110 static_assert(std::is_floating_point<ValueType>::value,
111 "level set measure is supported only for scalar, floating-point grids");
137 Real area(
bool useWorldUnits =
true);
142 Real volume(
bool useWorldUnits =
true);
147 Real totMeanCurvature(
bool useWorldUnits =
true);
152 Real totGaussianCurvature(
bool useWorldUnits =
true);
157 Real avgMeanCurvature(
bool useWorldUnits =
true) {
return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
162 Real avgGaussianCurvature(
bool useWorldUnits =
true) {
return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
166 int eulerCharacteristic();
171 int genus() {
return 1 - this->eulerCharacteristic()/2;}
175 using LeafT =
typename TreeType::LeafNodeType;
176 using VoxelCIterT =
typename LeafT::ValueOnCIter;
177 using LeafRange =
typename ManagerType::LeafRange;
178 using LeafIterT =
typename LeafRange::Iterator;
179 using ManagerPtr = std::unique_ptr<ManagerType>;
180 using BufferPtr = std::unique_ptr<double[]>;
186 const GridType *mGrid;
189 InterruptT *mInterrupter;
190 double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
192 bool mUpdateArea, mUpdateCurvature;
195 bool checkInterrupter();
199 MeasureArea(
LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
201 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring area and volume of level set");
202 if (parent->mGrainSize>0) {
203 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
205 (*this)(parent->mLeafs->leafRange());
207 tbb::parallel_invoke([&](){parent->mArea = parent->reduce(0);},
208 [&](){parent->mVolume = parent->reduce(1)/3.0;});
209 parent->mUpdateArea =
false;
210 if (parent->mInterrupter) parent->mInterrupter->end();
212 MeasureArea(
const MeasureArea& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
213 void operator()(
const LeafRange& range)
const;
214 LevelSetMeasure* mParent;
215 mutable math::GradStencil<GridT, false> mStencil;
218 struct MeasureCurvatures
220 MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
222 if (parent->mInterrupter) parent->mInterrupter->start(
"Measuring curvatures of level set");
223 if (parent->mGrainSize>0) {
224 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *
this);
226 (*this)(parent->mLeafs->leafRange());
228 tbb::parallel_invoke([&](){parent->mTotMeanCurvature = parent->reduce(0);},
229 [&](){parent->mTotGausCurvature = parent->reduce(1);});
230 parent->mUpdateCurvature =
false;
231 if (parent->mInterrupter) parent->mInterrupter->end();
233 MeasureCurvatures(
const MeasureCurvatures& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
234 void operator()(
const LeafRange& range)
const;
235 LevelSetMeasure* mParent;
236 mutable math::CurvatureStencil<GridT, false> mStencil;
239 double reduce(
int offset)
241 double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
242 tbb::parallel_sort(first, last);
244 while(first != last) sum += *first++;
251 template<
typename Gr
idT,
typename InterruptT>
254 : mInterrupter(interrupt)
260 template<
typename Gr
idT,
typename InterruptT>
264 if (!grid.hasUniformVoxels()) {
266 "The transform must have uniform scale for the LevelSetMeasure to function");
270 "LevelSetMeasure only supports level sets;"
271 " try setting the grid class to \"level set\"");
275 "LevelSetMeasure does not support empty grids;");
278 mDx = grid.voxelSize()[0];
279 mLeafs = std::make_unique<ManagerType>(mGrid->tree());
280 mBuffer = std::make_unique<double[]>(2*mLeafs->leafCount());
281 mUpdateArea = mUpdateCurvature =
true;
284 template<
typename Gr
idT,
typename InterruptT>
288 if (mUpdateArea) {MeasureArea m(
this);};
294 template<
typename Gr
idT,
typename InterruptT>
298 if (mUpdateArea) {MeasureArea m(
this);};
299 double volume = mVolume;
300 if (useWorldUnits) volume *=
math::Pow3(mDx) ;
304 template<
typename Gr
idT,
typename InterruptT>
308 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
309 return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
312 template<
typename Gr
idT,
typename InterruptT>
316 if (mUpdateCurvature) {MeasureCurvatures m(
this);};
317 return mTotGausCurvature;
320 template<
typename Gr
idT,
typename InterruptT>
324 const Real x = this->totGaussianCurvature(
true) / (2.0*math::pi<Real>());
330 template<
typename Gr
idT,
typename InterruptT>
335 thread::cancelGroupExecution();
341 template<
typename Gr
idT,
typename InterruptT>
343 LevelSetMeasure<GridT, InterruptT>::
344 MeasureArea::operator()(
const LeafRange& range)
const
348 mParent->checkInterrupter();
349 const Real invDx = 1.0/mParent->mDx;
350 const DiracDelta<Real> DD(1.5);
351 const size_t leafCount = mParent->mLeafs->leafCount();
352 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
353 Real sumA = 0, sumV = 0;
354 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
355 const Real dd = DD(invDx * (*voxelIter));
357 mStencil.moveTo(voxelIter);
358 const Coord& p = mStencil.getCenterCoord();
359 const Vec3T g = mStencil.gradient();
360 sumA += dd*g.length();
361 sumV += dd*(g[0]*
Real(p[0]) + g[1]*
Real(p[1]) + g[2]*
Real(p[2]));
364 double* ptr = mParent->mBuffer.get() + leafIter.pos();
371 template<
typename Gr
idT,
typename InterruptT>
373 LevelSetMeasure<GridT, InterruptT>::
374 MeasureCurvatures::operator()(
const LeafRange& range)
const
376 using Vec3T = math::Vec3<ValueType>;
378 mParent->checkInterrupter();
379 const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
380 const DiracDelta<Real> DD(1.5);
381 ValueType mean, gauss;
382 const size_t leafCount = mParent->mLeafs->leafCount();
383 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
384 Real sumM = 0, sumG = 0;
385 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
386 const Real dd = DD(invDx * (*voxelIter));
388 mStencil.moveTo(voxelIter);
389 const Vec3T g = mStencil.gradient();
390 const Real dA = dd*g.length();
391 mStencil.curvatures(mean, gauss);
393 sumG += dA*gauss*dx2;
396 double* ptr = mParent->mBuffer.get() + leafIter.pos();
408 template<
class Gr
idT>
410 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
411 doLevelSetArea(
const GridT& grid,
bool useWorldUnits)
413 LevelSetMeasure<GridT> m(grid);
414 return m.area(useWorldUnits);
417 template<
class Gr
idT>
419 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
420 doLevelSetArea(
const GridT&,
bool)
423 "level set area is supported only for scalar, floating-point grids");
429 template<
class Gr
idT>
433 return doLevelSetArea<GridT>(grid, useWorldUnits);
441 template<
class Gr
idT>
443 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
444 doLevelSetVolume(
const GridT& grid,
bool useWorldUnits)
446 LevelSetMeasure<GridT> m(grid);
447 return m.volume(useWorldUnits);
450 template<
class Gr
idT>
452 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
Real>::type
453 doLevelSetVolume(
const GridT&,
bool)
456 "level set volume is supported only for scalar, floating-point grids");
462 template<
class Gr
idT>
466 return doLevelSetVolume<GridT>(grid, useWorldUnits);
474 template<
class Gr
idT>
476 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
477 doLevelSetEulerCharacteristic(
const GridT& grid)
479 LevelSetMeasure<GridT> m(grid);
480 return m.eulerCharacteristic();
483 template<
class Gr
idT>
485 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
486 doLevelSetEulerCharacteristic(
const GridT&)
489 "level set euler characteristic is supported only for scalar, floating-point grids");
496 template<
class Gr
idT>
500 return doLevelSetEulerCharacteristic(grid);
508 template<
class Gr
idT>
510 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
511 doLevelSetEuler(
const GridT& grid)
513 LevelSetMeasure<GridT> m(grid);
514 return m.eulerCharacteristics();
518 template<
class Gr
idT>
520 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value,
int>::type
521 doLevelSetGenus(
const GridT& grid)
523 LevelSetMeasure<GridT> m(grid);
527 template<
class Gr
idT>
529 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value,
int>::type
530 doLevelSetGenus(
const GridT&)
533 "level set genus is supported only for scalar, floating-point grids");
539 template<
class Gr
idT>
543 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: 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
Type Pow3(Type x)
Return x3.
Definition: Math.h:555
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:822
constexpr T pi()
Pi constant taken from Boost to match old behaviour.
Definition: Math.h:118
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
double Real
Definition: Types.h:60
@ GRID_LEVEL_SET
Definition: Types.h:337
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:180