HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevelSetMeasure.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2012-2018 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
29 ///////////////////////////////////////////////////////////////////////////
30 //
31 /// @author Ken Museth
32 ///
33 /// @file LevelSetMeasure.h
34 
35 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
37 
38 #include <openvdb/math/Math.h>
39 #include <openvdb/Types.h>
40 #include <openvdb/Grid.h>
44 #include <openvdb/math/Operators.h>
46 #include <hboost/math/constants/constants.hpp>//for Pi
47 #include <tbb/parallel_for.h>
48 #include <tbb/parallel_sort.h>
49 #include <type_traits>
50 
51 
52 namespace openvdb {
54 namespace OPENVDB_VERSION_NAME {
55 namespace tools {
56 
57 /// @brief Return the surface area of a narrow-band level set.
58 ///
59 /// @param grid a scalar, floating-point grid with one or more disjoint,
60 /// closed isosurfaces at the given @a isovalue
61 /// @param useWorldSpace if true the area is computed in
62 /// world space units, else in voxel units.
63 ///
64 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set.
65 template<class GridType>
66 inline Real
67 levelSetArea(const GridType& grid, bool useWorldSpace = true);
68 
69 /// @brief Return the volume of a narrow-band level set surface.
70 ///
71 /// @param grid a scalar, floating-point grid with one or more disjoint,
72 /// closed isosurfaces at the given @a isovalue
73 /// @param useWorldSpace if true the volume is computed in
74 /// world space units, else in voxel units.
75 ///
76 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set.
77 template<class GridType>
78 inline Real
79 levelSetVolume(const GridType& grid, bool useWorldSpace = true);
80 
81 /// @brief Compute the surface area and volume of a narrow-band level set.
82 ///
83 /// @param grid a scalar, floating-point grid with one or more disjoint,
84 /// closed isosurfaces at the given @a isovalue
85 /// @param area surface area of the level set
86 /// @param volume volume of the level set surface
87 /// @param useWorldSpace if true the area and volume are computed in
88 /// world space units, else in voxel units.
89 ///
90 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set.
91 template<class GridType>
92 inline void
93 levelSetMeasure(const GridType& grid, Real& area, Real& volume, bool useWorldSpace = true);
94 
95 /// @brief Compute the surface area and volume of a narrow-band level set.
96 ///
97 /// @param grid a scalar, floating-point grid with one or more disjoint,
98 /// closed isosurfaces at the given @a isovalue
99 /// @param area surface area of the level set
100 /// @param volume volume of the level set surface
101 /// @param avgCurvature average mean curvature of the level set surface
102 /// @param useWorldSpace if true the area, volume and curvature are computed in
103 /// world space units, else in voxel units.
104 ///
105 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set.
106 template<class GridType>
107 inline void
108 levelSetMeasure(const GridType& grid, Real& area, Real& volume, Real& avgCurvature,
109  bool useWorldSpace = true);
110 
111 /// @brief Smeared-out and continuous Dirac Delta function.
112 template<typename RealT>
114 {
115 public:
116  DiracDelta(RealT eps) : mC(0.5/eps), mD(2*hboost::math::constants::pi<RealT>()*mC), mE(eps) {}
117  inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
118 private:
119  const RealT mC, mD, mE;
120 };
121 
122 
123 /// @brief Multi-threaded computation of surface area, volume and
124 /// average mean-curvature for narrow band level sets.
125 ///
126 /// @details To reduce the risk of round-off errors (primarily due to
127 /// catastrophic cancellation) and guarantee determinism during
128 /// multi-threading this class is implemented using parallel_for, and
129 /// delayed reduction of a sorted list.
130 template<typename GridT, typename InterruptT = util::NullInterrupter>
132 {
133 public:
134  using GridType = GridT;
135  using TreeType = typename GridType::TreeType;
136  using ValueType = typename TreeType::ValueType;
139  "level set measure is supported only for scalar, floating-point grids");
140 
141  /// @brief Main constructor from a grid
142  /// @param grid The level set to be measured.
143  /// @param interrupt Optional interrupter.
144  /// @throw RuntimeError if the grid is not a level set.
145  LevelSetMeasure(const GridType& grid, InterruptT* interrupt = nullptr);
146 
147  LevelSetMeasure(ManagerType& leafs, Real Dx, InterruptT* interrupt);
148 
149  /// @brief Re-initialize using the specified grid.
150  void reinit(const GridType& grid);
151 
152  /// @brief Re-initialize using the specified LeafManager and voxelSize.
153  void reinit(ManagerType& leafs, Real dx);
154 
155  /// @brief Destructor
156  virtual ~LevelSetMeasure() {}
157 
158  /// @return the grain-size used for multi-threading
159  int getGrainSize() const { return mGrainSize; }
160 
161  /// @brief Set the grain-size used for multi-threading.
162  /// @note A grain size of 0 or less disables multi-threading!
163  void setGrainSize(int grainsize) { mGrainSize = grainsize; }
164 
165  /// @brief Compute the surface area and volume of the level
166  /// set. Use the last argument to specify the result in world or
167  /// voxel units.
168  /// @note This method is faster (about 3x) then the measure method
169  /// below that also computes the average mean-curvature.
170  void measure(Real& area, Real& volume, bool useWorldUnits = true);
171 
172  /// @brief Compute the surface area, volume, and average
173  /// mean-curvature of the level set. Use the last argument to
174  /// specify the result in world or voxel units.
175  /// @note This method is slower (about 3x) then the measure method
176  /// above that only computes the area and volume.
177  void measure(Real& area, Real& volume, Real& avgMeanCurvature, bool useWorldUnits = true);
178 
179 private:
180  // disallow copy construction and copy by assignment!
181  LevelSetMeasure(const LevelSetMeasure&);// not implemented
182  LevelSetMeasure& operator=(const LevelSetMeasure&);// not implemented
183 
184  const TreeType* mTree;
185  ManagerType* mLeafs;
186  InterruptT* mInterrupter;
187  double mDx;
188  double* mArray;
189  int mGrainSize;
190 
191  // @brief Return false if the process was interrupted
192  bool checkInterrupter();
193 
194  using LeafT = typename TreeType::LeafNodeType;
195  using VoxelCIterT = typename LeafT::ValueOnCIter;
196  using LeafRange = typename ManagerType::LeafRange;
197  using LeafIterT = typename LeafRange::Iterator;
198 
199  struct Measure2
200  {
201  Measure2(LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
202  {
203  if (parent->mGrainSize>0) {
204  tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
205  } else {
206  (*this)(parent->mLeafs->leafRange());
207  }
208  }
209  Measure2(const Measure2& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
210  void operator()(const LeafRange& range) const;
211  LevelSetMeasure* mParent;
212  typename GridT::ConstAccessor mAcc;
213  };
214  struct Measure3
215  {
216  Measure3(LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
217  {
218  if (parent->mGrainSize>0) {
219  tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
220  } else {
221  (*this)(parent->mLeafs->leafRange());
222  }
223  }
224  Measure3(const Measure3& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
225  void operator()(const LeafRange& range) const;
226  LevelSetMeasure* mParent;
227  typename GridT::ConstAccessor mAcc;
228  };
229  inline double reduce(double* first, double scale)
230  {
231  double* last = first + mLeafs->leafCount();
232  tbb::parallel_sort(first, last);//reduces catastrophic cancellation
233  Real sum = 0.0;
234  while(first != last) sum += *first++;
235  return scale * sum;
236  }
237 
238 }; // end of LevelSetMeasure class
239 
240 
241 template<typename GridT, typename InterruptT>
242 inline
244  : mTree(&(grid.tree()))
245  , mLeafs(nullptr)
246  , mInterrupter(interrupt)
247  , mDx(grid.voxelSize()[0])
248  , mArray(nullptr)
249  , mGrainSize(1)
250 {
251  if (!grid.hasUniformVoxels()) {
252  OPENVDB_THROW(RuntimeError,
253  "The transform must have uniform scale for the LevelSetMeasure to function");
254  }
255  if (grid.getGridClass() != GRID_LEVEL_SET) {
256  OPENVDB_THROW(RuntimeError,
257  "LevelSetMeasure only supports level sets;"
258  " try setting the grid class to \"level set\"");
259  }
260 }
261 
262 
263 template<typename GridT, typename InterruptT>
264 inline
266  ManagerType& leafs, Real dx, InterruptT* interrupt)
267  : mTree(&(leafs.tree()))
268  , mLeafs(&leafs)
269  , mInterrupter(interrupt)
270  , mDx(dx)
271  , mArray(nullptr)
272  , mGrainSize(1)
273 {
274 }
275 
276 template<typename GridT, typename InterruptT>
277 inline void
279 {
280  if (!grid.hasUniformVoxels()) {
281  OPENVDB_THROW(RuntimeError,
282  "The transform must have uniform scale for the LevelSetMeasure to function");
283  }
284  if (grid.getGridClass() != GRID_LEVEL_SET) {
285  OPENVDB_THROW(RuntimeError,
286  "LevelSetMeasure only supports level sets;"
287  " try setting the grid class to \"level set\"");
288  }
289  mTree = &(grid.tree());
290  mLeafs = nullptr;
291  mDx = grid.voxelSize()[0];
292 }
293 
294 
295 template<typename GridT, typename InterruptT>
296 inline void
298 {
299  mLeafs = &leafs;
300  mTree = &(leafs.tree());
301  mDx = dx;
302 }
303 
304 ////////////////////////////////////////
305 
306 
307 template<typename GridT, typename InterruptT>
308 inline void
310 {
311  if (mInterrupter) mInterrupter->start("Measuring level set");
312 
313 
314  const bool newLeafs = mLeafs == nullptr;
315  if (newLeafs) mLeafs = new ManagerType(*mTree);
316  const size_t leafCount = mLeafs->leafCount();
317  if (leafCount == 0) {
318  area = volume = 0;
319  return;
320  }
321  mArray = new double[2*leafCount];
322 
323  Measure2 m(this);
324 
325  const double dx = useWorldUnits ? mDx : 1.0;
326  area = this->reduce(mArray, math::Pow2(dx));
327  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
328 
329  if (newLeafs) {
330  delete mLeafs;
331  mLeafs = nullptr;
332  }
333  delete [] mArray;
334 
335  if (mInterrupter) mInterrupter->end();
336 }
337 
338 
339 template<typename GridT, typename InterruptT>
340 inline void
342  Real& avgMeanCurvature,
343  bool useWorldUnits)
344 {
345  if (mInterrupter) mInterrupter->start("Measuring level set");
346 
347  const bool newLeafs = mLeafs == nullptr;
348  if (newLeafs) mLeafs = new ManagerType(*mTree);
349  const size_t leafCount = mLeafs->leafCount();
350  if (leafCount == 0) {
351  area = volume = avgMeanCurvature = 0;
352  return;
353  }
354  mArray = new double[3*leafCount];
355 
356  Measure3 m(this);
357 
358  const double dx = useWorldUnits ? mDx : 1.0;
359  area = this->reduce(mArray, math::Pow2(dx));
360  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
361  avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
362 
363  if (newLeafs) {
364  delete mLeafs;
365  mLeafs = nullptr;
366  }
367  delete [] mArray;
368 
369  if (mInterrupter) mInterrupter->end();
370 }
371 
372 
373 ///////////////////////// PRIVATE METHODS //////////////////////
374 
375 
376 template<typename GridT, typename InterruptT>
377 inline bool
379 {
380  if (util::wasInterrupted(mInterrupter)) {
381  tbb::task::self().cancel_group_execution();
382  return false;
383  }
384  return true;
385 }
386 
387 template<typename GridT, typename InterruptT>
388 inline void
389 LevelSetMeasure<GridT, InterruptT>::
390 Measure2::operator()(const LeafRange& range) const
391 {
392  using Vec3T = math::Vec3<ValueType>;
393  using Grad = math::ISGradient<math::CD_2ND>;
394  mParent->checkInterrupter();
395  const Real invDx = 1.0/mParent->mDx;
396  const DiracDelta<Real> DD(1.5);
397  const size_t leafCount = mParent->mLeafs->leafCount();
398  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
399  Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation
400  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
401  const Real dd = DD(invDx * (*voxelIter));
402  if (dd > 0.0) {
403  const Coord p = voxelIter.getCoord();
404  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
405  sumA += dd * g.dot(g);
406  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
407  }
408  }
409  double* v = mParent->mArray + leafIter.pos();
410  *v = sumA;
411  v += leafCount;
412  *v = sumV;
413  }
414 }
415 
416 template<typename GridT, typename InterruptT>
417 inline void
418 LevelSetMeasure<GridT, InterruptT>::
419 Measure3::operator()(const LeafRange& range) const
420 {
421  using Vec3T = math::Vec3<ValueType>;
422  using Grad = math::ISGradient<math::CD_2ND>;
423  using Curv = math::ISMeanCurvature<math::CD_SECOND, math::CD_2ND>;
424  mParent->checkInterrupter();
425  const Real invDx = 1.0/mParent->mDx;
426  const DiracDelta<Real> DD(1.5);
427  ValueType alpha, beta;
428  const size_t leafCount = mParent->mLeafs->leafCount();
429  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
430  Real sumA = 0, sumV = 0, sumC = 0;//reduce risk of catastrophic cancellation
431  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
432  const Real dd = DD(invDx * (*voxelIter));
433  if (dd > 0.0) {
434  const Coord p = voxelIter.getCoord();
435  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
436  const Real dA = dd * g.dot(g);
437  sumA += dA;
438  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
439  Curv::result(mAcc, p, alpha, beta);
440  sumC += dA * alpha/(2*math::Pow2(beta))*invDx;
441  }
442  }
443  double* v = mParent->mArray + leafIter.pos();
444  *v = sumA;
445  v += leafCount;
446  *v = sumV;
447  v += leafCount;
448  *v = sumC;
449  }
450 }
451 
452 
453 ////////////////////////////////////////
454 
455 
456 //{
457 /// @cond OPENVDB_LEVEL_SET_MEASURE_INTERNAL
458 
459 template<class GridT>
460 inline
462 doLevelSetArea(const GridT& grid, bool useWorldSpace)
463 {
464  Real area, volume;
465  LevelSetMeasure<GridT> m(grid);
466  m.measure(area, volume, useWorldSpace);
467  return area;
468 }
469 
470 template<class GridT>
471 inline
473 doLevelSetArea(const GridT&, bool)
474 {
475  OPENVDB_THROW(TypeError,
476  "level set area is supported only for scalar, floating-point grids");
477 }
478 
479 /// @endcond
480 //}
481 
482 
483 template<class GridT>
484 inline Real
485 levelSetArea(const GridT& grid, bool useWorldSpace)
486 {
487  return doLevelSetArea<GridT>(grid, useWorldSpace);
488 }
489 
490 
491 ////////////////////////////////////////
492 
493 
494 //{
495 /// @cond OPENVDB_LEVEL_SET_MEASURE_INTERNAL
496 
497 template<class GridT>
498 inline
500 doLevelSetVolume(const GridT& grid, bool useWorldSpace)
501 {
502  Real area, volume;
503  LevelSetMeasure<GridT> m(grid);
504  m.measure(area, volume, useWorldSpace);
505  return volume;
506 }
507 
508 template<class GridT>
509 inline
511 doLevelSetVolume(const GridT&, bool)
512 {
513  OPENVDB_THROW(TypeError,
514  "level set volume is supported only for scalar, floating-point grids");
515 }
516 
517 /// @endcond
518 //}
519 
520 
521 template<class GridT>
522 inline Real
523 levelSetVolume(const GridT& grid, bool useWorldSpace)
524 {
525  return doLevelSetVolume<GridT>(grid, useWorldSpace);
526 }
527 
528 
529 ////////////////////////////////////////
530 
531 
532 //{
533 /// @cond OPENVDB_LEVEL_SET_MEASURE_INTERNAL
534 
535 template<class GridT>
536 inline
538 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
539 {
540  LevelSetMeasure<GridT> m(grid);
541  m.measure(area, volume, useWorldSpace);
542 }
543 
544 template<class GridT>
545 inline
547 doLevelSetMeasure(const GridT&, Real&, Real&, bool)
548 {
549  OPENVDB_THROW(TypeError,
550  "level set measure is supported only for scalar, floating-point grids");
551 }
552 
553 /// @endcond
554 //}
555 
556 
557 template<class GridT>
558 inline void
559 levelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
560 {
561  doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
562 }
563 
564 
565 ////////////////////////////////////////
566 
567 
568 //{
569 /// @cond OPENVDB_LEVEL_SET_MEASURE_INTERNAL
570 
571 template<class GridT>
572 inline
574 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature,
575  bool useWorldSpace)
576 {
577  LevelSetMeasure<GridT> m(grid);
578  m.measure(area, volume, avgCurvature, useWorldSpace);
579 }
580 
581 template<class GridT>
582 inline
584 doLevelSetMeasure(const GridT&, Real&, Real&, Real&, bool)
585 {
586  OPENVDB_THROW(TypeError,
587  "level set measure is supported only for scalar, floating-point grids");
588 }
589 
590 /// @endcond
591 //}
592 
593 
594 template<class GridT>
595 inline void
596 levelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature, bool useWorldSpace)
597 {
598  doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
599 }
600 
601 } // namespace tools
602 } // namespace OPENVDB_VERSION_NAME
603 } // namespace openvdb
604 
605 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
606 
607 // Copyright (c) 2012-2018 DreamWorks Animation LLC
608 // All rights reserved. This software is distributed under the
609 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
GLint first
Definition: glcorearb.h:404
SYS_API double cos(double x)
Definition: SYS_FPUMath.h:69
GLenum GLint * range
Definition: glcorearb.h:1924
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.
Type Pow2(Type x)
Return x2.
Definition: Math.h:502
const GLdouble * v
Definition: glcorearb.h:836
GLboolean GLboolean g
Definition: glcorearb.h:1221
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
GA_API const UT_StringHolder scale
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
void reinit(const GridType &grid)
Re-initialize using the specified grid.
Coord Abs(const Coord &xyz)
Definition: Coord.h:513
Real levelSetVolume(const GridType &grid, bool useWorldSpace=true)
Return the volume of a narrow-band level set surface.
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:110
Type Pow3(Type x)
Return x3.
Definition: Math.h:506
GLfloat GLfloat GLfloat alpha
Definition: glcorearb.h:111
void levelSetMeasure(const GridType &grid, Real &area, Real &volume, bool useWorldSpace=true)
Compute the surface area and volume of a narrow-band level set.
GLsizei const GLfloat * value
Definition: glcorearb.h:823
void measure(Real &area, Real &volume, bool useWorldUnits=true)
Compute the surface area and volume of the level set. Use the last argument to specify the result in ...
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
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:135
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
GA_API const UT_StringHolder area