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;
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>
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...
ValueAccessors are designed to help accelerate accesses into the OpenVDB Tree structures by storing c...
Definition Exceptions.h:63
Definition Exceptions.h:64
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
Definition Stencils.h:1232
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
constexpr T pi()
Pi constant taken from Boost to match old behaviour.
Definition Math.h:119
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:455
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