HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevelSetMeasure.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
4 /// @author Ken Museth
5 ///
6 /// @file LevelSetMeasure.h
7 
8 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
9 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
10 
11 #include "openvdb/Types.h"
12 #include "openvdb/Grid.h"
15 #include "openvdb/math/Math.h"
17 #include "openvdb/math/Operators.h"
18 #include "openvdb/math/Stencils.h"
21 
22 #include <tbb/parallel_for.h>
23 #include <tbb/parallel_sort.h>
24 #include <tbb/parallel_invoke.h>
25 
26 #include <type_traits>
27 
28 namespace openvdb {
30 namespace OPENVDB_VERSION_NAME {
31 namespace tools {
32 
33 /// @brief Return the surface area of a narrow-band level set.
34 ///
35 /// @param grid a scalar, floating-point grid with one or more disjoint,
36 /// closed level set surfaces
37 /// @param useWorldSpace if true the area is computed in
38 /// world space units, else in voxel units.
39 ///
40 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
41 template<class GridType>
42 inline Real
43 levelSetArea(const GridType& grid, bool useWorldSpace = true);
44 
45 /// @brief Return the volume of a narrow-band level set surface.
46 ///
47 /// @param grid a scalar, floating-point grid with one or more disjoint,
48 /// closed level set surfaces
49 /// @param useWorldSpace if true the volume is computed in
50 /// world space units, else in voxel units.
51 ///
52 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
53 template<class GridType>
54 inline Real
55 levelSetVolume(const GridType& grid, bool useWorldSpace = true);
56 
57 /// @brief Return the Euler Characteristics of a narrow-band level set surface (possibly disconnected).
58 ///
59 /// @param grid a scalar, floating-point grid with one or more disjoint,
60 /// closed level set surfaces
61 ///
62 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
63 template<class GridType>
64 inline int
65 levelSetEulerCharacteristic(const GridType& grid);
66 
67 /// @brief Return the genus of a narrow-band level set surface.
68 ///
69 /// @param grid a scalar, floating-point grid with one or more disjoint,
70 /// closed level set surfaces
71 /// @warning The genus is only well defined for a single connected surface
72 ///
73 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
74 template<class GridType>
75 inline int
76 levelSetGenus(const GridType& grid);
77 
78 ////////////////////////////////////////////////////////////////////////////////////////
79 
80 /// @brief Smeared-out and continuous Dirac Delta function.
81 template<typename RealT>
83 {
84 public:
85  // eps is the half-width of the dirac delta function in units of phi
86  DiracDelta(RealT eps) : mC(0.5/eps), mD(2*math::pi<RealT>()*mC), mE(eps) {}
87  // values of the dirac delta function are in units of one over the units of phi
88  inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
89 private:
90  const RealT mC, mD, mE;
91 };// DiracDelta functor
92 
93 
94 /// @brief Multi-threaded computation of surface area, volume and
95 /// average mean-curvature for narrow band level sets.
96 ///
97 /// @details To reduce the risk of round-off errors (primarily due to
98 /// catastrophic cancellation) and guarantee determinism during
99 /// multi-threading this class is implemented using parallel_for, and
100 /// delayed reduction of a sorted list.
101 template<typename GridT, typename InterruptT = util::NullInterrupter>
103 {
104 public:
105  using GridType = GridT;
106  using TreeType = typename GridType::TreeType;
107  using ValueType = typename TreeType::ValueType;
109 
111  "level set measure is supported only for scalar, floating-point grids");
112 
113  /// @brief Main constructor from a grid
114  /// @param grid The level set to be measured.
115  /// @param interrupt Optional interrupter.
116  /// @throw RuntimeError if the grid is not a level set or if it's empty.
117  LevelSetMeasure(const GridType& grid, InterruptT* interrupt = nullptr);
118 
119  /// @brief Re-initialize using the specified grid.
120  /// @param grid The level set to be measured.
121  /// @throw RuntimeError if the grid is not a level set or if it's empty.
122  void init(const GridType& grid);
123 
124  /// @brief Destructor
125  virtual ~LevelSetMeasure() {}
126 
127  /// @return the grain-size used for multi-threading
128  int getGrainSize() const { return mGrainSize; }
129 
130  /// @brief Set the grain-size used for multi-threading.
131  /// @note A grain size of 0 or less disables multi-threading!
132  void setGrainSize(int grainsize) { mGrainSize = grainsize; }
133 
134  /// @brief Compute the surface area of the level set.
135  /// @param useWorldUnits Specifies if the result is in world or voxel units.
136  /// @note Performs internal caching so only the initial call incurs actual computation.
137  Real area(bool useWorldUnits = true);
138 
139  /// @brief Compute the volume of the level set surface.
140  /// @param useWorldUnits Specifies if the result is in world or voxel units.
141  /// @note Performs internal caching so only the initial call incurs actual computation.
142  Real volume(bool useWorldUnits = true);
143 
144  /// @brief Compute the total mean curvature of the level set surface.
145  /// @param useWorldUnits Specifies if the result is in world or voxel units.
146  /// @note Performs internal caching so only the initial call incurs actual computation.
147  Real totMeanCurvature(bool useWorldUnits = true);
148 
149  /// @brief Compute the total gaussian curvature of the level set surface.
150  /// @param useWorldUnits Specifies if the result is in world or voxel units.
151  /// @note Performs internal caching so only the initial call incurs actual computation.
152  Real totGaussianCurvature(bool useWorldUnits = true);
153 
154  /// @brief Compute the average mean curvature of the level set surface.
155  /// @param useWorldUnits Specifies if the result is in world or voxel units.
156  /// @note Performs internal caching so only the initial call incurs actual computation.
157  Real avgMeanCurvature(bool useWorldUnits = true) {return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
158 
159  /// @brief Compute the average gaussian curvature of the level set surface.
160  /// @param useWorldUnits Specifies if the result is in world or voxel units.
161  /// @note Performs internal caching so only the initial call incurs actual computation.
162  Real avgGaussianCurvature(bool useWorldUnits = true) {return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
163 
164  /// @brief Compute the Euler characteristic of the level set surface.
165  /// @note Performs internal caching so only the initial call incurs actual computation.
166  int eulerCharacteristic();
167 
168  /// @brief Compute the genus of the level set surface.
169  /// @warning The genus is only well defined for a single connected surface.
170  /// @note Performs internal caching so only the initial call incurs actual computation.
171  int genus() { return 1 - this->eulerCharacteristic()/2;}
172 
173 private:
174 
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[]>;
181 
182  // disallow copy construction and copy by assignment!
183  LevelSetMeasure(const LevelSetMeasure&);// not implemented
184  LevelSetMeasure& operator=(const LevelSetMeasure&);// not implemented
185 
186  const GridType *mGrid;
187  ManagerPtr mLeafs;
188  BufferPtr mBuffer;
189  InterruptT *mInterrupter;
190  double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
191  int mGrainSize;
192  bool mUpdateArea, mUpdateCurvature;
193 
194  // @brief Return false if the process was interrupted
195  bool checkInterrupter();
196 
197  struct MeasureArea
198  {
199  MeasureArea(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
200  {
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);
204  } else {
205  (*this)(parent->mLeafs->leafRange());
206  }
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();
211  }
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;
216  };// MeasureArea
217 
218  struct MeasureCurvatures
219  {
220  MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
221  {
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);
225  } else {
226  (*this)(parent->mLeafs->leafRange());
227  }
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();
232  }
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;
237  };// MeasureCurvatures
238 
239  double reduce(int offset)
240  {
241  double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
242  tbb::parallel_sort(first, last);// mitigates catastrophic cancellation
243  Real sum = 0.0;
244  while(first != last) sum += *first++;
245  return sum;
246  }
247 
248 }; // end of LevelSetMeasure class
249 
250 
251 template<typename GridT, typename InterruptT>
252 inline
254  : mInterrupter(interrupt)
255  , mGrainSize(1)
256 {
257  this->init(grid);
258 }
259 
260 template<typename GridT, typename InterruptT>
261 inline void
263 {
264  if (!grid.hasUniformVoxels()) {
265  OPENVDB_THROW(RuntimeError,
266  "The transform must have uniform scale for the LevelSetMeasure to function");
267  }
268  if (grid.getGridClass() != GRID_LEVEL_SET) {
269  OPENVDB_THROW(RuntimeError,
270  "LevelSetMeasure only supports level sets;"
271  " try setting the grid class to \"level set\"");
272  }
273  if (grid.empty()) {
274  OPENVDB_THROW(RuntimeError,
275  "LevelSetMeasure does not support empty grids;");
276  }
277  mGrid = &grid;
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;
282 }
283 
284 template<typename GridT, typename InterruptT>
285 inline Real
287 {
288  if (mUpdateArea) {MeasureArea m(this);};
289  double area = mArea;
290  if (useWorldUnits) area *= math::Pow2(mDx);
291  return area;
292 }
293 
294 template<typename GridT, typename InterruptT>
295 inline Real
297 {
298  if (mUpdateArea) {MeasureArea m(this);};
299  double volume = mVolume;
300  if (useWorldUnits) volume *= math::Pow3(mDx) ;
301  return volume;
302 }
303 
304 template<typename GridT, typename InterruptT>
305 inline Real
307 {
308  if (mUpdateCurvature) {MeasureCurvatures m(this);};
309  return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
310 }
311 
312 template<typename GridT, typename InterruptT>
313 inline Real
315 {
316  if (mUpdateCurvature) {MeasureCurvatures m(this);};
317  return mTotGausCurvature;
318 }
319 
320 template<typename GridT, typename InterruptT>
321 inline int
323 {
324  const Real x = this->totGaussianCurvature(true) / (2.0*math::pi<Real>());
325  return int(math::Round( x ));
326 }
327 
328 ///////////////////////// PRIVATE METHODS //////////////////////
329 
330 template<typename GridT, typename InterruptT>
331 inline bool
333 {
334  if (util::wasInterrupted(mInterrupter)) {
336  return false;
337  }
338  return true;
339 }
340 
341 template<typename GridT, typename InterruptT>
342 inline void
343 LevelSetMeasure<GridT, InterruptT>::
344 MeasureArea::operator()(const LeafRange& range) const
345 {
346  using Vec3T = math::Vec3<ValueType>;
347  // computations are performed in index space where dV = 1
348  mParent->checkInterrupter();
349  const Real invDx = 1.0/mParent->mDx;
350  const DiracDelta<Real> DD(1.5);// dirac delta function is 3 voxel units wide
351  const size_t leafCount = mParent->mLeafs->leafCount();
352  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
353  Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation
354  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
355  const Real dd = DD(invDx * (*voxelIter));
356  if (dd > 0.0) {
357  mStencil.moveTo(voxelIter);
358  const Coord& p = mStencil.getCenterCoord();// in voxel units
359  const Vec3T g = mStencil.gradient();// in world units
360  sumA += dd*g.length();// \delta(\phi)*|\nabla\phi|
361  sumV += dd*(g[0]*Real(p[0]) + g[1]*Real(p[1]) + g[2]*Real(p[2]));// \delta(\phi)\vec{x}\cdot\nabla\phi
362  }
363  }
364  double* ptr = mParent->mBuffer.get() + leafIter.pos();
365  *ptr = sumA;
366  ptr += leafCount;
367  *ptr = sumV;
368  }
369 }
370 
371 template<typename GridT, typename InterruptT>
372 inline void
373 LevelSetMeasure<GridT, InterruptT>::
374 MeasureCurvatures::operator()(const LeafRange& range) const
375 {
376  using Vec3T = math::Vec3<ValueType>;
377  // computations are performed in index space where dV = 1
378  mParent->checkInterrupter();
379  const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
380  const DiracDelta<Real> DD(1.5);// dirac delta function is 3 voxel units wide
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;//reduce risk of catastrophic cancellation
385  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
386  const Real dd = DD(invDx * (*voxelIter));
387  if (dd > 0.0) {
388  mStencil.moveTo(voxelIter);
389  const Vec3T g = mStencil.gradient();
390  const Real dA = dd*g.length();// \delta(\phi)*\delta(\phi)
391  mStencil.curvatures(mean, gauss);
392  sumM += dA*mean*dx;// \delta(\phi)*\delta(\phi)*MeanCurvature
393  sumG += dA*gauss*dx2;// \delta(\phi)*\delta(\phi)*GaussCurvature
394  }
395  }
396  double* ptr = mParent->mBuffer.get() + leafIter.pos();
397  *ptr = sumM;
398  ptr += leafCount;
399  *ptr = sumG;
400  }
401 }
402 
403 ////////////////////////////////////////
404 
405 //{
406 /// @cond OPENVDB_DOCS_INTERNAL
407 
408 template<class GridT>
409 inline
411 doLevelSetArea(const GridT& grid, bool useWorldUnits)
412 {
413  LevelSetMeasure<GridT> m(grid);
414  return m.area(useWorldUnits);
415 }
416 
417 template<class GridT>
418 inline
420 doLevelSetArea(const GridT&, bool)
421 {
422  OPENVDB_THROW(TypeError,
423  "level set area is supported only for scalar, floating-point grids");
424 }
425 
426 /// @endcond
427 //}
428 
429 template<class GridT>
430 inline Real
431 levelSetArea(const GridT& grid, bool useWorldUnits)
432 {
433  return doLevelSetArea<GridT>(grid, useWorldUnits);
434 }
435 
436 ////////////////////////////////////////
437 
438 //{
439 /// @cond OPENVDB_DOCS_INTERNAL
440 
441 template<class GridT>
442 inline
444 doLevelSetVolume(const GridT& grid, bool useWorldUnits)
445 {
446  LevelSetMeasure<GridT> m(grid);
447  return m.volume(useWorldUnits);
448 }
449 
450 template<class GridT>
451 inline
453 doLevelSetVolume(const GridT&, bool)
454 {
455  OPENVDB_THROW(TypeError,
456  "level set volume is supported only for scalar, floating-point grids");
457 }
458 
459 /// @endcond
460 //}
461 
462 template<class GridT>
463 inline Real
464 levelSetVolume(const GridT& grid, bool useWorldUnits)
465 {
466  return doLevelSetVolume<GridT>(grid, useWorldUnits);
467 }
468 
469 ////////////////////////////////////////
470 
471 //{
472 /// @cond OPENVDB_DOCS_INTERNAL
473 
474 template<class GridT>
475 inline
477 doLevelSetEulerCharacteristic(const GridT& grid)
478 {
479  LevelSetMeasure<GridT> m(grid);
480  return m.eulerCharacteristic();
481 }
482 
483 template<class GridT>
484 inline
486 doLevelSetEulerCharacteristic(const GridT&)
487 {
488  OPENVDB_THROW(TypeError,
489  "level set euler characteristic is supported only for scalar, floating-point grids");
490 }
491 
492 /// @endcond
493 //}
494 
495 
496 template<class GridT>
497 inline int
498 levelSetEulerCharacteristic(const GridT& grid)
499 {
500  return doLevelSetEulerCharacteristic(grid);
501 }
502 
503 ////////////////////////////////////////
504 
505 //{
506 /// @cond OPENVDB_DOCS_INTERNAL
507 
508 template<class GridT>
509 inline
511 doLevelSetEuler(const GridT& grid)
512 {
513  LevelSetMeasure<GridT> m(grid);
514  return m.eulerCharacteristics();
515 
516 }
517 
518 template<class GridT>
519 inline
521 doLevelSetGenus(const GridT& grid)
522 {
523  LevelSetMeasure<GridT> m(grid);
524  return m.genus();
525 }
526 
527 template<class GridT>
528 inline
530 doLevelSetGenus(const GridT&)
531 {
532  OPENVDB_THROW(TypeError,
533  "level set genus is supported only for scalar, floating-point grids");
534 }
535 
536 /// @endcond
537 //}
538 
539 template<class GridT>
540 inline int
541 levelSetGenus(const GridT& grid)
542 {
543  return doLevelSetGenus(grid);
544 }
545 
546 } // namespace tools
547 } // namespace OPENVDB_VERSION_NAME
548 } // namespace openvdb
549 
550 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
Real totMeanCurvature(bool useWorldUnits=true)
Compute the total mean curvature of the level set surface.
void parallel_for(int64_t start, int64_t end, std::function< void(int64_t index)> &&task, parallel_options opt=parallel_options(0, Split_Y, 1))
Definition: parallel.h:127
SYS_API double cos(double x)
Smeared-out and continuous Dirac Delta function.
Real levelSetArea(const GridType &grid, bool useWorldSpace=true)
Return the surface area of a narrow-band level set.
GLint first
Definition: glcorearb.h:404
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
int levelSetGenus(const GridType &grid)
Return the genus of a narrow-band level set surface.
GLboolean GLboolean g
Definition: glcorearb.h:1221
int eulerCharacteristic()
Compute the Euler characteristic of the level set surface.
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:178
GLint GLenum GLint x
Definition: glcorearb.h:408
Real avgMeanCurvature(bool useWorldUnits=true)
Compute the average mean curvature of the level set surface.
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
Real totGaussianCurvature(bool useWorldUnits=true)
Compute the total gaussian curvature of the level set surface.
Real levelSetVolume(const GridType &grid, bool useWorldSpace=true)
Return the volume of a narrow-band level set surface.
float Round(float x)
Return x rounded to the nearest integer.
Definition: Math.h:822
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
typename tree::LeafManager< const TreeType > ManagerType
This class manages a linear array of pointers to a given tree's leaf nodes, as well as optional auxil...
Definition: LeafManager.h:84
Type Pow3(Type x)
Return x3.
Definition: Math.h:555
GLfloat GLfloat p
Definition: glew.h:16656
typedef int(WINAPI *PFNWGLRELEASEPBUFFERDCARBPROC)(HPBUFFERARB hPbuffer
Real volume(bool useWorldUnits=true)
Compute the volume of the level set surface.
const void * ptr(const T *p)
Definition: format.h:3603
int levelSetEulerCharacteristic(const GridType &grid)
Return the Euler Characteristics of a narrow-band level set surface (possibly disconnected).
Real avgGaussianCurvature(bool useWorldUnits=true)
Compute the average gaussian curvature of the level set surface.
void init(const GridType &grid)
Re-initialize using the specified grid.
int genus()
Compute the genus of the level set surface.
const GLdouble * m
Definition: glew.h:9166
GLsizei const GLfloat * value
Definition: glcorearb.h:823
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
constexpr T pi()
Pi constant taken from Boost to match old behaviour.
Definition: Math.h:118
GLenum GLint * range
Definition: glcorearb.h:1924
GLintptr offset
Definition: glcorearb.h:664
Real area(bool useWorldUnits=true)
Compute the surface area of the level set.
A LeafManager manages a linear array of pointers to a given tree's leaf nodes, as well as optional au...
LevelSetMeasure(const GridType &grid, InterruptT *interrupt=nullptr)
Main constructor from a grid.
bool wasInterrupted(T *i, int percent=-1)
Multi-threaded computation of surface area, volume and average mean-curvature for narrow band level s...
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:114
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
GA_API const UT_StringHolder area