OpenVDB  3.1.0
LevelSetMeasure.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2015 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 //
30 //
34 
35 #ifndef OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
36 #define OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
37 
38 #include <tbb/parallel_for.h>
39 #include <tbb/parallel_sort.h>
40 #include <boost/bind.hpp>
41 #include <boost/function.hpp>
42 #include <boost/type_traits/is_floating_point.hpp>
43 #include <boost/utility/enable_if.hpp>
44 #include <boost/math/constants/constants.hpp>//for Pi
45 #include <openvdb/math/Math.h>
46 #include <openvdb/Types.h>
47 #include <openvdb/Grid.h>
48 #include <openvdb/tree/LeafManager.h>
49 #include <openvdb/tree/ValueAccessor.h>
50 #include <openvdb/math/FiniteDifference.h>
51 #include <openvdb/math/Operators.h>
52 #include <openvdb/util/NullInterrupter.h>
53 
54 namespace openvdb {
56 namespace OPENVDB_VERSION_NAME {
57 namespace tools {
58 
67 template<class GridType>
68 inline Real
69 levelSetArea(const GridType& grid, bool useWorldSpace = true);
70 
79 template<class GridType>
80 inline Real
81 levelSetVolume(const GridType& grid, bool useWorldSpace = true);
82 
93 template<class GridType>
94 inline void
95 levelSetMeasure(const GridType& grid, Real& area, Real& volume, bool useWorldSpace = true);
96 
108 template<class GridType>
109 inline void
110 levelSetMeasure(const GridType& grid, Real& area, Real& volume, Real& avgCurvature,
111  bool useWorldSpace = true);
112 
114 template<typename RealT>
116 {
117 public:
118  DiracDelta(RealT eps) : mC(0.5/eps), mD(2*boost::math::constants::pi<RealT>()*mC), mE(eps) {}
119  inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
120 private:
121  const RealT mC, mD, mE;
122 };
123 
124 
132 template<typename GridT,
133  typename InterruptT = util::NullInterrupter>
135 {
136 public:
137  typedef GridT GridType;
138  typedef typename GridType::TreeType TreeType;
139  typedef typename TreeType::ValueType ValueType;
141  BOOST_STATIC_ASSERT(boost::is_floating_point<ValueType>::value);
142 
147  LevelSetMeasure(const GridType& grid, InterruptT* interrupt = NULL);
148 
149  LevelSetMeasure(ManagerType& leafs, Real Dx, InterruptT* interrupt);
150 
152  void reinit(const GridType& grid);
153 
155  void reinit(ManagerType& leafs, Real dx);
156 
158  virtual ~LevelSetMeasure() {}
159 
161  int getGrainSize() const { return mGrainSize; }
162 
165  void setGrainSize(int grainsize) { mGrainSize = grainsize; }
166 
172  void measure(Real& area, Real& volume, bool useWorldUnits = true);
173 
179  void measure(Real& area, Real& volume, Real& avgMeanCurvature, bool useWorldUnits = true);
180 
181 private:
182  // disallow copy construction and copy by assignment!
183  LevelSetMeasure(const LevelSetMeasure&);// not implemented
184  LevelSetMeasure& operator=(const LevelSetMeasure&);// not implemented
185 
186  const TreeType* mTree;
187  ManagerType* mLeafs;
188  InterruptT* mInterrupter;
189  double mDx;
190  double* mArray;
191  int mGrainSize;
192 
193  // @brief Return false if the process was interrupted
194  bool checkInterrupter();
195 
196  typedef typename TreeType::LeafNodeType LeafT;
197  typedef typename LeafT::ValueOnCIter VoxelCIterT;
198  typedef typename ManagerType::LeafRange LeafRange;
199  typedef typename LeafRange::Iterator LeafIterT;
200 
201  struct Measure2
202  {
203  Measure2(LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
204  {
205  if (parent->mGrainSize>0) {
206  tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
207  } else {
208  (*this)(parent->mLeafs->leafRange());
209  }
210  }
211  Measure2(const Measure2& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
212  void operator()(const LeafRange& range) const;
213  LevelSetMeasure* mParent;
214  typename GridT::ConstAccessor mAcc;
215  };
216  struct Measure3
217  {
218  Measure3(LevelSetMeasure* parent) : mParent(parent), mAcc(*mParent->mTree)
219  {
220  if (parent->mGrainSize>0) {
221  tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
222  } else {
223  (*this)(parent->mLeafs->leafRange());
224  }
225  }
226  Measure3(const Measure3& other) : mParent(other.mParent), mAcc(*mParent->mTree) {}
227  void operator()(const LeafRange& range) const;
228  LevelSetMeasure* mParent;
229  typename GridT::ConstAccessor mAcc;
230  };
231  inline double reduce(double* first, double scale)
232  {
233  double* last = first + mLeafs->leafCount();
234  tbb::parallel_sort(first, last);//reduces catastrophic cancellation
235  Real sum = 0.0;
236  while(first != last) sum += *first++;
237  return scale * sum;
238  }
239 
240 }; // end of LevelSetMeasure class
241 
242 
243 template<typename GridT, typename InterruptT>
244 inline
246  : mTree(&(grid.tree()))
247  , mLeafs(NULL)
248  , mInterrupter(interrupt)
249  , mDx(grid.voxelSize()[0])
250  , mArray(NULL)
251  , mGrainSize(1)
252 {
253  if (!grid.hasUniformVoxels()) {
255  "The transform must have uniform scale for the LevelSetMeasure to function");
256  }
257  if (grid.getGridClass() != GRID_LEVEL_SET) {
259  "LevelSetMeasure only supports level sets;"
260  " try setting the grid class to \"level set\"");
261  }
262 }
263 
264 
265 template<typename GridT, typename InterruptT>
266 inline
268  ManagerType& leafs, Real dx, InterruptT* interrupt)
269  : mTree(&(leafs.tree()))
270  , mLeafs(&leafs)
271  , mInterrupter(interrupt)
272  , mDx(dx)
273  , mArray(NULL)
274  , mGrainSize(1)
275 {
276 }
277 
278 template<typename GridT, typename InterruptT>
279 inline void
281 {
282  if (!grid.hasUniformVoxels()) {
284  "The transform must have uniform scale for the LevelSetMeasure to function");
285  }
286  if (grid.getGridClass() != GRID_LEVEL_SET) {
288  "LevelSetMeasure only supports level sets;"
289  " try setting the grid class to \"level set\"");
290  }
291  mTree = &(grid.tree());
292  mLeafs = NULL;
293  mDx = grid.voxelSize()[0];
294 }
295 
296 
297 template<typename GridT, typename InterruptT>
298 inline void
300 {
301  mLeafs = &leafs;
302  mTree = &(leafs.tree());
303  mDx = dx;
304 }
305 
307 
308 
309 template<typename GridT, typename InterruptT>
310 inline void
311 LevelSetMeasure<GridT, InterruptT>::measure(Real& area, Real& volume, bool useWorldUnits)
312 {
313  if (mInterrupter) mInterrupter->start("Measuring level set");
314 
315 
316  const bool newLeafs = mLeafs == NULL;
317  if (newLeafs) mLeafs = new ManagerType(*mTree);
318  const size_t leafCount = mLeafs->leafCount();
319  if (leafCount == 0) {
320  area = volume = 0;
321  return;
322  }
323  mArray = new double[2*leafCount];
324 
325  Measure2 m(this);
326 
327  const double dx = useWorldUnits ? mDx : 1.0;
328  area = this->reduce(mArray, math::Pow2(dx));
329  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
330 
331  if (newLeafs) {
332  delete mLeafs;
333  mLeafs = NULL;
334  }
335  delete [] mArray;
336 
337  if (mInterrupter) mInterrupter->end();
338 }
339 
340 
341 template<typename GridT, typename InterruptT>
342 inline void
344  Real& avgMeanCurvature,
345  bool useWorldUnits)
346 {
347  if (mInterrupter) mInterrupter->start("Measuring level set");
348 
349  const bool newLeafs = mLeafs == NULL;
350  if (newLeafs) mLeafs = new ManagerType(*mTree);
351  const size_t leafCount = mLeafs->leafCount();
352  if (leafCount == 0) {
353  area = volume = avgMeanCurvature = 0;
354  return;
355  }
356  mArray = new double[3*leafCount];
357 
358  Measure3 m(this);
359 
360  const double dx = useWorldUnits ? mDx : 1.0;
361  area = this->reduce(mArray, math::Pow2(dx));
362  volume = this->reduce(mArray + leafCount, math::Pow3(dx) / 3.0);
363  avgMeanCurvature = this->reduce(mArray + 2*leafCount, dx/area);
364 
365  if (newLeafs) {
366  delete mLeafs;
367  mLeafs = NULL;
368  }
369  delete [] mArray;
370 
371  if (mInterrupter) mInterrupter->end();
372 }
373 
374 
376 
377 
378 template<typename GridT, typename InterruptT>
379 inline bool
381 {
382  if (util::wasInterrupted(mInterrupter)) {
383  tbb::task::self().cancel_group_execution();
384  return false;
385  }
386  return true;
387 }
388 
389 template<typename GridT, typename InterruptT>
390 inline void
392 Measure2::operator()(const LeafRange& range) const
393 {
394  typedef math::Vec3<ValueType> Vec3T;
395  typedef math::ISGradient<math::CD_2ND> Grad;
396  mParent->checkInterrupter();
397  const Real invDx = 1.0/mParent->mDx;
398  const DiracDelta<Real> DD(1.5);
399  const size_t leafCount = mParent->mLeafs->leafCount();
400  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
401  Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation
402  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
403  const Real dd = DD(invDx * (*voxelIter));
404  if (dd > 0.0) {
405  const Coord p = voxelIter.getCoord();
406  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
407  sumA += dd * g.dot(g);
408  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
409  }
410  }
411  double* v = mParent->mArray + leafIter.pos();
412  *v = sumA;
413  v += leafCount;
414  *v = sumV;
415  }
416 }
417 
418 template<typename GridT, typename InterruptT>
419 inline void
421 Measure3::operator()(const LeafRange& range) const
422 {
423  typedef math::Vec3<ValueType> Vec3T;
424  typedef math::ISGradient<math::CD_2ND> Grad;
426  mParent->checkInterrupter();
427  const Real invDx = 1.0/mParent->mDx;
428  const DiracDelta<Real> DD(1.5);
429  ValueType alpha, beta;
430  const size_t leafCount = mParent->mLeafs->leafCount();
431  for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
432  Real sumA = 0, sumV = 0, sumC = 0;//reduce risk of catastrophic cancellation
433  for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
434  const Real dd = DD(invDx * (*voxelIter));
435  if (dd > 0.0) {
436  const Coord p = voxelIter.getCoord();
437  const Vec3T g = invDx*Grad::result(mAcc, p);//voxel units
438  const Real dA = dd * g.dot(g);
439  sumA += dA;
440  sumV += dd * (g[0]*p[0]+g[1]*p[1]+g[2]*p[2]);
441  Curv::result(mAcc, p, alpha, beta);
442  sumC += dA * alpha/(2*math::Pow2(beta))*invDx;
443  }
444  }
445  double* v = mParent->mArray + leafIter.pos();
446  *v = sumA;
447  v += leafCount;
448  *v = sumV;
449  v += leafCount;
450  *v = sumC;
451  }
452 }
453 
455 
456 template<class GridT>
457 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
458 doLevelSetArea(const GridT& grid, bool useWorldSpace)
459 {
460  Real area, volume;
461  LevelSetMeasure<GridT> m(grid);
462  m.measure(area, volume, useWorldSpace);
463  return area;
464 }
465 
466 template<class GridT>
467 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
468 doLevelSetArea(const GridT&, bool)
469 {
471  "level set area is supported only for scalar, floating-point grids");
472 }
473 
474 template<class GridT>
475 inline Real
476 levelSetArea(const GridT& grid, bool useWorldSpace)
477 {
478  return doLevelSetArea<GridT>(grid, useWorldSpace);
479 }
480 
482 
483 template<class GridT>
484 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
485 doLevelSetVolume(const GridT& grid, bool useWorldSpace)
486 {
487  Real area, volume;
488  LevelSetMeasure<GridT> m(grid);
489  m.measure(area, volume, useWorldSpace);
490  return volume;
491 }
492 
493 template<class GridT>
494 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType>, Real>::type
495 doLevelSetVolume(const GridT&, bool)
496 {
498  "level set volume is supported only for scalar, floating-point grids");
499 }
500 
501 template<class GridT>
502 inline Real
503 levelSetVolume(const GridT& grid, bool useWorldSpace)
504 {
505  return doLevelSetVolume<GridT>(grid, useWorldSpace);
506 }
507 
509 
510 template<class GridT>
511 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
512 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
513 {
514  LevelSetMeasure<GridT> m(grid);
515  m.measure(area, volume, useWorldSpace);
516 }
517 
518 template<class GridT>
519 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
520 doLevelSetMeasure(const GridT&, Real&, Real&, bool)
521 {
523  "level set measure is supported only for scalar, floating-point grids");
524 }
525 
526 template<class GridT>
527 inline void
528 levelSetMeasure(const GridT& grid, Real& area, Real& volume, bool useWorldSpace)
529 {
530  doLevelSetMeasure<GridT>(grid, area, volume, useWorldSpace);
531 }
532 
534 
535 template<class GridT>
536 inline typename boost::enable_if<boost::is_floating_point<typename GridT::ValueType> >::type
537 doLevelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature,
538  bool useWorldSpace)
539 {
540  LevelSetMeasure<GridT> m(grid);
541  m.measure(area, volume, avgCurvature, useWorldSpace);
542 }
543 
544 template<class GridT>
545 inline typename boost::disable_if<boost::is_floating_point<typename GridT::ValueType> >::type
546 doLevelSetMeasure(const GridT&, Real&, Real&, Real&, bool)
547 {
549  "level set measure is supported only for scalar, floating-point grids");
550 }
551 
552 template<class GridT>
553 inline void
554 levelSetMeasure(const GridT& grid, Real& area, Real& volume, Real& avgCurvature, bool useWorldSpace)
555 {
556  doLevelSetMeasure<GridT>(grid, area, volume, avgCurvature, useWorldSpace);
557 }
558 
559 } // namespace tools
560 } // namespace OPENVDB_VERSION_NAME
561 } // namespace openvdb
562 
563 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
564 
565 // Copyright (c) 2012-2015 DreamWorks Animation LLC
566 // All rights reserved. This software is distributed under the
567 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:76
void reinit(const GridType &grid)
Re-initialize using the specified grid.
Definition: LevelSetMeasure.h:280
void setGrainSize(int grainsize)
Set the grain-size used for multi-threading.
Definition: LevelSetMeasure.h:165
TreeType::ValueType ValueType
Definition: LevelSetMeasure.h:139
int32_t Abs(int32_t i)
Return the absolute value of the given quantity.
Definition: Math.h:293
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 ...
Definition: LevelSetMeasure.h:311
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
LeafRange leafRange(size_t grainsize=1) const
Return a TBB-compatible LeafRange.
Definition: LeafManager.h:350
Multi-threaded computation of surface area, volume and average mean-curvature for narrow band level s...
Definition: LevelSetMeasure.h:134
This class manages a linear array of pointers to a given tree&#39;s leaf nodes, as well as optional auxil...
Definition: LeafManager.h:114
void levelSetMeasure(const GridT &grid, Real &area, Real &volume, Real &avgCurvature, bool useWorldSpace)
Definition: LevelSetMeasure.h:554
LevelSetMeasure(const GridType &grid, InterruptT *interrupt=NULL)
Main constructor from a grid.
Definition: LevelSetMeasure.h:245
Smeared-out and continuous Dirac Delta function.
Definition: LevelSetMeasure.h:115
RealT operator()(RealT phi) const
Definition: LevelSetMeasure.h:119
Iterator begin() const
Definition: LeafManager.h:181
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Real levelSetArea(const GridType &grid, bool useWorldSpace=true)
Return the surface area of a narrow-band level set.
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
Real levelSetVolume(const GridType &grid, bool useWorldSpace=true)
Return the volume of a narrow-band level set surface.
const TreeType & tree() const
Return a const reference to tree associated with this manager.
Definition: LeafManager.h:307
Definition: Mat.h:146
Definition: Exceptions.h:86
void levelSetMeasure(const GridType &grid, Real &area, Real &volume, bool useWorldSpace=true)
Compute the surface area and volume of a narrow-band level set.
Definition: Exceptions.h:39
GridType::TreeType TreeType
Definition: LevelSetMeasure.h:138
size_t leafCount() const
Return the number of leaf nodes.
Definition: LeafManager.h:304
Definition: Types.h:206
Type Pow2(Type x)
Return .
Definition: Math.h:514
Gradient operators defined in index space of various orders.
Definition: Operators.h:122
GridT GridType
Definition: LevelSetMeasure.h:137
boost::enable_if< boost::is_floating_point< typename GridT::ValueType >, Real >::type doLevelSetArea(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:458
boost::enable_if< boost::is_floating_point< typename GridT::ValueType >, Real >::type doLevelSetVolume(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:485
double Real
Definition: Types.h:64
virtual ~LevelSetMeasure()
Destructor.
Definition: LevelSetMeasure.h:158
Type Pow3(Type x)
Return .
Definition: Math.h:518
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
int getGrainSize() const
Definition: LevelSetMeasure.h:161
MatType scale(const Vec3< typename MatType::value_type > &s)
Return a matrix that scales by s.
Definition: Mat.h:594
tree::LeafManager< const TreeType > ManagerType
Definition: LevelSetMeasure.h:140
Definition: Exceptions.h:87
boost::enable_if< boost::is_floating_point< typename GridT::ValueType > >::type doLevelSetMeasure(const GridT &grid, Real &area, Real &volume, bool useWorldSpace)
Definition: LevelSetMeasure.h:512
Real levelSetArea(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:476
DiracDelta(RealT eps)
Definition: LevelSetMeasure.h:118
Real levelSetVolume(const GridT &grid, bool useWorldSpace)
Definition: LevelSetMeasure.h:503
Compute the mean curvature in index space.
Definition: Operators.h:555