GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/tools/LevelSetMeasure.h
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 108 123 87.8%
Functions: 20 56 35.7%
Branches: 62 126 49.2%
Line Branch Exec Source
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"
13 #include "openvdb/tree/LeafManager.h"
14 #include "openvdb/tree/ValueAccessor.h"
15 #include "openvdb/math/Math.h"
16 #include "openvdb/math/FiniteDifference.h"
17 #include "openvdb/math/Operators.h"
18 #include "openvdb/math/Stencils.h"
19 #include "openvdb/util/NullInterrupter.h"
20 #include "openvdb/thread/Threading.h"
21 #include <openvdb/openvdb.h>
22
23 #include <tbb/parallel_for.h>
24 #include <tbb/parallel_sort.h>
25 #include <tbb/parallel_invoke.h>
26
27 #include <type_traits>
28
29 namespace openvdb {
30 OPENVDB_USE_VERSION_NAMESPACE
31 namespace OPENVDB_VERSION_NAME {
32 namespace tools {
33
34 /// @brief Return the surface area of a narrow-band level set.
35 ///
36 /// @param grid a scalar, floating-point grid with one or more disjoint,
37 /// closed level set surfaces
38 /// @param useWorldSpace if true the area is computed in
39 /// world space units, else in voxel units.
40 ///
41 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
42 template<class GridType>
43 Real
44 levelSetArea(const GridType& grid, bool useWorldSpace = true);
45
46 /// @brief Return the volume of a narrow-band level set surface.
47 ///
48 /// @param grid a scalar, floating-point grid with one or more disjoint,
49 /// closed level set surfaces
50 /// @param useWorldSpace if true the volume is computed in
51 /// world space units, else in voxel units.
52 ///
53 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
54 template<class GridType>
55 Real
56 levelSetVolume(const GridType& grid, bool useWorldSpace = true);
57
58 /// @brief Return the Euler Characteristics of a narrow-band level set surface (possibly disconnected).
59 ///
60 /// @param grid a scalar, floating-point grid with one or more disjoint,
61 /// closed level set surfaces
62 ///
63 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
64 template<class GridType>
65 int
66 levelSetEulerCharacteristic(const GridType& grid);
67
68 /// @brief Return the genus of a narrow-band level set surface.
69 ///
70 /// @param grid a scalar, floating-point grid with one or more disjoint,
71 /// closed level set surfaces
72 /// @warning The genus is only well defined for a single connected surface
73 ///
74 /// @throw TypeError if @a grid is not scalar or not floating-point or not a level set or empty.
75 template<class GridType>
76 int
77 levelSetGenus(const GridType& grid);
78
79 ////////////////////////////////////////////////////////////////////////////////////////
80
81 /// @brief Smeared-out and continuous Dirac Delta function.
82 template<typename RealT>
83 class DiracDelta
84 {
85 public:
86 // eps is the half-width of the dirac delta function in units of phi
87 DiracDelta(RealT eps) : mC(0.5/eps), mD(2*math::pi<RealT>()*mC), mE(eps) {}
88 // values of the dirac delta function are in units of one over the units of phi
89
4/8
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2441734 times.
✓ Branch 5 taken 2445978 times.
✓ Branch 6 taken 9640106 times.
✓ Branch 7 taken 9475644 times.
 24003462 inline RealT operator()(RealT phi) const { return math::Abs(phi) > mE ? 0 : mC*(1+cos(mD*phi)); }
90 private:
91 const RealT mC, mD, mE;
92 };// DiracDelta functor
93
94
95 /// @brief Multi-threaded computation of surface area, volume and
96 /// average mean-curvature for narrow band level sets.
97 ///
98 /// @details To reduce the risk of round-off errors (primarily due to
99 /// catastrophic cancellation) and guarantee determinism during
100 /// multi-threading this class is implemented using parallel_for, and
101 /// delayed reduction of a sorted list.
102 template<typename GridT, typename InterruptT = util::NullInterrupter>
103 class LevelSetMeasure
104 {
105 public:
106 using GridType = GridT;
107 using TreeType = typename GridType::TreeType;
108 using ValueType = typename TreeType::ValueType;
109 using ManagerType = typename tree::LeafManager<const TreeType>;
110
111 static_assert(std::is_floating_point<ValueType>::value,
112 "level set measure is supported only for scalar, floating-point grids");
113
114 /// @brief Main constructor from a grid
115 /// @param grid The level set to be measured.
116 /// @param interrupt Optional interrupter.
117 /// @throw RuntimeError if the grid is not a level set or if it's empty.
118 LevelSetMeasure(const GridType& grid, InterruptT* interrupt = nullptr);
119
120 /// @brief Re-initialize using the specified grid.
121 /// @param grid The level set to be measured.
122 /// @throw RuntimeError if the grid is not a level set or if it's empty.
123 void init(const GridType& grid);
124
125 /// @brief Destructor
126
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 64 virtual ~LevelSetMeasure() {}
127
128 /// @return the grain-size used for multi-threading
129 ✗ int getGrainSize() const { return mGrainSize; }
130
131 /// @brief Set the grain-size used for multi-threading.
132 /// @note A grain size of 0 or less disables multi-threading!
133 ✗ void setGrainSize(int grainsize) { mGrainSize = grainsize; }
134
135 /// @brief Compute the surface area of the level set.
136 /// @param useWorldUnits Specifies if the result is in world or voxel units.
137 /// @note Performs internal caching so only the initial call incurs actual computation.
138 Real area(bool useWorldUnits = true);
139
140 /// @brief Compute the volume of the level set surface.
141 /// @param useWorldUnits Specifies if the result is in world or voxel units.
142 /// @note Performs internal caching so only the initial call incurs actual computation.
143 Real volume(bool useWorldUnits = true);
144
145 /// @brief Compute the total mean curvature of the level set surface.
146 /// @param useWorldUnits Specifies if the result is in world or voxel units.
147 /// @note Performs internal caching so only the initial call incurs actual computation.
148 Real totMeanCurvature(bool useWorldUnits = true);
149
150 /// @brief Compute the total gaussian curvature of the level set surface.
151 /// @param useWorldUnits Specifies if the result is in world or voxel units.
152 /// @note Performs internal caching so only the initial call incurs actual computation.
153 Real totGaussianCurvature(bool useWorldUnits = true);
154
155 /// @brief Compute the average mean curvature of the level set surface.
156 /// @param useWorldUnits Specifies if the result is in world or voxel units.
157 /// @note Performs internal caching so only the initial call incurs actual computation.
158 ✗ Real avgMeanCurvature(bool useWorldUnits = true) {return this->totMeanCurvature(useWorldUnits) / this->area(useWorldUnits);}
159
160 /// @brief Compute the average gaussian curvature of the level set surface.
161 /// @param useWorldUnits Specifies if the result is in world or voxel units.
162 /// @note Performs internal caching so only the initial call incurs actual computation.
163 ✗ Real avgGaussianCurvature(bool useWorldUnits = true) {return this->totGaussianCurvature(useWorldUnits) / this->area(useWorldUnits); }
164
165 /// @brief Compute the Euler characteristic of the level set surface.
166 /// @note Performs internal caching so only the initial call incurs actual computation.
167 int eulerCharacteristic();
168
169 /// @brief Compute the genus of the level set surface.
170 /// @warning The genus is only well defined for a single connected surface.
171 /// @note Performs internal caching so only the initial call incurs actual computation.
172
6/12
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
 8 int genus() { return 1 - this->eulerCharacteristic()/2;}
173
174 private:
175
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[]>;
182
183 // disallow copy construction and copy by assignment!
184 LevelSetMeasure(const LevelSetMeasure&);// not implemented
185 LevelSetMeasure& operator=(const LevelSetMeasure&);// not implemented
186
187 const GridType *mGrid;
188 ManagerPtr mLeafs;
189 BufferPtr mBuffer;
190 InterruptT *mInterrupter;
191 double mDx, mArea, mVolume, mTotMeanCurvature, mTotGausCurvature;
192 int mGrainSize;
193 bool mUpdateArea, mUpdateCurvature;
194
195 // @brief Return false if the process was interrupted
196 bool checkInterrupter();
197
198 ✗ struct MeasureArea
199 {
200 8 MeasureArea(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
201 {
202
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 8 if (parent->mInterrupter) parent->mInterrupter->start("Measuring area and volume of level set");
203
1/2
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
 8 if (parent->mGrainSize>0) {
204
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 8 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
205 } else {
206 ✗ (*this)(parent->mLeafs->leafRange());
207 }
208
1/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 12 tbb::parallel_invoke([&](){parent->mArea = parent->reduce(0);},
209 4 [&](){parent->mVolume = parent->reduce(1)/3.0;});
210 8 parent->mUpdateArea = false;
211
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 8 if (parent->mInterrupter) parent->mInterrupter->end();
212 }
213 110 MeasureArea(const MeasureArea& other) : mParent(other.mParent), mStencil(*mParent->mGrid) {}
214 void operator()(const LeafRange& range) const;
215 LevelSetMeasure* mParent;
216 mutable math::GradStencil<GridT, false> mStencil;
217 };// MeasureArea
218
219 ✗ struct MeasureCurvatures
220 {
221 32 MeasureCurvatures(LevelSetMeasure* parent) : mParent(parent), mStencil(*mParent->mGrid)
222 {
223
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 32 if (parent->mInterrupter) parent->mInterrupter->start("Measuring curvatures of level set");
224
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 32 if (parent->mGrainSize>0) {
225
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 32 tbb::parallel_for(parent->mLeafs->leafRange(parent->mGrainSize), *this);
226 } else {
227 ✗ (*this)(parent->mLeafs->leafRange());
228 }
229
1/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 48 tbb::parallel_invoke([&](){parent->mTotMeanCurvature = parent->reduce(0);},
230 16 [&](){parent->mTotGausCurvature = parent->reduce(1);});
231 32 parent->mUpdateCurvature = false;
232
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 32 if (parent->mInterrupter) parent->mInterrupter->end();
233 }
234 441 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;
238 };// MeasureCurvatures
239
240 80 double reduce(int offset)
241 {
242 80 double *first = mBuffer.get() + offset*mLeafs->leafCount(), *last = first + mLeafs->leafCount();
243 tbb::parallel_sort(first, last);// mitigates catastrophic cancellation
244 Real sum = 0.0;
245
2/2
✓ Branch 0 taken 256018 times.
✓ Branch 1 taken 40 times.
 512116 while(first != last) sum += *first++;
246 80 return sum;
247 }
248
249 }; // end of LevelSetMeasure class
250
251
252 template<typename GridT, typename InterruptT>
253 inline
254 38 LevelSetMeasure<GridT, InterruptT>::LevelSetMeasure(const GridType& grid, InterruptT* interrupt)
255 : mInterrupter(interrupt)
256
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 3 times.
 44 , mGrainSize(1)
257 {
258
2/2
✓ Branch 1 taken 16 times.
✓ Branch 2 taken 3 times.
 38 this->init(grid);
259 }
260
261 template<typename GridT, typename InterruptT>
262 inline void
263 42 LevelSetMeasure<GridT, InterruptT>::init(const GridType& grid)
264 {
265
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 21 times.
 42 if (!grid.hasUniformVoxels()) {
266 ✗ OPENVDB_THROW(RuntimeError,
267 "The transform must have uniform scale for the LevelSetMeasure to function");
268 }
269
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 21 times.
 42 if (grid.getGridClass() != GRID_LEVEL_SET) {
270 ✗ OPENVDB_THROW(RuntimeError,
271 "LevelSetMeasure only supports level sets;"
272 " try setting the grid class to \"level set\"");
273 }
274
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 18 times.
 42 if (grid.empty()) {
275
2/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 24 OPENVDB_THROW(RuntimeError,
276 "LevelSetMeasure does not support empty grids;");
277 }
278 36 mGrid = &grid;
279 36 mDx = grid.voxelSize()[0];
280 36 mLeafs = std::make_unique<ManagerType>(mGrid->tree());
281 36 mBuffer = std::make_unique<double[]>(2*mLeafs->leafCount());
282 36 mUpdateArea = mUpdateCurvature = true;
283 }
284
285 template<typename GridT, typename InterruptT>
286 inline Real
287 36 LevelSetMeasure<GridT, InterruptT>::area(bool useWorldUnits)
288 {
289
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 15 times.
 36 if (mUpdateArea) {MeasureArea m(this);};
290 36 double area = mArea;
291
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 9 times.
 36 if (useWorldUnits) area *= math::Pow2(mDx);
292 36 return area;
293 }
294
295 template<typename GridT, typename InterruptT>
296 inline Real
297 16 LevelSetMeasure<GridT, InterruptT>::volume(bool useWorldUnits)
298 {
299
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 7 times.
 16 if (mUpdateArea) {MeasureArea m(this);};
300 16 double volume = mVolume;
301
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 4 times.
 16 if (useWorldUnits) volume *= math::Pow3(mDx) ;
302 16 return volume;
303 }
304
305 template<typename GridT, typename InterruptT>
306 inline Real
307 ✗ LevelSetMeasure<GridT, InterruptT>::totMeanCurvature(bool useWorldUnits)
308 {
309 ✗ if (mUpdateCurvature) {MeasureCurvatures m(this);};
310 ✗ return mTotMeanCurvature * (useWorldUnits ? mDx : 1);
311 }
312
313 template<typename GridT, typename InterruptT>
314 inline Real
315 28 LevelSetMeasure<GridT, InterruptT>::totGaussianCurvature(bool)
316 {
317
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 28 if (mUpdateCurvature) {MeasureCurvatures m(this);};
318 28 return mTotGausCurvature;
319 }
320
321 template<typename GridT, typename InterruptT>
322 inline int
323 28 LevelSetMeasure<GridT, InterruptT>::eulerCharacteristic()
324 {
325 28 const Real x = this->totGaussianCurvature(true) / (2.0*math::pi<Real>());
326 28 return int(math::Round( x ));
327 }
328
329 ///////////////////////// PRIVATE METHODS //////////////////////
330
331 template<typename GridT, typename InterruptT>
332 inline bool
333 5354 LevelSetMeasure<GridT, InterruptT>::checkInterrupter()
334 {
335
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2677 times.
 5354 if (util::wasInterrupted(mInterrupter)) {
336 ✗ thread::cancelGroupExecution();
337 ✗ return false;
338 }
339 return true;
340 }
341
342 template<typename GridT, typename InterruptT>
343 inline void
344 1070 LevelSetMeasure<GridT, InterruptT>::
345 MeasureArea::operator()(const LeafRange& range) const
346 {
347 using Vec3T = math::Vec3<ValueType>;
348 // computations are performed in index space where dV = 1
349 1070 mParent->checkInterrupter();
350 1070 const Real invDx = 1.0/mParent->mDx;
351 const DiracDelta<Real> DD(1.5);// dirac delta function is 3 voxel units wide
352 const size_t leafCount = mParent->mLeafs->leafCount();
353
2/2
✓ Branch 1 taken 26406 times.
✓ Branch 2 taken 535 times.
 53882 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
354 Real sumA = 0, sumV = 0;//reduce risk of catastrophic cancellation
355
2/2
✓ Branch 0 taken 4887712 times.
✓ Branch 1 taken 26406 times.
 9828236 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
356
2/2
✓ Branch 0 taken 2441734 times.
✓ Branch 1 taken 2445978 times.
 9775424 const Real dd = DD(invDx * (*voxelIter));
357
2/2
✓ Branch 0 taken 2441734 times.
✓ Branch 1 taken 2445978 times.
 9775424 if (dd > 0.0) {
358 4883468 mStencil.moveTo(voxelIter);
359 const Coord& p = mStencil.getCenterCoord();// in voxel units
360 4883468 const Vec3T g = mStencil.gradient();// in world units
361 4883468 sumA += dd*g.length();// \delta(\phi)*|\nabla\phi|
362 4883468 sumV += dd*(g[0]*Real(p[0]) + g[1]*Real(p[1]) + g[2]*Real(p[2]));// \delta(\phi)\vec{x}\cdot\nabla\phi
363 }
364 }
365 52812 double* ptr = mParent->mBuffer.get() + leafIter.pos();
366 52812 *ptr = sumA;
367 52812 ptr += leafCount;
368 52812 *ptr = sumV;
369 }
370 }
371
372 template<typename GridT, typename InterruptT>
373 inline void
374 4284 LevelSetMeasure<GridT, InterruptT>::
375 MeasureCurvatures::operator()(const LeafRange& range) const
376 {
377 using Vec3T = math::Vec3<ValueType>;
378 // computations are performed in index space where dV = 1
379 4284 mParent->checkInterrupter();
380 4284 const Real dx = mParent->mDx, dx2=dx*dx, invDx = 1.0/dx;
381 const DiracDelta<Real> DD(1.5);// dirac delta function is 3 voxel units wide
382 ValueType mean, gauss;
383 const size_t leafCount = mParent->mLeafs->leafCount();
384
2/2
✓ Branch 1 taken 101603 times.
✓ Branch 2 taken 2142 times.
 207490 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
385 Real sumM = 0, sumG = 0;//reduce risk of catastrophic cancellation
386
2/2
✓ Branch 0 taken 19115750 times.
✓ Branch 1 taken 101603 times.
 38434706 for (VoxelCIterT voxelIter = leafIter->cbeginValueOn(); voxelIter; ++voxelIter) {
387
2/2
✓ Branch 0 taken 9640106 times.
✓ Branch 1 taken 9475644 times.
 38231500 const Real dd = DD(invDx * (*voxelIter));
388
2/2
✓ Branch 0 taken 9640106 times.
✓ Branch 1 taken 9475644 times.
 38231500 if (dd > 0.0) {
389 19280212 mStencil.moveTo(voxelIter);
390 19280212 const Vec3T g = mStencil.gradient();
391 19280212 const Real dA = dd*g.length();// \delta(\phi)*\delta(\phi)
392 19280212 mStencil.curvatures(mean, gauss);
393 19280212 sumM += dA*mean*dx;// \delta(\phi)*\delta(\phi)*MeanCurvature
394 19280212 sumG += dA*gauss*dx2;// \delta(\phi)*\delta(\phi)*GaussCurvature
395 }
396 }
397 203206 double* ptr = mParent->mBuffer.get() + leafIter.pos();
398 203206 *ptr = sumM;
399 203206 ptr += leafCount;
400 203206 *ptr = sumG;
401 }
402 }
403
404 ////////////////////////////////////////
405
406 //{
407 /// @cond OPENVDB_DOCS_INTERNAL
408
409 template<class GridT>
410 inline
411 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, Real>::type
412 2 doLevelSetArea(const GridT& grid, bool useWorldUnits)
413 {
414 4 LevelSetMeasure<GridT> m(grid);
415
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 4 return m.area(useWorldUnits);
416 }
417
418 template<class GridT>
419 inline
420 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, Real>::type
421 doLevelSetArea(const GridT&, bool)
422 {
423 OPENVDB_THROW(TypeError,
424 "level set area is supported only for scalar, floating-point grids");
425 }
426
427 /// @endcond
428 //}
429
430 template<class GridT>
431 Real
432 2 levelSetArea(const GridT& grid, bool useWorldUnits)
433 {
434 2 return doLevelSetArea<GridT>(grid, useWorldUnits);
435 }
436
437 ////////////////////////////////////////
438
439 //{
440 /// @cond OPENVDB_DOCS_INTERNAL
441
442 template<class GridT>
443 inline
444 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, Real>::type
445 2 doLevelSetVolume(const GridT& grid, bool useWorldUnits)
446 {
447 4 LevelSetMeasure<GridT> m(grid);
448
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 4 return m.volume(useWorldUnits);
449 }
450
451 template<class GridT>
452 inline
453 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, Real>::type
454 doLevelSetVolume(const GridT&, bool)
455 {
456 OPENVDB_THROW(TypeError,
457 "level set volume is supported only for scalar, floating-point grids");
458 }
459
460 /// @endcond
461 //}
462
463 template<class GridT>
464 Real
465 2 levelSetVolume(const GridT& grid, bool useWorldUnits)
466 {
467 2 return doLevelSetVolume<GridT>(grid, useWorldUnits);
468 }
469
470 ////////////////////////////////////////
471
472 //{
473 /// @cond OPENVDB_DOCS_INTERNAL
474
475 template<class GridT>
476 inline
477 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
478 18 doLevelSetEulerCharacteristic(const GridT& grid)
479 {
480 36 LevelSetMeasure<GridT> m(grid);
481
1/2
✓ Branch 1 taken 9 times.
✗ Branch 2 not taken.
 36 return m.eulerCharacteristic();
482 }
483
484 template<class GridT>
485 inline
486 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, int>::type
487 doLevelSetEulerCharacteristic(const GridT&)
488 {
489 OPENVDB_THROW(TypeError,
490 "level set euler characteristic is supported only for scalar, floating-point grids");
491 }
492
493 /// @endcond
494 //}
495
496
497 template<class GridT>
498 int
499 18 levelSetEulerCharacteristic(const GridT& grid)
500 {
501 18 return doLevelSetEulerCharacteristic(grid);
502 }
503
504 ////////////////////////////////////////
505
506 //{
507 /// @cond OPENVDB_DOCS_INTERNAL
508
509 template<class GridT>
510 inline
511 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
512 doLevelSetEuler(const GridT& grid)
513 {
514 LevelSetMeasure<GridT> m(grid);
515 return m.eulerCharacteristics();
516
517 }
518
519 template<class GridT>
520 inline
521 typename std::enable_if<std::is_floating_point<typename GridT::ValueType>::value, int>::type
522 16 doLevelSetGenus(const GridT& grid)
523 {
524 26 LevelSetMeasure<GridT> m(grid);
525 10 return m.genus();
526 }
527
528 template<class GridT>
529 inline
530 typename std::enable_if<!std::is_floating_point<typename GridT::ValueType>::value, int>::type
531 doLevelSetGenus(const GridT&)
532 {
533 OPENVDB_THROW(TypeError,
534 "level set genus is supported only for scalar, floating-point grids");
535 }
536
537 /// @endcond
538 //}
539
540 template<class GridT>
541 int
542 8 levelSetGenus(const GridT& grid)
543 {
544 8 return doLevelSetGenus(grid);
545 }
546
547
548 ////////////////////////////////////////
549
550
551 // Explicit Template Instantiation
552
553 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
554
555 #ifdef OPENVDB_INSTANTIATE_LEVELSETMEASURE
556 #include <openvdb/util/ExplicitInstantiation.h>
557 #endif
558
559 #define _FUNCTION(TreeT) \
560 Real levelSetArea(const Grid<TreeT>&, bool)
561 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
562 #undef _FUNCTION
563
564 #define _FUNCTION(TreeT) \
565 Real levelSetVolume(const Grid<TreeT>&, bool)
566 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
567 #undef _FUNCTION
568
569 #define _FUNCTION(TreeT) \
570 int levelSetEulerCharacteristic(const Grid<TreeT>&)
571 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
572 #undef _FUNCTION
573
574 #define _FUNCTION(TreeT) \
575 int levelSetGenus(const Grid<TreeT>&)
576 OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
577 #undef _FUNCTION
578
579 OPENVDB_INSTANTIATE_CLASS LevelSetMeasure<FloatGrid, util::NullInterrupter>;
580 OPENVDB_INSTANTIATE_CLASS LevelSetMeasure<DoubleGrid, util::NullInterrupter>;
581
582 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
583
584
585 } // namespace tools
586 } // namespace OPENVDB_VERSION_NAME
587 } // namespace openvdb
588
589 #endif // OPENVDB_TOOLS_LEVELSETMEASURE_HAS_BEEN_INCLUDED
590