HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LevelSetUtil.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 /// @file tools/LevelSetUtil.h
32 ///
33 /// @brief Miscellaneous utility methods that operate primarily
34 /// or exclusively on level set grids.
35 ///
36 /// @author Mihai Alden
37 
38 #ifndef OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
39 #define OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
40 
41 #include "MeshToVolume.h" // for traceExteriorBoundaries
42 #include "SignedFloodFill.h" // for signedFloodFillWithValues
43 
44 #include <openvdb/Types.h>
45 #include <openvdb/Grid.h>
46 #include <hboost/mpl/at.hpp>
47 #include <hboost/mpl/int.hpp>
48 #include <tbb/blocked_range.h>
49 #include <tbb/parallel_for.h>
50 #include <tbb/parallel_reduce.h>
51 #include <tbb/parallel_sort.h>
52 #include <algorithm>
53 #include <cmath>
54 #include <cstdlib>
55 #include <deque>
56 #include <limits>
57 #include <memory>
58 #include <set>
59 #include <vector>
60 
61 
62 namespace openvdb {
64 namespace OPENVDB_VERSION_NAME {
65 namespace tools {
66 
67 // MS Visual C++ requires this extra level of indirection in order to compile
68 // THIS MUST EXIST IN AN UNNAMED NAMESPACE IN ORDER TO COMPILE ON WINDOWS
69 namespace {
70 
71 template<typename GridType>
72 inline typename GridType::ValueType lsutilGridMax()
73 {
75 }
76 
77 template<typename GridType>
78 inline typename GridType::ValueType lsutilGridZero()
79 {
80  return zeroVal<typename GridType::ValueType>();
81 }
82 
83 } // unnamed namespace
84 
85 
86 ////////////////////////////////////////
87 
88 
89 /// @brief Threaded method to convert a sparse level set/SDF into a sparse fog volume
90 ///
91 /// @details For a level set, the active and negative-valued interior half of the
92 /// narrow band becomes a linear ramp from 0 to 1; the inactive interior becomes
93 /// active with a constant value of 1; and the exterior, including the background
94 /// and the active exterior half of the narrow band, becomes inactive with a constant
95 /// value of 0. The interior, though active, remains sparse.
96 /// @details For a generic SDF, a specified cutoff distance determines the width
97 /// of the ramp, but otherwise the result is the same as for a level set.
98 ///
99 /// @param grid level set/SDF grid to transform
100 /// @param cutoffDistance optional world space cutoff distance for the ramp
101 /// (automatically clamped if greater than the interior
102 /// narrow band width)
103 template<class GridType>
104 inline void
106  GridType& grid,
107  typename GridType::ValueType cutoffDistance = lsutilGridMax<GridType>());
108 
109 
110 /// @brief Threaded method to construct a boolean mask that represents interior regions
111 /// in a signed distance field.
112 ///
113 /// @return A shared pointer to either a boolean grid or tree with the same tree
114 /// configuration and potentially transform as the input @c volume and whose active
115 /// and @c true values correspond to the interior of the input signed distance field.
116 ///
117 /// @param volume Signed distance field / level set volume.
118 /// @param isovalue Threshold below which values are considered part of the
119 /// interior region.
120 template<class GridOrTreeType>
121 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
123  const GridOrTreeType& volume,
124  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>());
125 
126 
127 /// @brief Extracts the interior regions of a signed distance field and topologically enclosed
128 /// (watertight) regions of value greater than the @a isovalue (cavities) that can arise
129 /// as the result of CSG union operations between different shapes where at least one of
130 /// the shapes has a concavity that is capped.
131 ///
132 /// For example the enclosed region of a capped bottle would include the walls and
133 /// the interior cavity.
134 ///
135 /// @return A shared pointer to either a boolean grid or tree with the same tree configuration
136 /// and potentially transform as the input @c volume and whose active and @c true values
137 /// correspond to the interior and enclosed regions in the input signed distance field.
138 ///
139 /// @param volume Signed distance field / level set volume.
140 /// @param isovalue Threshold below which values are considered part of the interior region.
141 /// @param fillMask Optional boolean tree, when provided enclosed cavity regions that are not
142 /// completely filled by this mask are ignored.
143 ///
144 /// For instance if the fill mask does not completely fill the bottle in the
145 /// previous example only the walls and cap are returned and the interior
146 /// cavity will be ignored.
147 template<typename GridOrTreeType>
148 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
150  const GridOrTreeType& volume,
151  typename GridOrTreeType::ValueType isovalue = lsutilGridZero<GridOrTreeType>(),
152  const typename TreeAdapter<GridOrTreeType>::TreeType::template ValueConverter<bool>::Type*
153  fillMask = nullptr);
154 
155 
156 /// @brief Return a mask of the voxels that intersect the implicit surface with
157 /// the given @a isovalue.
158 ///
159 /// @param volume Signed distance field / level set volume.
160 /// @param isovalue The crossing point that is considered the surface.
161 template<typename GridOrTreeType>
162 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
163 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue);
164 
165 
166 /// @brief Return a mask for each connected component of the given grid's active voxels.
167 ///
168 /// @param volume Input grid or tree
169 /// @param masks Output set of disjoint active topology masks sorted in descending order
170 /// based on the active voxel count.
171 template<typename GridOrTreeType>
172 inline void
173 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
174  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks);
175 
176 
177 /// @brief Separates disjoint active topology components into distinct grids or trees.
178 ///
179 /// @details Supports volumes with active tiles.
180 ///
181 /// @param volume Input grid or tree
182 /// @param segments Output set of disjoint active topology components sorted in
183 /// descending order based on the active voxel count.
184 template<typename GridOrTreeType>
185 inline void
186 segmentActiveVoxels(const GridOrTreeType& volume,
187  std::vector<typename GridOrTreeType::Ptr>& segments);
188 
189 
190 /// @brief Separates disjoint SDF surfaces into distinct grids or trees.
191 ///
192 /// @details Supports asymmetric interior / exterior narrowband widths and
193 /// SDF volumes with dense interior regions.
194 ///
195 /// @param volume Input signed distance field / level set volume
196 /// @param segments Output set of disjoint SDF surfaces found in @a volume sorted in
197 /// descending order based on the surface intersecting voxel count.
198 template<typename GridOrTreeType>
199 inline void
200 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments);
201 
202 
203 ////////////////////////////////////////////////////////////////////////////////
204 ////////////////////////////////////////////////////////////////////////////////
205 
206 // Internal utility objects and implementation details
207 
208 
209 namespace level_set_util_internal {
210 
211 
212 template<typename LeafNodeType>
214 
215  using ValueType = typename LeafNodeType::ValueType;
217 
219  ValueType isovalue, const LeafNodeType ** nodes, BoolLeafNodeType ** maskNodes)
220  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
221  {
222  }
223 
224  void operator()(const tbb::blocked_range<size_t>& range) const {
225 
226  BoolLeafNodeType * maskNodePt = nullptr;
227 
228  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
229 
230  mMaskNodes[n] = nullptr;
231  const LeafNodeType& node = *mNodes[n];
232 
233  if (!maskNodePt) {
234  maskNodePt = new BoolLeafNodeType(node.origin(), false);
235  } else {
236  maskNodePt->setOrigin(node.origin());
237  }
238 
239  const ValueType* values = &node.getValue(0);
240  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
241  if (values[i] < mIsovalue) maskNodePt->setValueOn(i, true);
242  }
243 
244  if (maskNodePt->onVoxelCount() > 0) {
245  mMaskNodes[n] = maskNodePt;
246  maskNodePt = nullptr;
247  }
248  }
249 
250  if (maskNodePt) delete maskNodePt;
251  }
252 
253  LeafNodeType const * const * const mNodes;
256 }; // MaskInteriorVoxels
257 
258 
259 template<typename TreeType, typename InternalNodeType>
261 
262  using ValueType = typename TreeType::ValueType;
263 
264  MaskInteriorTiles(ValueType isovalue, const TreeType& tree, InternalNodeType ** maskNodes)
265  : mTree(&tree), mMaskNodes(maskNodes), mIsovalue(isovalue) { }
266 
267  void operator()(const tbb::blocked_range<size_t>& range) const {
269  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
270  typename InternalNodeType::ValueAllIter it = mMaskNodes[n]->beginValueAll();
271  for (; it; ++it) {
272  if (acc.getValue(it.getCoord()) < mIsovalue) {
273  it.setValue(true);
274  it.setValueOn(true);
275  }
276  }
277  }
278  }
279 
280  TreeType const * const mTree;
281  InternalNodeType ** const mMaskNodes;
283 }; // MaskInteriorTiles
284 
285 
286 template<typename TreeType>
287 struct PopulateTree {
288 
289  using ValueType = typename TreeType::ValueType;
290  using LeafNodeType = typename TreeType::LeafNodeType;
291 
292  PopulateTree(TreeType& tree, LeafNodeType** leafnodes,
293  const size_t * nodexIndexMap, ValueType background)
294  : mNewTree(background)
295  , mTreePt(&tree)
296  , mNodes(leafnodes)
297  , mNodeIndexMap(nodexIndexMap)
298  {
299  }
300 
301  PopulateTree(PopulateTree& rhs, tbb::split)
302  : mNewTree(rhs.mNewTree.background())
303  , mTreePt(&mNewTree)
304  , mNodes(rhs.mNodes)
305  , mNodeIndexMap(rhs.mNodeIndexMap)
306  {
307  }
308 
309  void operator()(const tbb::blocked_range<size_t>& range) {
310 
311  tree::ValueAccessor<TreeType> acc(*mTreePt);
312 
313  if (mNodeIndexMap) {
314  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
315  for (size_t i = mNodeIndexMap[n], I = mNodeIndexMap[n + 1]; i < I; ++i) {
316  if (mNodes[i] != nullptr) acc.addLeaf(mNodes[i]);
317  }
318  }
319  } else {
320  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
321  acc.addLeaf(mNodes[n]);
322  }
323  }
324  }
325 
326  void join(PopulateTree& rhs) { mTreePt->merge(*rhs.mTreePt); }
327 
328 private:
329  TreeType mNewTree;
330  TreeType * const mTreePt;
331  LeafNodeType ** const mNodes;
332  size_t const * const mNodeIndexMap;
333 }; // PopulateTree
334 
335 
336 /// @brief Negative active values are set @c 0, everything else is set to @c 1.
337 template<typename LeafNodeType>
339 
340  using ValueType = typename LeafNodeType::ValueType;
342 
344  ValueType isovalue, const LeafNodeType ** nodes, CharLeafNodeType ** maskNodes)
345  : mNodes(nodes), mMaskNodes(maskNodes), mIsovalue(isovalue)
346  {
347  }
348 
349  void operator()(const tbb::blocked_range<size_t>& range) const {
350 
351  CharLeafNodeType * maskNodePt = nullptr;
352 
353  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
354 
355  mMaskNodes[n] = nullptr;
356  const LeafNodeType& node = *mNodes[n];
357 
358  if (!maskNodePt) {
359  maskNodePt = new CharLeafNodeType(node.origin(), 1);
360  } else {
361  maskNodePt->setOrigin(node.origin());
362  }
363 
364  typename LeafNodeType::ValueOnCIter it;
365  for (it = node.cbeginValueOn(); it; ++it) {
366  maskNodePt->setValueOn(it.pos(), ((*it - mIsovalue) < 0.0) ? 0 : 1);
367  }
368 
369  if (maskNodePt->onVoxelCount() > 0) {
370  mMaskNodes[n] = maskNodePt;
371  maskNodePt = nullptr;
372  }
373  }
374 
375  if (maskNodePt) delete maskNodePt;
376  }
377 
378  LeafNodeType const * const * const mNodes;
381 }; // LabelBoundaryVoxels
382 
383 
384 template<typename LeafNodeType>
386  using ValueType = typename LeafNodeType::ValueType;
387 
388  FlipRegionSign(LeafNodeType ** nodes) : mNodes(nodes) { }
389 
390  void operator()(const tbb::blocked_range<size_t>& range) const {
391  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
392  ValueType* values = const_cast<ValueType*>(&mNodes[n]->getValue(0));
393  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
394  values[i] = values[i] < 0 ? 1 : -1;
395  }
396  }
397  }
398 
399  LeafNodeType ** const mNodes;
400 }; // FlipRegionSign
401 
402 
403 template<typename LeafNodeType>
405 
406  using ValueType = typename LeafNodeType::ValueType;
407 
408  FindMinVoxelValue(LeafNodeType const * const * const leafnodes)
409  : minValue(std::numeric_limits<ValueType>::max())
410  , mNodes(leafnodes)
411  {
412  }
413 
415  : minValue(std::numeric_limits<ValueType>::max())
416  , mNodes(rhs.mNodes)
417  {
418  }
419 
420  void operator()(const tbb::blocked_range<size_t>& range) {
421  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
422  const ValueType* data = mNodes[n]->buffer().data();
423  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
424  minValue = std::min(minValue, data[i]);
425  }
426  }
427  }
428 
430 
432 
433  LeafNodeType const * const * const mNodes;
434 }; // FindMinVoxelValue
435 
436 
437 template<typename InternalNodeType>
439 
440  using ValueType = typename InternalNodeType::ValueType;
441 
442  FindMinTileValue(InternalNodeType const * const * const nodes)
443  : minValue(std::numeric_limits<ValueType>::max())
444  , mNodes(nodes)
445  {
446  }
447 
449  : minValue(std::numeric_limits<ValueType>::max())
450  , mNodes(rhs.mNodes)
451  {
452  }
453 
454  void operator()(const tbb::blocked_range<size_t>& range) {
455  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
456  typename InternalNodeType::ValueAllCIter it = mNodes[n]->beginValueAll();
457  for (; it; ++it) {
458  minValue = std::min(minValue, *it);
459  }
460  }
461  }
462 
464 
466 
467  InternalNodeType const * const * const mNodes;
468 }; // FindMinTileValue
469 
470 
471 template<typename LeafNodeType>
473 
474  using ValueType = typename LeafNodeType::ValueType;
475 
476  SDFVoxelsToFogVolume(LeafNodeType ** nodes, ValueType cutoffDistance)
477  : mNodes(nodes), mWeight(ValueType(1.0) / cutoffDistance)
478  {
479  }
480 
481  void operator()(const tbb::blocked_range<size_t>& range) const {
482 
483  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
484 
485  LeafNodeType& node = *mNodes[n];
486  node.setValuesOff();
487 
488  ValueType* values = node.buffer().data();
489  for (Index i = 0; i < LeafNodeType::SIZE; ++i) {
490  values[i] = values[i] > ValueType(0.0) ? ValueType(0.0) : values[i] * mWeight;
491  if (values[i] > ValueType(0.0)) node.setValueOn(i);
492  }
493 
494  if (node.onVoxelCount() == 0) {
495  delete mNodes[n];
496  mNodes[n] = nullptr;
497  }
498  }
499  }
500 
501  LeafNodeType ** const mNodes;
503 }; // SDFVoxelsToFogVolume
504 
505 
506 template<typename TreeType, typename InternalNodeType>
508 
509  SDFTilesToFogVolume(const TreeType& tree, InternalNodeType ** nodes)
510  : mTree(&tree), mNodes(nodes) { }
511 
512  void operator()(const tbb::blocked_range<size_t>& range) const {
513 
514  using ValueType = typename TreeType::ValueType;
516 
517  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
518  typename InternalNodeType::ValueAllIter it = mNodes[n]->beginValueAll();
519  for (; it; ++it) {
520  if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
521  it.setValue(ValueType(1.0));
522  it.setValueOn(true);
523  }
524  }
525  }
526  }
527 
528  TreeType const * const mTree;
529  InternalNodeType ** const mNodes;
530 }; // SDFTilesToFogVolume
531 
532 
533 template<typename TreeType>
535 
536  using ValueType = typename TreeType::ValueType;
537  using LeafNodeType = typename TreeType::LeafNodeType;
538  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
539  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
540 
541  FillMaskBoundary(const TreeType& tree, ValueType isovalue, const BoolTreeType& fillMask,
542  const BoolLeafNodeType ** fillNodes, BoolLeafNodeType ** newNodes)
543  : mTree(&tree)
544  , mFillMask(&fillMask)
545  , mFillNodes(fillNodes)
546  , mNewNodes(newNodes)
547  , mIsovalue(isovalue)
548  {
549  }
550 
551  void operator()(const tbb::blocked_range<size_t>& range) const {
552 
553  tree::ValueAccessor<const BoolTreeType> maskAcc(*mFillMask);
554  tree::ValueAccessor<const TreeType> distAcc(*mTree);
555 
556  std::unique_ptr<char[]> valueMask(new char[BoolLeafNodeType::SIZE]);
557 
558  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
559 
560  mNewNodes[n] = nullptr;
561  const BoolLeafNodeType& node = *mFillNodes[n];
562  const Coord& origin = node.origin();
563 
564  const bool denseNode = node.isDense();
565 
566  // possible early out if the fill mask is dense
567  if (denseNode) {
568 
569  int denseNeighbors = 0;
570 
571  const BoolLeafNodeType* neighborNode =
572  maskAcc.probeConstLeaf(origin.offsetBy(-1, 0, 0));
573  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
574 
575  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(BoolLeafNodeType::DIM, 0, 0));
576  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
577 
578  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, -1, 0));
579  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
580 
581  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, BoolLeafNodeType::DIM, 0));
582  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
583 
584  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, -1));
585  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
586 
587  neighborNode = maskAcc.probeConstLeaf(origin.offsetBy(0, 0, BoolLeafNodeType::DIM));
588  if (neighborNode && neighborNode->isDense()) ++denseNeighbors;
589 
590  if (denseNeighbors == 6) continue;
591  }
592 
593  // rest value mask
594  memset(valueMask.get(), 0, sizeof(char) * BoolLeafNodeType::SIZE);
595 
596  const typename TreeType::LeafNodeType* distNode = distAcc.probeConstLeaf(origin);
597 
598  // check internal voxel neighbors
599 
600  bool earlyTermination = false;
601 
602  if (!denseNode) {
603  if (distNode) {
604  evalInternalNeighborsP(valueMask.get(), node, *distNode);
605  evalInternalNeighborsN(valueMask.get(), node, *distNode);
606  } else if (distAcc.getValue(origin) > mIsovalue) {
607  earlyTermination = evalInternalNeighborsP(valueMask.get(), node);
608  if (!earlyTermination) {
609  earlyTermination = evalInternalNeighborsN(valueMask.get(), node);
610  }
611  }
612  }
613 
614  // check external voxel neighbors
615 
616  if (!earlyTermination) {
617  evalExternalNeighborsX<true>(valueMask.get(), node, maskAcc, distAcc);
618  evalExternalNeighborsX<false>(valueMask.get(), node, maskAcc, distAcc);
619  evalExternalNeighborsY<true>(valueMask.get(), node, maskAcc, distAcc);
620  evalExternalNeighborsY<false>(valueMask.get(), node, maskAcc, distAcc);
621  evalExternalNeighborsZ<true>(valueMask.get(), node, maskAcc, distAcc);
622  evalExternalNeighborsZ<false>(valueMask.get(), node, maskAcc, distAcc);
623  }
624 
625  // Export marked boundary voxels.
626 
627  int numBoundaryValues = 0;
628  for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
629  numBoundaryValues += valueMask[i] == 1;
630  }
631 
632  if (numBoundaryValues > 0) {
633  mNewNodes[n] = new BoolLeafNodeType(origin, false);
634  for (Index i = 0, I = BoolLeafNodeType::SIZE; i < I; ++i) {
635  if (valueMask[i] == 1) mNewNodes[n]->setValueOn(i);
636  }
637  }
638  }
639  }
640 
641 private:
642  // Check internal voxel neighbors in positive {x, y, z} directions.
643  void evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node,
644  const LeafNodeType& distNode) const
645  {
646  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
647  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
648  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
649  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
650  for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
651  const Index pos = yPos + z;
652 
653  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
654 
655  if (!node.isValueOn(pos + 1) && distNode.getValue(pos + 1) > mIsovalue) {
656  valueMask[pos] = 1;
657  }
658  }
659  }
660  }
661 
662  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
663  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
664  for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
665  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
666  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
667  const Index pos = yPos + z;
668 
669  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
670 
671  if (!node.isValueOn(pos + BoolLeafNodeType::DIM) &&
672  distNode.getValue(pos + BoolLeafNodeType::DIM) > mIsovalue) {
673  valueMask[pos] = 1;
674  }
675  }
676  }
677  }
678 
679  for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
680  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
681  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
682  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
683  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
684  const Index pos = yPos + z;
685 
686  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
687 
688  if (!node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
689  (distNode.getValue(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)
690  > mIsovalue))
691  {
692  valueMask[pos] = 1;
693  }
694  }
695  }
696  }
697  }
698 
699  bool evalInternalNeighborsP(char* valueMask, const BoolLeafNodeType& node) const {
700 
701  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
702  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
703  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
704  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
705  for (Index z = 0; z < BoolLeafNodeType::DIM - 1; ++z) {
706  const Index pos = yPos + z;
707 
708  if (node.isValueOn(pos) && !node.isValueOn(pos + 1)) {
709  valueMask[pos] = 1;
710  return true;
711  }
712  }
713  }
714  }
715 
716  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
717  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
718  for (Index y = 0; y < BoolLeafNodeType::DIM - 1; ++y) {
719  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
720  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
721  const Index pos = yPos + z;
722 
723  if (node.isValueOn(pos) && !node.isValueOn(pos + BoolLeafNodeType::DIM)) {
724  valueMask[pos] = 1;
725  return true;
726  }
727  }
728  }
729  }
730 
731  for (Index x = 0; x < BoolLeafNodeType::DIM - 1; ++x) {
732  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
733  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
734  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
735  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
736  const Index pos = yPos + z;
737 
738  if (node.isValueOn(pos) &&
739  !node.isValueOn(pos + BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
740  valueMask[pos] = 1;
741  return true;
742  }
743  }
744  }
745  }
746 
747  return false;
748  }
749 
750  // Check internal voxel neighbors in negative {x, y, z} directions.
751 
752  void evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node,
753  const LeafNodeType& distNode) const
754  {
755  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
756  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
757  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
758  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
759  for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
760  const Index pos = yPos + z;
761 
762  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
763 
764  if (!node.isValueOn(pos - 1) && distNode.getValue(pos - 1) > mIsovalue) {
765  valueMask[pos] = 1;
766  }
767  }
768  }
769  }
770 
771  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
772  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
773  for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
774  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
775  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
776  const Index pos = yPos + z;
777 
778  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
779 
780  if (!node.isValueOn(pos - BoolLeafNodeType::DIM) &&
781  distNode.getValue(pos - BoolLeafNodeType::DIM) > mIsovalue) {
782  valueMask[pos] = 1;
783  }
784  }
785  }
786  }
787 
788  for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
789  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
790  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
791  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
792  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
793  const Index pos = yPos + z;
794 
795  if (valueMask[pos] != 0 || !node.isValueOn(pos)) continue;
796 
797  if (!node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM) &&
798  (distNode.getValue(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)
799  > mIsovalue))
800  {
801  valueMask[pos] = 1;
802  }
803  }
804  }
805  }
806  }
807 
808 
809  bool evalInternalNeighborsN(char* valueMask, const BoolLeafNodeType& node) const {
810 
811  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
812  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
813  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
814  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
815  for (Index z = 1; z < BoolLeafNodeType::DIM; ++z) {
816  const Index pos = yPos + z;
817 
818  if (node.isValueOn(pos) && !node.isValueOn(pos - 1)) {
819  valueMask[pos] = 1;
820  return true;
821  }
822  }
823  }
824  }
825 
826  for (Index x = 0; x < BoolLeafNodeType::DIM; ++x) {
827  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
828  for (Index y = 1; y < BoolLeafNodeType::DIM; ++y) {
829  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
830  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
831  const Index pos = yPos + z;
832 
833  if (node.isValueOn(pos) && !node.isValueOn(pos - BoolLeafNodeType::DIM)) {
834  valueMask[pos] = 1;
835  return true;
836  }
837  }
838  }
839  }
840 
841  for (Index x = 1; x < BoolLeafNodeType::DIM; ++x) {
842  const Index xPos = x << (2 * BoolLeafNodeType::LOG2DIM);
843  for (Index y = 0; y < BoolLeafNodeType::DIM; ++y) {
844  const Index yPos = xPos + (y << BoolLeafNodeType::LOG2DIM);
845  for (Index z = 0; z < BoolLeafNodeType::DIM; ++z) {
846  const Index pos = yPos + z;
847 
848  if (node.isValueOn(pos) &&
849  !node.isValueOn(pos - BoolLeafNodeType::DIM * BoolLeafNodeType::DIM)) {
850  valueMask[pos] = 1;
851  return true;
852  }
853  }
854  }
855  }
856 
857  return false;
858  }
859 
860 
861  // Check external voxel neighbors
862 
863  // If UpWind is true check the X+ oriented node face, else the X- oriented face.
864  template<bool UpWind>
865  void evalExternalNeighborsX(char* valueMask, const BoolLeafNodeType& node,
866  const tree::ValueAccessor<const BoolTreeType>& maskAcc,
867  const tree::ValueAccessor<const TreeType>& distAcc) const {
868 
869  const Coord& origin = node.origin();
870  Coord ijk(0, 0, 0), nijk;
871  int step = -1;
872 
873  if (UpWind) {
874  step = 1;
875  ijk[0] = int(BoolLeafNodeType::DIM) - 1;
876  }
877 
878  const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
879 
880  for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
881  const Index yPos = xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
882 
883  for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
884  const Index pos = yPos + ijk[2];
885 
886  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
887 
888  nijk = origin + ijk.offsetBy(step, 0, 0);
889 
890  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
891  valueMask[pos] = 1;
892  }
893  }
894  }
895  }
896  }
897 
898  // If UpWind is true check the Y+ oriented node face, else the Y- oriented face.
899  template<bool UpWind>
900  void evalExternalNeighborsY(char* valueMask, const BoolLeafNodeType& node,
901  const tree::ValueAccessor<const BoolTreeType>& maskAcc,
902  const tree::ValueAccessor<const TreeType>& distAcc) const {
903 
904  const Coord& origin = node.origin();
905  Coord ijk(0, 0, 0), nijk;
906  int step = -1;
907 
908  if (UpWind) {
909  step = 1;
910  ijk[1] = int(BoolLeafNodeType::DIM) - 1;
911  }
912 
913  const Index yPos = ijk[1] << int(BoolLeafNodeType::LOG2DIM);
914 
915  for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
916  const Index xPos = yPos + (ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM)));
917 
918  for (ijk[2] = 0; ijk[2] < int(BoolLeafNodeType::DIM); ++ijk[2]) {
919  const Index pos = xPos + ijk[2];
920 
921  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
922 
923  nijk = origin + ijk.offsetBy(0, step, 0);
924  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
925  valueMask[pos] = 1;
926  }
927  }
928  }
929  }
930  }
931 
932  // If UpWind is true check the Z+ oriented node face, else the Z- oriented face.
933  template<bool UpWind>
934  void evalExternalNeighborsZ(char* valueMask, const BoolLeafNodeType& node,
935  const tree::ValueAccessor<const BoolTreeType>& maskAcc,
936  const tree::ValueAccessor<const TreeType>& distAcc) const {
937 
938  const Coord& origin = node.origin();
939  Coord ijk(0, 0, 0), nijk;
940  int step = -1;
941 
942  if (UpWind) {
943  step = 1;
944  ijk[2] = int(BoolLeafNodeType::DIM) - 1;
945  }
946 
947  for (ijk[0] = 0; ijk[0] < int(BoolLeafNodeType::DIM); ++ijk[0]) {
948  const Index xPos = ijk[0] << (2 * int(BoolLeafNodeType::LOG2DIM));
949 
950  for (ijk[1] = 0; ijk[1] < int(BoolLeafNodeType::DIM); ++ijk[1]) {
951  const Index pos = ijk[2] + xPos + (ijk[1] << int(BoolLeafNodeType::LOG2DIM));
952 
953  if (valueMask[pos] == 0 && node.isValueOn(pos)) {
954 
955  nijk = origin + ijk.offsetBy(0, 0, step);
956  if (!maskAcc.isValueOn(nijk) && distAcc.getValue(nijk) > mIsovalue) {
957  valueMask[pos] = 1;
958  }
959  }
960  }
961  }
962  }
963 
964  //////////
965 
966  TreeType const * const mTree;
967  BoolTreeType const * const mFillMask;
968  BoolLeafNodeType const * const * const mFillNodes;
969  BoolLeafNodeType ** const mNewNodes;
970  ValueType const mIsovalue;
971 }; // FillMaskBoundary
972 
973 
974 /// @brief Constructs a memory light char tree that represents the exterior region with @c +1
975 /// and the interior regions with @c -1.
976 template <class TreeType>
977 inline typename TreeType::template ValueConverter<char>::Type::Ptr
978 computeEnclosedRegionMask(const TreeType& tree, typename TreeType::ValueType isovalue,
979  const typename TreeType::template ValueConverter<bool>::Type* fillMask)
980 {
981  using LeafNodeType = typename TreeType::LeafNodeType;
982  using RootNodeType = typename TreeType::RootNodeType;
983  using NodeChainType = typename RootNodeType::NodeChainType;
984  using InternalNodeType = typename hboost::mpl::at<NodeChainType, hboost::mpl::int_<1>>::type;
985 
986  using CharTreeType = typename TreeType::template ValueConverter<char>::Type;
987  using CharLeafNodeType = typename CharTreeType::LeafNodeType;
988 
989  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
990  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
991 
992  const TreeType* treePt = &tree;
993 
994  size_t numLeafNodes = 0, numInternalNodes = 0;
995 
996  std::vector<const LeafNodeType*> nodes;
997  std::vector<size_t> leafnodeCount;
998 
999  {
1000  // compute the prefix sum of the leafnode count in each internal node.
1001  std::vector<const InternalNodeType*> internalNodes;
1002  treePt->getNodes(internalNodes);
1003 
1004  numInternalNodes = internalNodes.size();
1005 
1006  leafnodeCount.push_back(0);
1007  for (size_t n = 0; n < numInternalNodes; ++n) {
1008  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
1009  }
1010 
1011  numLeafNodes = leafnodeCount.back();
1012 
1013  // extract all leafnodes
1014  nodes.reserve(numLeafNodes);
1015 
1016  for (size_t n = 0; n < numInternalNodes; ++n) {
1017  internalNodes[n]->getNodes(nodes);
1018  }
1019  }
1020 
1021  // create mask leafnodes
1022  std::unique_ptr<CharLeafNodeType*[]> maskNodes(new CharLeafNodeType*[numLeafNodes]);
1023 
1024  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1025  LabelBoundaryVoxels<LeafNodeType>(isovalue, &nodes[0], maskNodes.get()));
1026 
1027  // create mask grid
1028  typename CharTreeType::Ptr maskTree(new CharTreeType(1));
1029 
1030  PopulateTree<CharTreeType> populate(*maskTree, maskNodes.get(), &leafnodeCount[0], 1);
1031  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1032 
1033  // optionally evaluate the fill mask
1034 
1035  std::vector<CharLeafNodeType*> extraMaskNodes;
1036 
1037  if (fillMask) {
1038 
1039  std::vector<const BoolLeafNodeType*> fillMaskNodes;
1040  fillMask->getNodes(fillMaskNodes);
1041 
1042  std::unique_ptr<BoolLeafNodeType*[]> boundaryMaskNodes(
1043  new BoolLeafNodeType*[fillMaskNodes.size()]);
1044 
1045  tbb::parallel_for(tbb::blocked_range<size_t>(0, fillMaskNodes.size()),
1046  FillMaskBoundary<TreeType>(tree, isovalue, *fillMask, &fillMaskNodes[0],
1047  boundaryMaskNodes.get()));
1048 
1049  tree::ValueAccessor<CharTreeType> maskAcc(*maskTree);
1050 
1051  for (size_t n = 0, N = fillMaskNodes.size(); n < N; ++n) {
1052 
1053  if (boundaryMaskNodes[n] == nullptr) continue;
1054 
1055  const BoolLeafNodeType& boundaryNode = *boundaryMaskNodes[n];
1056  const Coord& origin = boundaryNode.origin();
1057 
1058  CharLeafNodeType* maskNodePt = maskAcc.probeLeaf(origin);
1059 
1060  if (!maskNodePt) {
1061  maskNodePt = maskAcc.touchLeaf(origin);
1062  extraMaskNodes.push_back(maskNodePt);
1063  }
1064 
1065  char* data = maskNodePt->buffer().data();
1066 
1067  typename BoolLeafNodeType::ValueOnCIter it = boundaryNode.cbeginValueOn();
1068  for (; it; ++it) {
1069  if (data[it.pos()] != 0) data[it.pos()] = -1;
1070  }
1071 
1072  delete boundaryMaskNodes[n];
1073  }
1074  }
1075 
1076  // eliminate enclosed regions
1077  tools::traceExteriorBoundaries(*maskTree);
1078 
1079  // flip voxel sign to negative inside and positive outside.
1080  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1081  FlipRegionSign<CharLeafNodeType>(maskNodes.get()));
1082 
1083  if (!extraMaskNodes.empty()) {
1084  tbb::parallel_for(tbb::blocked_range<size_t>(0, extraMaskNodes.size()),
1085  FlipRegionSign<CharLeafNodeType>(&extraMaskNodes[0]));
1086  }
1087 
1088  // propagate sign information into tile region
1089  tools::signedFloodFill(*maskTree);
1090 
1091  return maskTree;
1092 } // computeEnclosedRegionMask()
1093 
1094 
1095 template <class TreeType>
1096 inline typename TreeType::template ValueConverter<bool>::Type::Ptr
1097 computeInteriorMask(const TreeType& tree, typename TreeType::ValueType iso)
1098 {
1099  using ValueType = typename TreeType::ValueType;
1100  using LeafNodeType = typename TreeType::LeafNodeType;
1101  using RootNodeType = typename TreeType::RootNodeType;
1102  using NodeChainType = typename RootNodeType::NodeChainType;
1103  using InternalNodeType = typename hboost::mpl::at<NodeChainType, hboost::mpl::int_<1> >::type;
1104 
1105  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1106  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1107  using BoolRootNodeType = typename BoolTreeType::RootNodeType;
1108  using BoolNodeChainType = typename BoolRootNodeType::NodeChainType;
1109  using BoolInternalNodeType =
1110  typename hboost::mpl::at<BoolNodeChainType, hboost::mpl::int_<1>>::type;
1111 
1112  /////
1113 
1114  // Clamp the isovalue to the level set's background value minus epsilon.
1115  // (In a valid narrow-band level set, all voxels, including background voxels,
1116  // have values less than or equal to the background value, so an isovalue
1117  // greater than or equal to the background value would produce a mask with
1118  // effectively infinite extent.)
1119  iso = std::min(iso,
1120  static_cast<ValueType>(tree.background() - math::Tolerance<ValueType>::value()));
1121 
1122  size_t numLeafNodes = 0, numInternalNodes = 0;
1123 
1124  std::vector<const LeafNodeType*> nodes;
1125  std::vector<size_t> leafnodeCount;
1126 
1127  {
1128  // compute the prefix sum of the leafnode count in each internal node.
1129  std::vector<const InternalNodeType*> internalNodes;
1130  tree.getNodes(internalNodes);
1131 
1132  numInternalNodes = internalNodes.size();
1133 
1134  leafnodeCount.push_back(0);
1135  for (size_t n = 0; n < numInternalNodes; ++n) {
1136  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
1137  }
1138 
1139  numLeafNodes = leafnodeCount.back();
1140 
1141  // extract all leafnodes
1142  nodes.reserve(numLeafNodes);
1143 
1144  for (size_t n = 0; n < numInternalNodes; ++n) {
1145  internalNodes[n]->getNodes(nodes);
1146  }
1147  }
1148 
1149  // create mask leafnodes
1150  std::unique_ptr<BoolLeafNodeType*[]> maskNodes(new BoolLeafNodeType*[numLeafNodes]);
1151 
1152  tbb::parallel_for(tbb::blocked_range<size_t>(0, numLeafNodes),
1153  MaskInteriorVoxels<LeafNodeType>(iso, &nodes[0], maskNodes.get()));
1154 
1155 
1156  // create mask grid
1157  typename BoolTreeType::Ptr maskTree(new BoolTreeType(false));
1158 
1159  PopulateTree<BoolTreeType> populate(*maskTree, maskNodes.get(), &leafnodeCount[0], false);
1160  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
1161 
1162 
1163  // evaluate tile values
1164  std::vector<BoolInternalNodeType*> internalMaskNodes;
1165  maskTree->getNodes(internalMaskNodes);
1166 
1167  tbb::parallel_for(tbb::blocked_range<size_t>(0, internalMaskNodes.size()),
1168  MaskInteriorTiles<TreeType, BoolInternalNodeType>(iso, tree, &internalMaskNodes[0]));
1169 
1171 
1172  typename BoolTreeType::ValueAllIter it(*maskTree);
1173  it.setMaxDepth(BoolTreeType::ValueAllIter::LEAF_DEPTH - 2);
1174 
1175  for ( ; it; ++it) {
1176  if (acc.getValue(it.getCoord()) < iso) {
1177  it.setValue(true);
1178  it.setActiveState(true);
1179  }
1180  }
1181 
1182  return maskTree;
1183 } // computeInteriorMask()
1184 
1185 
1186 template<typename InputTreeType>
1188 {
1189  using InputValueType = typename InputTreeType::ValueType;
1190  using InputLeafNodeType = typename InputTreeType::LeafNodeType;
1191  using BoolTreeType = typename InputTreeType::template ValueConverter<bool>::Type;
1192  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1193 
1195  const InputTreeType& inputTree,
1196  const std::vector<const InputLeafNodeType*>& inputLeafNodes,
1197  BoolTreeType& maskTree,
1198  InputValueType iso)
1199  : mInputAccessor(inputTree)
1200  , mInputNodes(!inputLeafNodes.empty() ? &inputLeafNodes.front() : nullptr)
1201  , mMaskTree(false)
1202  , mMaskAccessor(maskTree)
1203  , mIsovalue(iso)
1204  {
1205  }
1206 
1208  : mInputAccessor(rhs.mInputAccessor.tree())
1209  , mInputNodes(rhs.mInputNodes)
1210  , mMaskTree(false)
1211  , mMaskAccessor(mMaskTree)
1212  , mIsovalue(rhs.mIsovalue)
1213  {
1214  }
1215 
1216  void operator()(const tbb::blocked_range<size_t>& range) {
1217 
1218  const InputValueType iso = mIsovalue;
1219  Coord ijk(0, 0, 0);
1220 
1221  BoolLeafNodeType* maskNodePt = nullptr;
1222 
1223  for (size_t n = range.begin(); mInputNodes && (n != range.end()); ++n) {
1224 
1225  const InputLeafNodeType& node = *mInputNodes[n];
1226 
1227  if (!maskNodePt) maskNodePt = new BoolLeafNodeType(node.origin(), false);
1228  else maskNodePt->setOrigin(node.origin());
1229 
1230  bool collectedData = false;
1231 
1232  for (typename InputLeafNodeType::ValueOnCIter it = node.cbeginValueOn(); it; ++it) {
1233 
1234  bool isUnder = *it < iso;
1235 
1236  ijk = it.getCoord();
1237 
1238  ++ijk[2];
1239  bool signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +z edge
1240  --ijk[2];
1241 
1242  if (!signChange) {
1243  --ijk[2];
1244  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -z edge
1245  ++ijk[2];
1246  }
1247 
1248  if (!signChange) {
1249  ++ijk[1];
1250  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +y edge
1251  --ijk[1];
1252  }
1253 
1254  if (!signChange) {
1255  --ijk[1];
1256  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -y edge
1257  ++ijk[1];
1258  }
1259 
1260  if (!signChange) {
1261  ++ijk[0];
1262  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // +x edge
1263  --ijk[0];
1264  }
1265 
1266  if (!signChange) {
1267  --ijk[0];
1268  signChange = isUnder != (mInputAccessor.getValue(ijk) < iso); // -x edge
1269  ++ijk[0];
1270  }
1271 
1272  if (signChange) {
1273  collectedData = true;
1274  maskNodePt->setValueOn(it.pos(), true);
1275  }
1276  }
1277 
1278  if (collectedData) {
1279  mMaskAccessor.addLeaf(maskNodePt);
1280  maskNodePt = nullptr;
1281  }
1282  }
1283 
1284  if (maskNodePt) delete maskNodePt;
1285  }
1286 
1288  mMaskAccessor.tree().merge(rhs.mMaskAccessor.tree());
1289  }
1290 
1291 private:
1293  InputLeafNodeType const * const * const mInputNodes;
1294 
1295  BoolTreeType mMaskTree;
1296  tree::ValueAccessor<BoolTreeType> mMaskAccessor;
1297 
1298  InputValueType mIsovalue;
1299 }; // MaskIsovalueCrossingVoxels
1300 
1301 
1302 ////////////////////////////////////////
1303 
1304 
1305 template<typename NodeType>
1307 {
1309  using NodeMaskType = typename NodeType::NodeMaskType;
1310 
1311  NodeMaskSegment() : connections(), mask(false), origin(0,0,0), visited(false) {}
1312 
1313  std::vector<NodeMaskSegment*> connections;
1315  Coord origin;
1316  bool visited;
1317 }; // struct NodeMaskSegment
1318 
1319 
1320 template<typename NodeType>
1321 inline void
1322 nodeMaskSegmentation(const NodeType& node,
1323  std::vector<typename NodeMaskSegment<NodeType>::Ptr>& segments)
1324 {
1325  using NodeMaskType = typename NodeType::NodeMaskType;
1326  using NodeMaskSegmentType = NodeMaskSegment<NodeType>;
1327  using NodeMaskSegmentTypePtr = typename NodeMaskSegmentType::Ptr;
1328 
1329  NodeMaskType nodeMask(node.getValueMask());
1330  std::deque<Index> indexList;
1331 
1332  while (!nodeMask.isOff()) {
1333 
1334  NodeMaskSegmentTypePtr segment(new NodeMaskSegmentType());
1335  segment->origin = node.origin();
1336 
1337  NodeMaskType& mask = segment->mask;
1338 
1339  indexList.push_back(nodeMask.findFirstOn());
1340  nodeMask.setOff(indexList.back()); // mark as visited
1341  Coord ijk(0, 0, 0);
1342 
1343  while (!indexList.empty()) {
1344 
1345  const Index pos = indexList.back();
1346  indexList.pop_back();
1347 
1348  if (mask.isOn(pos)) continue;
1349  mask.setOn(pos);
1350 
1351  ijk = NodeType::offsetToLocalCoord(pos);
1352 
1353  Index npos = pos - 1;
1354  if (ijk[2] != 0 && nodeMask.isOn(npos)) {
1355  nodeMask.setOff(npos);
1356  indexList.push_back(npos);
1357  }
1358 
1359  npos = pos + 1;
1360  if (ijk[2] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1361  nodeMask.setOff(npos);
1362  indexList.push_back(npos);
1363  }
1364 
1365  npos = pos - NodeType::DIM;
1366  if (ijk[1] != 0 && nodeMask.isOn(npos)) {
1367  nodeMask.setOff(npos);
1368  indexList.push_back(npos);
1369  }
1370 
1371  npos = pos + NodeType::DIM;
1372  if (ijk[1] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1373  nodeMask.setOff(npos);
1374  indexList.push_back(npos);
1375  }
1376 
1377  npos = pos - NodeType::DIM * NodeType::DIM;
1378  if (ijk[0] != 0 && nodeMask.isOn(npos)) {
1379  nodeMask.setOff(npos);
1380  indexList.push_back(npos);
1381  }
1382 
1383  npos = pos + NodeType::DIM * NodeType::DIM;
1384  if (ijk[0] != (NodeType::DIM - 1) && nodeMask.isOn(npos)) {
1385  nodeMask.setOff(npos);
1386  indexList.push_back(npos);
1387  }
1388 
1389  }
1390 
1391  segments.push_back(segment);
1392  }
1393 }
1394 
1395 
1396 template<typename NodeType>
1398 {
1401  using NodeMaskSegmentVector = typename std::vector<NodeMaskSegmentTypePtr>;
1402 
1403  SegmentNodeMask(std::vector<NodeType*>& nodes, NodeMaskSegmentVector* nodeMaskArray)
1404  : mNodes(!nodes.empty() ? &nodes.front() : nullptr)
1405  , mNodeMaskArray(nodeMaskArray)
1406  {
1407  }
1408 
1409  void operator()(const tbb::blocked_range<size_t>& range) const {
1410  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1411  NodeType& node = *mNodes[n];
1413 
1414  // hack origin data to store array offset
1415  Coord& origin = const_cast<Coord&>(node.origin());
1416  origin[0] = static_cast<int>(n);
1417  }
1418  }
1419 
1420  NodeType * const * const mNodes;
1422 }; // struct SegmentNodeMask
1423 
1424 
1425 template<typename TreeType, typename NodeType>
1427 {
1428  using NodeMaskType = typename NodeType::NodeMaskType;
1431  using NodeMaskSegmentVector = typename std::vector<NodeMaskSegmentTypePtr>;
1432 
1433  ConnectNodeMaskSegments(const TreeType& tree, NodeMaskSegmentVector* nodeMaskArray)
1434  : mTree(&tree)
1435  , mNodeMaskArray(nodeMaskArray)
1436  {
1437  }
1438 
1439  void operator()(const tbb::blocked_range<size_t>& range) const {
1440 
1442 
1443  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1444 
1445  NodeMaskSegmentVector& segments = mNodeMaskArray[n];
1446  if (segments.empty()) continue;
1447 
1448  std::vector<std::set<NodeMaskSegmentType*> > connections(segments.size());
1449 
1450  Coord ijk = segments[0]->origin;
1451 
1452  const NodeType* node = acc.template probeConstNode<NodeType>(ijk);
1453  if (!node) continue;
1454 
1455  // get neighbour nodes
1456 
1457  ijk[2] += NodeType::DIM;
1458  const NodeType* nodeZUp = acc.template probeConstNode<NodeType>(ijk);
1459  ijk[2] -= (NodeType::DIM + NodeType::DIM);
1460  const NodeType* nodeZDown = acc.template probeConstNode<NodeType>(ijk);
1461  ijk[2] += NodeType::DIM;
1462 
1463  ijk[1] += NodeType::DIM;
1464  const NodeType* nodeYUp = acc.template probeConstNode<NodeType>(ijk);
1465  ijk[1] -= (NodeType::DIM + NodeType::DIM);
1466  const NodeType* nodeYDown = acc.template probeConstNode<NodeType>(ijk);
1467  ijk[1] += NodeType::DIM;
1468 
1469  ijk[0] += NodeType::DIM;
1470  const NodeType* nodeXUp = acc.template probeConstNode<NodeType>(ijk);
1471  ijk[0] -= (NodeType::DIM + NodeType::DIM);
1472  const NodeType* nodeXDown = acc.template probeConstNode<NodeType>(ijk);
1473  ijk[0] += NodeType::DIM;
1474 
1475  const Index startPos = node->getValueMask().findFirstOn();
1476  for (Index pos = startPos; pos < NodeMaskType::SIZE; ++pos) {
1477 
1478  if (!node->isValueOn(pos)) continue;
1479 
1480  ijk = NodeType::offsetToLocalCoord(pos);
1481 
1482 #ifdef _MSC_FULL_VER
1483  #if _MSC_FULL_VER >= 190000000 && _MSC_FULL_VER < 190024210
1484  // Visual Studio 2015 had a codegen bug that wasn't fixed until Update 3
1485  volatile Index npos = 0;
1486  #else
1487  Index npos = 0;
1488  #endif
1489 #else
1490  Index npos = 0;
1491 #endif
1492 
1493  if (ijk[2] == 0) {
1494  npos = pos + (NodeType::DIM - 1);
1495  if (nodeZDown && nodeZDown->isValueOn(npos)) {
1496  NodeMaskSegmentType* nsegment =
1497  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZDown)], npos);
1498  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1499  connections[idx].insert(nsegment);
1500  }
1501  } else if (ijk[2] == (NodeType::DIM - 1)) {
1502  npos = pos - (NodeType::DIM - 1);
1503  if (nodeZUp && nodeZUp->isValueOn(npos)) {
1504  NodeMaskSegmentType* nsegment =
1505  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeZUp)], npos);
1506  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1507  connections[idx].insert(nsegment);
1508  }
1509  }
1510 
1511  if (ijk[1] == 0) {
1512  npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1513  if (nodeYDown && nodeYDown->isValueOn(npos)) {
1514  NodeMaskSegmentType* nsegment =
1515  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYDown)], npos);
1516  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1517  connections[idx].insert(nsegment);
1518  }
1519  } else if (ijk[1] == (NodeType::DIM - 1)) {
1520  npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1521  if (nodeYUp && nodeYUp->isValueOn(npos)) {
1522  NodeMaskSegmentType* nsegment =
1523  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeYUp)], npos);
1524  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1525  connections[idx].insert(nsegment);
1526  }
1527  }
1528 
1529  if (ijk[0] == 0) {
1530  npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1531  if (nodeXDown && nodeXDown->isValueOn(npos)) {
1532  NodeMaskSegmentType* nsegment =
1533  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXDown)], npos);
1534  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1535  connections[idx].insert(nsegment);
1536  }
1537  } else if (ijk[0] == (NodeType::DIM - 1)) {
1538  npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1539  if (nodeXUp && nodeXUp->isValueOn(npos)) {
1540  NodeMaskSegmentType* nsegment =
1541  findNodeMaskSegment(mNodeMaskArray[getNodeOffset(*nodeXUp)], npos);
1542  const Index idx = findNodeMaskSegmentIndex(segments, pos);
1543  connections[idx].insert(nsegment);
1544  }
1545  }
1546  }
1547 
1548  for (size_t i = 0, I = connections.size(); i < I; ++i) {
1549 
1550  typename std::set<NodeMaskSegmentType*>::iterator
1551  it = connections[i].begin(), end = connections[i].end();
1552 
1553  std::vector<NodeMaskSegmentType*>& segmentConnections = segments[i]->connections;
1554  segmentConnections.reserve(connections.size());
1555  for (; it != end; ++it) {
1556  segmentConnections.push_back(*it);
1557  }
1558  }
1559  } // end range loop
1560  }
1561 
1562 private:
1563 
1564  static inline size_t getNodeOffset(const NodeType& node) {
1565  return static_cast<size_t>(node.origin()[0]);
1566  }
1567 
1568  static inline NodeMaskSegmentType*
1569  findNodeMaskSegment(NodeMaskSegmentVector& segments, Index pos)
1570  {
1571  NodeMaskSegmentType* segment = nullptr;
1572 
1573  for (size_t n = 0, N = segments.size(); n < N; ++n) {
1574  if (segments[n]->mask.isOn(pos)) {
1575  segment = segments[n].get();
1576  break;
1577  }
1578  }
1579 
1580  return segment;
1581  }
1582 
1583  static inline Index
1584  findNodeMaskSegmentIndex(NodeMaskSegmentVector& segments, Index pos)
1585  {
1586  for (Index n = 0, N = Index(segments.size()); n < N; ++n) {
1587  if (segments[n]->mask.isOn(pos)) return n;
1588  }
1589  return Index(-1);
1590  }
1591 
1592  TreeType const * const mTree;
1593  NodeMaskSegmentVector * const mNodeMaskArray;
1594 }; // struct ConnectNodeMaskSegments
1595 
1596 
1597 template<typename TreeType>
1599 {
1600  using LeafNodeType = typename TreeType::LeafNodeType;
1601  using TreeTypePtr = typename TreeType::Ptr;
1603 
1604  MaskSegmentGroup(const std::vector<NodeMaskSegmentType*>& segments)
1605  : mSegments(!segments.empty() ? &segments.front() : nullptr)
1606  , mTree(new TreeType(false))
1607  {
1608  }
1609 
1610  MaskSegmentGroup(const MaskSegmentGroup& rhs, tbb::split)
1611  : mSegments(rhs.mSegments)
1612  , mTree(new TreeType(false))
1613  {
1614  }
1615 
1616  TreeTypePtr& mask() { return mTree; }
1617 
1618  void join(MaskSegmentGroup& rhs) { mTree->merge(*rhs.mTree); }
1619 
1620  void operator()(const tbb::blocked_range<size_t>& range) {
1621 
1622  tree::ValueAccessor<TreeType> acc(*mTree);
1623 
1624  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1625  NodeMaskSegmentType& segment = *mSegments[n];
1626  LeafNodeType* node = acc.touchLeaf(segment.origin);
1627  node->getValueMask() |= segment.mask;
1628  }
1629  }
1630 
1631 private:
1632  NodeMaskSegmentType * const * const mSegments;
1633  TreeTypePtr mTree;
1634 }; // struct MaskSegmentGroup
1635 
1636 
1637 ////////////////////////////////////////
1638 
1639 
1640 template<typename TreeType>
1642 {
1643  using ValueType = typename TreeType::ValueType;
1644  using LeafNodeType = typename TreeType::LeafNodeType;
1645  using NodeMaskType = typename LeafNodeType::NodeMaskType;
1646 
1647  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1648  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1649 
1650  /////
1651 
1652  ExpandLeafNodeRegion(const TreeType& distTree, BoolTreeType& maskTree,
1653  std::vector<BoolLeafNodeType*>& maskNodes)
1654  : mDistTree(&distTree)
1655  , mMaskTree(&maskTree)
1656  , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
1657  , mNewMaskTree(false)
1658  {
1659  }
1660 
1662  : mDistTree(rhs.mDistTree)
1663  , mMaskTree(rhs.mMaskTree)
1664  , mMaskNodes(rhs.mMaskNodes)
1665  , mNewMaskTree(false)
1666  {
1667  }
1668 
1669  BoolTreeType& newMaskTree() { return mNewMaskTree; }
1670 
1671  void join(ExpandLeafNodeRegion& rhs) { mNewMaskTree.merge(rhs.mNewMaskTree); }
1672 
1673  void operator()(const tbb::blocked_range<size_t>& range) {
1674 
1675  using NodeType = LeafNodeType;
1676 
1677  tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1678  tree::ValueAccessor<const BoolTreeType> maskAcc(*mMaskTree);
1679  tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
1680 
1681  NodeMaskType maskZUp, maskZDown, maskYUp, maskYDown, maskXUp, maskXDown;
1682 
1683  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1684 
1685  BoolLeafNodeType& maskNode = *mMaskNodes[n];
1686  if (maskNode.isEmpty()) continue;
1687 
1688  Coord ijk = maskNode.origin(), nijk;
1689 
1690  const LeafNodeType* distNode = distAcc.probeConstLeaf(ijk);
1691  if (!distNode) continue;
1692 
1693  const ValueType *dataZUp = nullptr, *dataZDown = nullptr,
1694  *dataYUp = nullptr, *dataYDown = nullptr,
1695  *dataXUp = nullptr, *dataXDown = nullptr;
1696 
1697  ijk[2] += NodeType::DIM;
1698  getData(ijk, distAcc, maskAcc, maskZUp, dataZUp);
1699  ijk[2] -= (NodeType::DIM + NodeType::DIM);
1700  getData(ijk, distAcc, maskAcc, maskZDown, dataZDown);
1701  ijk[2] += NodeType::DIM;
1702 
1703  ijk[1] += NodeType::DIM;
1704  getData(ijk, distAcc, maskAcc, maskYUp, dataYUp);
1705  ijk[1] -= (NodeType::DIM + NodeType::DIM);
1706  getData(ijk, distAcc, maskAcc, maskYDown, dataYDown);
1707  ijk[1] += NodeType::DIM;
1708 
1709  ijk[0] += NodeType::DIM;
1710  getData(ijk, distAcc, maskAcc, maskXUp, dataXUp);
1711  ijk[0] -= (NodeType::DIM + NodeType::DIM);
1712  getData(ijk, distAcc, maskAcc, maskXDown, dataXDown);
1713  ijk[0] += NodeType::DIM;
1714 
1715  for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
1716 
1717  const Index pos = it.pos();
1718  const ValueType val = std::abs(distNode->getValue(pos));
1719 
1720  ijk = BoolLeafNodeType::offsetToLocalCoord(pos);
1721  nijk = ijk + maskNode.origin();
1722 
1723  if (dataZUp && ijk[2] == (BoolLeafNodeType::DIM - 1)) {
1724  const Index npos = pos - (NodeType::DIM - 1);
1725  if (maskZUp.isOn(npos) && std::abs(dataZUp[npos]) > val) {
1726  newMaskAcc.setValueOn(nijk.offsetBy(0, 0, 1));
1727  }
1728  } else if (dataZDown && ijk[2] == 0) {
1729  const Index npos = pos + (NodeType::DIM - 1);
1730  if (maskZDown.isOn(npos) && std::abs(dataZDown[npos]) > val) {
1731  newMaskAcc.setValueOn(nijk.offsetBy(0, 0, -1));
1732  }
1733  }
1734 
1735  if (dataYUp && ijk[1] == (BoolLeafNodeType::DIM - 1)) {
1736  const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM;
1737  if (maskYUp.isOn(npos) && std::abs(dataYUp[npos]) > val) {
1738  newMaskAcc.setValueOn(nijk.offsetBy(0, 1, 0));
1739  }
1740  } else if (dataYDown && ijk[1] == 0) {
1741  const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM;
1742  if (maskYDown.isOn(npos) && std::abs(dataYDown[npos]) > val) {
1743  newMaskAcc.setValueOn(nijk.offsetBy(0, -1, 0));
1744  }
1745  }
1746 
1747  if (dataXUp && ijk[0] == (BoolLeafNodeType::DIM - 1)) {
1748  const Index npos = pos - (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1749  if (maskXUp.isOn(npos) && std::abs(dataXUp[npos]) > val) {
1750  newMaskAcc.setValueOn(nijk.offsetBy(1, 0, 0));
1751  }
1752  } else if (dataXDown && ijk[0] == 0) {
1753  const Index npos = pos + (NodeType::DIM - 1) * NodeType::DIM * NodeType::DIM;
1754  if (maskXDown.isOn(npos) && std::abs(dataXDown[npos]) > val) {
1755  newMaskAcc.setValueOn(nijk.offsetBy(-1, 0, 0));
1756  }
1757  }
1758 
1759  } // end value on loop
1760  } // end range loop
1761  }
1762 
1763 private:
1764 
1765  static inline void
1766  getData(const Coord& ijk, tree::ValueAccessor<const TreeType>& distAcc,
1768  const ValueType*& data)
1769  {
1770  const LeafNodeType* node = distAcc.probeConstLeaf(ijk);
1771  if (node) {
1772  data = node->buffer().data();
1773  mask = node->getValueMask();
1774  const BoolLeafNodeType* maskNodePt = maskAcc.probeConstLeaf(ijk);
1775  if (maskNodePt) mask -= maskNodePt->getValueMask();
1776  }
1777  }
1778 
1779  TreeType const * const mDistTree;
1780  BoolTreeType * const mMaskTree;
1781  BoolLeafNodeType ** const mMaskNodes;
1782 
1783  BoolTreeType mNewMaskTree;
1784 }; // struct ExpandLeafNodeRegion
1785 
1786 
1787 template<typename TreeType>
1789 {
1790  using ValueType = typename TreeType::ValueType;
1791  using LeafNodeType = typename TreeType::LeafNodeType;
1792  using NodeMaskType = typename LeafNodeType::NodeMaskType;
1794 
1795  FillLeafNodeVoxels(const TreeType& tree, std::vector<BoolLeafNodeType*>& maskNodes)
1796  : mTree(&tree), mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
1797  {
1798  }
1799 
1800  void operator()(const tbb::blocked_range<size_t>& range) const {
1801 
1803 
1804  std::vector<Index> indexList;
1805  indexList.reserve(NodeMaskType::SIZE);
1806 
1807  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1808 
1809  BoolLeafNodeType& maskNode = *mMaskNodes[n];
1810 
1811  const LeafNodeType * distNode = distAcc.probeConstLeaf(maskNode.origin());
1812  if (!distNode) continue;
1813 
1814  NodeMaskType mask(distNode->getValueMask());
1815  NodeMaskType& narrowbandMask = maskNode.getValueMask();
1816 
1817  for (Index pos = narrowbandMask.findFirstOn(); pos < NodeMaskType::SIZE; ++pos) {
1818  if (narrowbandMask.isOn(pos)) indexList.push_back(pos);
1819  }
1820 
1821  mask -= narrowbandMask; // bitwise difference
1822  narrowbandMask.setOff();
1823 
1824  const ValueType* data = distNode->buffer().data();
1825  Coord ijk(0, 0, 0);
1826 
1827  while (!indexList.empty()) {
1828 
1829  const Index pos = indexList.back();
1830  indexList.pop_back();
1831 
1832  if (narrowbandMask.isOn(pos)) continue;
1833  narrowbandMask.setOn(pos);
1834 
1835  const ValueType dist = std::abs(data[pos]);
1836 
1837  ijk = LeafNodeType::offsetToLocalCoord(pos);
1838 
1839  Index npos = pos - 1;
1840  if (ijk[2] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1841  mask.setOff(npos);
1842  indexList.push_back(npos);
1843  }
1844 
1845  npos = pos + 1;
1846  if ((ijk[2] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1847  && std::abs(data[npos]) > dist)
1848  {
1849  mask.setOff(npos);
1850  indexList.push_back(npos);
1851  }
1852 
1853  npos = pos - LeafNodeType::DIM;
1854  if (ijk[1] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1855  mask.setOff(npos);
1856  indexList.push_back(npos);
1857  }
1858 
1859  npos = pos + LeafNodeType::DIM;
1860  if ((ijk[1] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1861  && std::abs(data[npos]) > dist)
1862  {
1863  mask.setOff(npos);
1864  indexList.push_back(npos);
1865  }
1866 
1867  npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1868  if (ijk[0] != 0 && mask.isOn(npos) && std::abs(data[npos]) > dist) {
1869  mask.setOff(npos);
1870  indexList.push_back(npos);
1871  }
1872 
1873  npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1874  if ((ijk[0] != (LeafNodeType::DIM - 1)) && mask.isOn(npos)
1875  && std::abs(data[npos]) > dist)
1876  {
1877  mask.setOff(npos);
1878  indexList.push_back(npos);
1879  }
1880  } // end flood fill loop
1881  } // end range loop
1882  }
1883 
1884  TreeType const * const mTree;
1886 }; // FillLeafNodeVoxels
1887 
1888 
1889 template<typename TreeType>
1891 {
1892  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
1893  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
1894  using BoolTreeTypePtr = typename BoolTreeType::Ptr;
1895 
1896  ExpandNarrowbandMask(const TreeType& tree, std::vector<BoolTreeTypePtr>& segments)
1897  : mTree(&tree), mSegments(!segments.empty() ? &segments.front() : nullptr)
1898  {
1899  }
1900 
1901  void operator()(const tbb::blocked_range<size_t>& range) const {
1902 
1903  const TreeType& distTree = *mTree;
1904  std::vector<BoolLeafNodeType*> nodes;
1905 
1906  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1907 
1908  BoolTreeType& narrowBandMask = *mSegments[n];
1909 
1910  BoolTreeType candidateMask(narrowBandMask, false, TopologyCopy());
1911 
1912  while (true) {
1913 
1914  nodes.clear();
1915  candidateMask.getNodes(nodes);
1916  if (nodes.empty()) break;
1917 
1918  const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
1919 
1920  tbb::parallel_for(nodeRange, FillLeafNodeVoxels<TreeType>(distTree, nodes));
1921 
1922  narrowBandMask.topologyUnion(candidateMask);
1923 
1924  ExpandLeafNodeRegion<TreeType> op(distTree, narrowBandMask, nodes);
1925  tbb::parallel_reduce(nodeRange, op);
1926 
1927  if (op.newMaskTree().empty()) break;
1928 
1929  candidateMask.clear();
1930  candidateMask.merge(op.newMaskTree());
1931  } // end expand loop
1932  } // end range loop
1933  }
1934 
1935  TreeType const * const mTree;
1937 }; // ExpandNarrowbandMask
1938 
1939 
1940 template<typename TreeType>
1942 {
1943  using TreeTypePtr = typename TreeType::Ptr;
1944  using ValueType = typename TreeType::ValueType;
1945  using LeafNodeType = typename TreeType::LeafNodeType;
1946  using RootNodeType = typename TreeType::RootNodeType;
1947  using NodeChainType = typename RootNodeType::NodeChainType;
1948  using InternalNodeType = typename hboost::mpl::at<NodeChainType, hboost::mpl::int_<1> >::type;
1949 
1950  FloodFillSign(const TreeType& tree, std::vector<TreeTypePtr>& segments)
1951  : mTree(&tree)
1952  , mSegments(!segments.empty() ? &segments.front() : nullptr)
1953  , mMinValue(ValueType(0.0))
1954  {
1956 
1957  {
1958  std::vector<const InternalNodeType*> nodes;
1959  tree.getNodes(nodes);
1960 
1961  if (!nodes.empty()) {
1962  FindMinTileValue<InternalNodeType> minOp(&nodes[0]);
1963  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1964  minSDFValue = std::min(minSDFValue, minOp.minValue);
1965  }
1966  }
1967 
1968  if (minSDFValue > ValueType(0.0)) {
1969  std::vector<const LeafNodeType*> nodes;
1970  tree.getNodes(nodes);
1971  if (!nodes.empty()) {
1972  FindMinVoxelValue<LeafNodeType> minOp(&nodes[0]);
1973  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
1974  minSDFValue = std::min(minSDFValue, minOp.minValue);
1975  }
1976  }
1977 
1978  mMinValue = minSDFValue;
1979  }
1980 
1981  void operator()(const tbb::blocked_range<size_t>& range) const {
1982  const ValueType interiorValue = -std::abs(mMinValue);
1983  const ValueType exteriorValue = std::abs(mTree->background());
1984  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1985  tools::signedFloodFillWithValues(*mSegments[n], exteriorValue, interiorValue);
1986  }
1987  }
1988 
1989 private:
1990 
1991  TreeType const * const mTree;
1992  TreeTypePtr * const mSegments;
1993  ValueType mMinValue;
1994 }; // FloodFillSign
1995 
1996 
1997 template<typename TreeType>
1999 {
2000  using TreeTypePtr = typename TreeType::Ptr;
2001  using ValueType = typename TreeType::ValueType;
2002  using LeafNodeType = typename TreeType::LeafNodeType;
2003 
2004  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2005  using BoolTreeTypePtr = typename BoolTreeType::Ptr;
2006  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2007 
2008  MaskedCopy(const TreeType& tree, std::vector<TreeTypePtr>& segments,
2009  std::vector<BoolTreeTypePtr>& masks)
2010  : mTree(&tree)
2011  , mSegments(!segments.empty() ? &segments.front() : nullptr)
2012  , mMasks(!masks.empty() ? &masks.front() : nullptr)
2013  {
2014  }
2015 
2016  void operator()(const tbb::blocked_range<size_t>& range) const {
2017 
2018  std::vector<const BoolLeafNodeType*> nodes;
2019 
2020  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2021 
2022  const BoolTreeType& mask = *mMasks[n];
2023 
2024  nodes.clear();
2025  mask.getNodes(nodes);
2026 
2027  Copy op(*mTree, nodes);
2028  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2029  mSegments[n] = op.outputTree();
2030  }
2031  }
2032 
2033 private:
2034 
2035  struct Copy {
2036  Copy(const TreeType& inputTree, std::vector<const BoolLeafNodeType*>& maskNodes)
2037  : mInputTree(&inputTree)
2038  , mMaskNodes(!maskNodes.empty() ? &maskNodes.front() : nullptr)
2039  , mOutputTreePtr(new TreeType(inputTree.background()))
2040  {
2041  }
2042 
2043  Copy(const Copy& rhs, tbb::split)
2044  : mInputTree(rhs.mInputTree)
2045  , mMaskNodes(rhs.mMaskNodes)
2046  , mOutputTreePtr(new TreeType(mInputTree->background()))
2047  {
2048  }
2049 
2050  TreeTypePtr& outputTree() { return mOutputTreePtr; }
2051 
2052  void join(Copy& rhs) { mOutputTreePtr->merge(*rhs.mOutputTreePtr); }
2053 
2054  void operator()(const tbb::blocked_range<size_t>& range) {
2055 
2056  tree::ValueAccessor<const TreeType> inputAcc(*mInputTree);
2057  tree::ValueAccessor<TreeType> outputAcc(*mOutputTreePtr);
2058 
2059  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2060 
2061  const BoolLeafNodeType& maskNode = *mMaskNodes[n];
2062  if (maskNode.isEmpty()) continue;
2063 
2064  const Coord& ijk = maskNode.origin();
2065 
2066  const LeafNodeType* inputNode = inputAcc.probeConstLeaf(ijk);
2067  if (inputNode) {
2068 
2069  LeafNodeType* outputNode = outputAcc.touchLeaf(ijk);
2070 
2071  for (typename BoolLeafNodeType::ValueOnCIter it = maskNode.cbeginValueOn();
2072  it; ++it)
2073  {
2074  const Index idx = it.pos();
2075  outputNode->setValueOn(idx, inputNode->getValue(idx));
2076  }
2077  } else {
2078  const int valueDepth = inputAcc.getValueDepth(ijk);
2079  if (valueDepth >= 0) {
2080  outputAcc.addTile(TreeType::RootNodeType::LEVEL - valueDepth,
2081  ijk, inputAcc.getValue(ijk), true);
2082  }
2083  }
2084  }
2085  }
2086 
2087  private:
2088  TreeType const * const mInputTree;
2089  BoolLeafNodeType const * const * const mMaskNodes;
2090  TreeTypePtr mOutputTreePtr;
2091  }; // struct Copy
2092 
2093  TreeType const * const mTree;
2094  TreeTypePtr * const mSegments;
2095  BoolTreeTypePtr * const mMasks;
2096 }; // MaskedCopy
2097 
2098 
2099 ////////////////////////////////////////
2100 
2101 
2102 template<typename VolumePtrType>
2104 {
2105  ComputeActiveVoxelCount(std::vector<VolumePtrType>& segments, size_t *countArray)
2106  : mSegments(!segments.empty() ? &segments.front() : nullptr)
2107  , mCountArray(countArray)
2108  {
2109  }
2110 
2111  void operator()(const tbb::blocked_range<size_t>& range) const {
2112  for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2113  mCountArray[n] = mSegments[n]->activeVoxelCount();
2114  }
2115  }
2116 
2117  VolumePtrType * const mSegments;
2118  size_t * const mCountArray;
2119 };
2120 
2121 
2123 {
2124  GreaterCount(const size_t *countArray) : mCountArray(countArray) {}
2125 
2126  inline bool operator() (const size_t& lhs, const size_t& rhs) const
2127  {
2128  return (mCountArray[lhs] > mCountArray[rhs]);
2129  }
2130 
2131  size_t const * const mCountArray;
2132 };
2133 
2134 ////////////////////////////////////////
2135 
2136 
2137 template<typename TreeType>
2139 {
2140  using TreeTypePtr = typename TreeType::Ptr;
2141  using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2142 
2143  static BoolTreePtrType constructMask(const TreeType&, BoolTreePtrType& maskTree)
2144  { return maskTree; }
2145  static TreeTypePtr construct(const TreeType&, TreeTypePtr& tree) { return tree; }
2146 };
2147 
2148 
2149 template<typename TreeType>
2150 struct GridOrTreeConstructor<Grid<TreeType> >
2151 {
2154  using TreeTypePtr = typename TreeType::Ptr;
2155 
2156  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2157  using BoolTreePtrType = typename BoolTreeType::Ptr;
2160 
2161  static BoolGridPtrType constructMask(const GridType& grid, BoolTreePtrType& maskTree) {
2162  BoolGridPtrType maskGrid(BoolGridType::create(maskTree));
2163  maskGrid->setTransform(grid.transform().copy());
2164  return maskGrid;
2165  }
2166 
2167  static GridTypePtr construct(const GridType& grid, TreeTypePtr& maskTree) {
2168  GridTypePtr maskGrid(GridType::create(maskTree));
2169  maskGrid->setTransform(grid.transform().copy());
2170  maskGrid->insertMeta(grid);
2171  return maskGrid;
2172  }
2173 };
2174 
2175 
2176 } // namespace level_set_util_internal
2177 
2178 
2179 ////////////////////////////////////////
2180 
2181 
2182 template <class GridType>
2183 inline void
2184 sdfToFogVolume(GridType& grid, typename GridType::ValueType cutoffDistance)
2185 {
2186  using ValueType = typename GridType::ValueType;
2187  using TreeType = typename GridType::TreeType;
2188  using LeafNodeType = typename TreeType::LeafNodeType;
2189  using RootNodeType = typename TreeType::RootNodeType;
2190  using NodeChainType = typename RootNodeType::NodeChainType;
2191  using InternalNodeType = typename hboost::mpl::at<NodeChainType, hboost::mpl::int_<1>>::type;
2192 
2193  //////////
2194 
2195  TreeType& tree = grid.tree();
2196 
2197  size_t numLeafNodes = 0, numInternalNodes = 0;
2198 
2199  std::vector<LeafNodeType*> nodes;
2200  std::vector<size_t> leafnodeCount;
2201 
2202  {
2203  // Compute the prefix sum of the leafnode count in each internal node.
2204  std::vector<InternalNodeType*> internalNodes;
2205  tree.getNodes(internalNodes);
2206 
2207  numInternalNodes = internalNodes.size();
2208 
2209  leafnodeCount.push_back(0);
2210  for (size_t n = 0; n < numInternalNodes; ++n) {
2211  leafnodeCount.push_back(leafnodeCount.back() + internalNodes[n]->leafCount());
2212  }
2213 
2214  numLeafNodes = leafnodeCount.back();
2215 
2216  // Steal all leafnodes (Removes them from the tree and transfers ownership.)
2217  nodes.reserve(numLeafNodes);
2218 
2219  for (size_t n = 0; n < numInternalNodes; ++n) {
2220  internalNodes[n]->stealNodes(nodes, tree.background(), false);
2221  }
2222 
2223  // Clamp cutoffDistance to min sdf value
2224  ValueType minSDFValue = std::numeric_limits<ValueType>::max();
2225 
2226  {
2228  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, internalNodes.size()), minOp);
2229  minSDFValue = std::min(minSDFValue, minOp.minValue);
2230  }
2231 
2232  if (minSDFValue > ValueType(0.0)) {
2234  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), minOp);
2235  minSDFValue = std::min(minSDFValue, minOp.minValue);
2236  }
2237 
2238  cutoffDistance = -std::abs(cutoffDistance);
2239  cutoffDistance = minSDFValue > cutoffDistance ? minSDFValue : cutoffDistance;
2240  }
2241 
2242  // Transform voxel values and delete leafnodes that are uniformly zero after the transformation.
2243  // (Positive values are set to zero with inactive state and negative values are remapped
2244  // from zero to one with active state.)
2245  tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
2247 
2248  // Populate a new tree with the remaining leafnodes
2249  typename TreeType::Ptr newTree(new TreeType(ValueType(0.0)));
2250 
2252  *newTree, &nodes[0], &leafnodeCount[0], 0);
2253  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, numInternalNodes), populate);
2254 
2255  // Transform tile values (Negative valued tiles are set to 1.0 with active state.)
2256  std::vector<InternalNodeType*> internalNodes;
2257  newTree->getNodes(internalNodes);
2258 
2259  tbb::parallel_for(tbb::blocked_range<size_t>(0, internalNodes.size()),
2261  tree, &internalNodes[0]));
2262 
2263  {
2265 
2266  typename TreeType::ValueAllIter it(*newTree);
2267  it.setMaxDepth(TreeType::ValueAllIter::LEAF_DEPTH - 2);
2268 
2269  for ( ; it; ++it) {
2270  if (acc.getValue(it.getCoord()) < ValueType(0.0)) {
2271  it.setValue(ValueType(1.0));
2272  it.setActiveState(true);
2273  }
2274  }
2275  }
2276 
2277  // Insert missing root level tiles. (The new tree is constructed from the remaining leafnodes
2278  // and will therefore not contain any root level tiles that may exist in the original tree.)
2279  {
2280  typename TreeType::ValueAllIter it(tree);
2281  it.setMaxDepth(TreeType::ValueAllIter::ROOT_DEPTH);
2282  for ( ; it; ++it) {
2283  if (it.getValue() < ValueType(0.0)) {
2284  newTree->addTile(TreeType::ValueAllIter::ROOT_LEVEL, it.getCoord(),
2285  ValueType(1.0), true);
2286  }
2287  }
2288  }
2289 
2290  grid.setTree(newTree);
2291  grid.setGridClass(GRID_FOG_VOLUME);
2292 }
2293 
2294 
2295 ////////////////////////////////////////
2296 
2297 
2298 template <class GridOrTreeType>
2299 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2300 sdfInteriorMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2301 {
2302  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2303  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2304 
2305  using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2306  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(tree, isovalue);
2307 
2309  volume, mask);
2310 }
2311 
2312 
2313 template<typename GridOrTreeType>
2314 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2315 extractEnclosedRegion(const GridOrTreeType& volume,
2316  typename GridOrTreeType::ValueType isovalue,
2318  fillMask)
2319 {
2320  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2321  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2322 
2323  using CharTreePtrType = typename TreeType::template ValueConverter<char>::Type::Ptr;
2324  CharTreePtrType regionMask = level_set_util_internal::computeEnclosedRegionMask(
2325  tree, isovalue, fillMask);
2326 
2327  using BoolTreePtrType = typename TreeType::template ValueConverter<bool>::Type::Ptr;
2328  BoolTreePtrType mask = level_set_util_internal::computeInteriorMask(*regionMask, 0);
2329 
2331  volume, mask);
2332 }
2333 
2334 
2335 ////////////////////////////////////////
2336 
2337 
2338 template<typename GridOrTreeType>
2339 inline typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr
2340 extractIsosurfaceMask(const GridOrTreeType& volume, typename GridOrTreeType::ValueType isovalue)
2341 {
2342  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2343  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2344 
2345  std::vector<const typename TreeType::LeafNodeType*> nodes;
2346  tree.getNodes(nodes);
2347 
2348  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2349  typename BoolTreeType::Ptr mask(new BoolTreeType(false));
2350 
2351  level_set_util_internal::MaskIsovalueCrossingVoxels<TreeType> op(tree, nodes, *mask, isovalue);
2352  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
2353 
2355  volume, mask);
2356 }
2357 
2358 
2359 ////////////////////////////////////////
2360 
2361 
2362 template<typename GridOrTreeType>
2363 inline void
2364 extractActiveVoxelSegmentMasks(const GridOrTreeType& volume,
2365  std::vector<typename GridOrTreeType::template ValueConverter<bool>::Type::Ptr>& masks)
2366 {
2367  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2368  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2369  using BoolTreePtrType = typename BoolTreeType::Ptr;
2370  using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2371 
2373  using NodeMaskSegmentPtrType = typename NodeMaskSegmentType::Ptr;
2374  using NodeMaskSegmentPtrVector = typename std::vector<NodeMaskSegmentPtrType>;
2375  using NodeMaskSegmentRawPtrVector = typename std::vector<NodeMaskSegmentType*>;
2376 
2377  /////
2378 
2379  const TreeType& tree = TreeAdapter<GridOrTreeType>::tree(volume);
2380 
2381  BoolTreeType topologyMask(tree, false, TopologyCopy());
2382 
2383  if (topologyMask.hasActiveTiles()) {
2384  topologyMask.voxelizeActiveTiles();
2385  }
2386 
2387  std::vector<BoolLeafNodeType*> leafnodes;
2388  topologyMask.getNodes(leafnodes);
2389 
2390  if (leafnodes.empty()) return;
2391 
2392  // 1. Split node masks into disjoint segments
2393  // Note: The LeafNode origin coord is modified to record the 'leafnodes' array offset.
2394 
2395  std::unique_ptr<NodeMaskSegmentPtrVector[]> nodeSegmentArray(
2396  new NodeMaskSegmentPtrVector[leafnodes.size()]);
2397 
2398  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2400  leafnodes, nodeSegmentArray.get()));
2401 
2402 
2403  // 2. Compute segment connectivity
2404 
2405  tbb::parallel_for(tbb::blocked_range<size_t>(0, leafnodes.size()),
2407  topologyMask, nodeSegmentArray.get()));
2408 
2409  topologyMask.clear();
2410 
2411  size_t nodeSegmentCount = 0;
2412  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2413  nodeSegmentCount += nodeSegmentArray[n].size();
2414  }
2415 
2416  // 3. Group connected segments
2417 
2418  std::deque<NodeMaskSegmentRawPtrVector> nodeSegmentGroups;
2419 
2420  NodeMaskSegmentType* nextSegment = nodeSegmentArray[0][0].get();
2421  while (nextSegment) {
2422 
2423  nodeSegmentGroups.push_back(NodeMaskSegmentRawPtrVector());
2424 
2425  std::vector<NodeMaskSegmentType*>& segmentGroup = nodeSegmentGroups.back();
2426  segmentGroup.reserve(nodeSegmentCount);
2427 
2428  std::deque<NodeMaskSegmentType*> segmentQueue;
2429  segmentQueue.push_back(nextSegment);
2430  nextSegment = nullptr;
2431 
2432  while (!segmentQueue.empty()) {
2433 
2434  NodeMaskSegmentType* segment = segmentQueue.back();
2435  segmentQueue.pop_back();
2436 
2437  if (segment->visited) continue;
2438  segment->visited = true;
2439 
2440  segmentGroup.push_back(segment);
2441 
2442  // queue connected segments
2443  std::vector<NodeMaskSegmentType*>& connections = segment->connections;
2444  for (size_t n = 0, N = connections.size(); n < N; ++n) {
2445  if (!connections[n]->visited) segmentQueue.push_back(connections[n]);
2446  }
2447  }
2448 
2449  // find first unvisited segment
2450  for (size_t n = 0, N = leafnodes.size(); n < N; ++n) {
2451  NodeMaskSegmentPtrVector& nodeSegments = nodeSegmentArray[n];
2452  for (size_t i = 0, I = nodeSegments.size(); i < I; ++i) {
2453  if (!nodeSegments[i]->visited) nextSegment = nodeSegments[i].get();
2454  }
2455  }
2456  }
2457 
2458  // 4. Mask segment groups
2459 
2460  if (nodeSegmentGroups.size() == 1) {
2461 
2462  BoolTreePtrType mask(new BoolTreeType(tree, false, TopologyCopy()));
2463 
2464  if (mask->hasActiveTiles()) {
2465  mask->voxelizeActiveTiles();
2466  }
2467 
2468  masks.push_back(
2470  volume, mask));
2471 
2472  } else if (nodeSegmentGroups.size() > 1) {
2473 
2474  for (size_t n = 0, N = nodeSegmentGroups.size(); n < N; ++n) {
2475 
2476  NodeMaskSegmentRawPtrVector& segmentGroup = nodeSegmentGroups[n];
2477 
2479  tbb::parallel_reduce(tbb::blocked_range<size_t>(0, segmentGroup.size()), op);
2480 
2481  masks.push_back(
2483  volume, op.mask()));
2484  }
2485  }
2486 
2487  // 5. Sort segments in descending order based on the active voxel count.
2488 
2489  if (masks.size() > 1) {
2490  const size_t segmentCount = masks.size();
2491 
2492  std::unique_ptr<size_t[]> segmentOrderArray(new size_t[segmentCount]);
2493  std::unique_ptr<size_t[]> voxelCountArray(new size_t[segmentCount]);
2494 
2495  for (size_t n = 0; n < segmentCount; ++n) {
2496  segmentOrderArray[n] = n;
2497  }
2498 
2499  tbb::parallel_for(tbb::blocked_range<size_t>(0, segmentCount),
2501  masks, voxelCountArray.get()));
2502 
2503  size_t *begin = segmentOrderArray.get();
2504  tbb::parallel_sort(begin, begin + masks.size(), level_set_util_internal::GreaterCount(
2505  voxelCountArray.get()));
2506 
2507  std::vector<BoolTreePtrType> orderedMasks;
2508  orderedMasks.reserve(masks.size());
2509 
2510  for (size_t n = 0; n < segmentCount; ++n) {
2511  orderedMasks.push_back(masks[segmentOrderArray[n]]);
2512  }
2513 
2514  masks.swap(orderedMasks);
2515  }
2516 
2517 } // extractActiveVoxelSegmentMasks()
2518 
2519 
2520 template<typename GridOrTreeType>
2521 inline void
2522 segmentActiveVoxels(const GridOrTreeType& volume,
2523  std::vector<typename GridOrTreeType::Ptr>& segments)
2524 {
2525  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2526  using TreePtrType = typename TreeType::Ptr;
2527  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2528  using BoolTreePtrType = typename BoolTreeType::Ptr;
2529 
2530  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2531 
2532  // 1. Segment active topology mask
2533  std::vector<BoolTreePtrType> maskSegmentArray;
2534  extractActiveVoxelSegmentMasks(inputTree, maskSegmentArray);
2535 
2536  const size_t numSegments = maskSegmentArray.size();
2537 
2538  if (numSegments < 2) {
2539  // single segment early-out
2540  TreePtrType segment(new TreeType(inputTree));
2541  segments.push_back(
2543  volume, segment));
2544  return;
2545  }
2546 
2547  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2548 
2549  // 2. Export segments
2550  std::vector<TreePtrType> outputSegmentArray(numSegments);
2551 
2552  tbb::parallel_for(segmentRange,
2553  level_set_util_internal::MaskedCopy<TreeType>(inputTree, outputSegmentArray,
2554  maskSegmentArray));
2555 
2556  for (size_t n = 0, N = numSegments; n < N; ++n) {
2557  segments.push_back(
2559  volume, outputSegmentArray[n]));
2560  }
2561 }
2562 
2563 
2564 template<typename GridOrTreeType>
2565 inline void
2566 segmentSDF(const GridOrTreeType& volume, std::vector<typename GridOrTreeType::Ptr>& segments)
2567 {
2568  using TreeType = typename TreeAdapter<GridOrTreeType>::TreeType;
2569  using TreePtrType = typename TreeType::Ptr;
2570  using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2571  using BoolTreePtrType = typename BoolTreeType::Ptr;
2572 
2573  const TreeType& inputTree = TreeAdapter<GridOrTreeType>::tree(volume);
2574 
2575  // 1. Mask zero crossing voxels
2576  BoolTreePtrType mask = extractIsosurfaceMask(inputTree, lsutilGridZero<GridOrTreeType>());
2577 
2578  // 2. Segment the zero crossing mask
2579  std::vector<BoolTreePtrType> maskSegmentArray;
2580  extractActiveVoxelSegmentMasks(*mask, maskSegmentArray);
2581 
2582  const size_t numSegments = maskSegmentArray.size();
2583 
2584  if (numSegments < 2) {
2585  // single segment early-out
2586  TreePtrType segment(new TreeType(inputTree));
2587  segments.push_back(
2589  volume, segment));
2590  return;
2591  }
2592 
2593  const tbb::blocked_range<size_t> segmentRange(0, numSegments);
2594 
2595 
2596  // 3. Expand zero crossing mask to capture sdf narrow band
2597  tbb::parallel_for(segmentRange,
2598  level_set_util_internal::ExpandNarrowbandMask<TreeType>(inputTree, maskSegmentArray));
2599 
2600  // 4. Export sdf segments
2601  std::vector<TreePtrType> outputSegmentArray(numSegments);
2602 
2603  tbb::parallel_for(segmentRange, level_set_util_internal::MaskedCopy<TreeType>(
2604  inputTree, outputSegmentArray, maskSegmentArray));
2605 
2606  tbb::parallel_for(segmentRange,
2607  level_set_util_internal::FloodFillSign<TreeType>(inputTree, outputSegmentArray));
2608 
2609 
2610  for (size_t n = 0, N = numSegments; n < N; ++n) {
2611  segments.push_back(
2613  volume, outputSegmentArray[n]));
2614  }
2615 }
2616 
2617 } // namespace tools
2618 } // namespace OPENVDB_VERSION_NAME
2619 } // namespace openvdb
2620 
2621 #endif // OPENVDB_TOOLS_LEVEL_SET_UTIL_HAS_BEEN_INCLUDED
2622 
2623 // Copyright (c) 2012-2018 DreamWorks Animation LLC
2624 // All rights reserved. This software is distributed under the
2625 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
typename InputTreeType::template ValueConverter< bool >::Type BoolTreeType
GA_API const UT_StringHolder dist
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
typename TreeType::template ValueConverter< bool >::Type::Ptr BoolTreePtrType
GLenum GLint * range
Definition: glcorearb.h:1924
MaskInteriorTiles(ValueType isovalue, const TreeType &tree, InternalNodeType **maskNodes)
Definition: LevelSetUtil.h:264
void sdfToFogVolume(GridType &grid, typename GridType::ValueType cutoffDistance=lsutilGridMax< GridType >())
Threaded method to convert a sparse level set/SDF into a sparse fog volume.
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:512
void setValueOn(const Coord &xyz)
Mark the voxel at the given coordinates as active but don't change its value.
Definition: LeafNode.h:440
MaskedCopy(const TreeType &tree, std::vector< TreeTypePtr > &segments, std::vector< BoolTreeTypePtr > &masks)
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:551
png_infop png_color_16p * background
Definition: png.h:2326
void setOff(Index32 n)
Set the nth bit off.
Definition: NodeMasks.h:459
GLdouble GLdouble GLdouble z
Definition: glcorearb.h:847
FillMaskBoundary(const TreeType &tree, ValueType isovalue, const BoolTreeType &fillMask, const BoolLeafNodeType **fillNodes, BoolLeafNodeType **newNodes)
Definition: LevelSetUtil.h:541
GLint GLuint mask
Definition: glcorearb.h:123
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:189
FloodFillSign(const TreeType &tree, std::vector< TreeTypePtr > &segments)
bool operator()(const size_t &lhs, const size_t &rhs) const
ExpandLeafNodeRegion(const TreeType &distTree, BoolTreeType &maskTree, std::vector< BoolLeafNodeType * > &maskNodes)
GLint y
Definition: glcorearb.h:102
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
Definition: LevelSetUtil.h:538
LabelBoundaryVoxels(ValueType isovalue, const LeafNodeType **nodes, CharLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:343
png_uint_32 i
Definition: png.h:2877
ComputeActiveVoxelCount(std::vector< VolumePtrType > &segments, size_t *countArray)
ConnectNodeMaskSegments(const TreeType &tree, NodeMaskSegmentVector *nodeMaskArray)
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractEnclosedRegion(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >(), const typename TreeAdapter< GridOrTreeType >::TreeType::template ValueConverter< bool >::Type *fillMask=nullptr)
Extracts the interior regions of a signed distance field and topologically enclosed (watertight) regi...
void operator()(const tbb::blocked_range< size_t > &range)
Definition: LevelSetUtil.h:309
This adapter allows code that is templated on a Tree type to accept either a Tree type or a Grid type...
Definition: Grid.h:944
void signedFloodFillWithValues(TreeOrLeafManagerT &tree, const typename TreeOrLeafManagerT::ValueType &outsideWidth, const typename TreeOrLeafManagerT::ValueType &insideWidth, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
typename std::vector< NodeMaskSegmentTypePtr > NodeMaskSegmentVector
void addLeaf(LeafNodeT *leaf)
Add the specified leaf to this tree, possibly creating a child branch in the process. If the leaf node already exists, replace it.
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
std::shared_ptr< T > SharedPtr
Definition: Types.h:139
static BoolTreePtrType constructMask(const TreeType &, BoolTreePtrType &maskTree)
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
TreeType::template ValueConverter< char >::Type::Ptr computeEnclosedRegionMask(const TreeType &tree, typename TreeType::ValueType isovalue, const typename TreeType::template ValueConverter< bool >::Type *fillMask)
Constructs a memory light char tree that represents the exterior region with +1 and the interior regi...
Definition: LevelSetUtil.h:978
SYS_FORCE_INLINE const_iterator end() const
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:481
GLdouble n
Definition: glcorearb.h:2007
void setOrigin(const Coord &origin)
Set the grid index coordinates of this node's local origin.
Definition: LeafNode.h:198
Negative active values are set 0, everything else is set to 1.
Definition: LevelSetUtil.h:338
const LeafNodeT * probeConstLeaf(const Coord &xyz) const
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
SegmentNodeMask(std::vector< NodeType * > &nodes, NodeMaskSegmentVector *nodeMaskArray)
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
Definition: ImathFun.h:55
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Index64 onVoxelCount() const
Return the number of voxels marked On.
Definition: LeafNode.h:166
ExpandNarrowbandMask(const TreeType &tree, std::vector< BoolTreeTypePtr > &segments)
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:133
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim...
Definition: LeafNode.h:64
void extractActiveVoxelSegmentMasks(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::template ValueConverter< bool >::Type::Ptr > &masks)
Return a mask for each connected component of the given grid's active voxels.
MaskSegmentGroup(const std::vector< NodeMaskSegmentType * > &segments)
GLuint GLuint end
Definition: glcorearb.h:474
void segmentActiveVoxels(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint active topology components into distinct grids or trees.
const NodeMaskType & getValueMask() const
Definition: LeafNode.h:889
LeafNodeT * touchLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists, create one, but preserve the values and active states of all voxels.
void nodeMaskSegmentation(const NodeType &node, std::vector< typename NodeMaskSegment< NodeType >::Ptr > &segments)
MaskIsovalueCrossingVoxels(const InputTreeType &inputTree, const std::vector< const InputLeafNodeType * > &inputLeafNodes, BoolTreeType &maskTree, InputValueType iso)
static BoolGridPtrType constructMask(const GridType &grid, BoolTreePtrType &maskTree)
TreeType::template ValueConverter< bool >::Type::Ptr computeInteriorMask(const TreeType &tree, typename TreeType::ValueType iso)
void operator()(const tbb::blocked_range< size_t > &range) const
static TreeType & tree(TreeType &t)
Definition: Grid.h:961
GLenum GLsizei GLsizei GLint * values
Definition: glcorearb.h:1601
GLboolean * data
Definition: glcorearb.h:130
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z), or nullptr if no such node exists...
const Coord & origin() const
Return the grid index coordinates of this node's local origin.
Definition: LeafNode.h:201
math::Transform & transform()
Return a reference to this grid's transform, which might be shared with other grids.
Definition: Grid.h:347
Container class that associates a tree with a transform and metadata.
Definition: Grid.h:55
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
void operator()(const tbb::blocked_range< size_t > &range) const
typedef int
Definition: png.h:1175
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:224
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
void operator()(const tbb::blocked_range< size_t > &range) const
FillLeafNodeVoxels(const TreeType &tree, std::vector< BoolLeafNodeType * > &maskNodes)
GLint GLenum GLint x
Definition: glcorearb.h:408
GLuint GLfloat * val
Definition: glcorearb.h:1607
void setValueOn(const Coord &xyz, const ValueType &value)
Set the value of the voxel at the given coordinates and mark the voxel as active. ...
GA_API const UT_StringHolder N
void segmentSDF(const GridOrTreeType &volume, std::vector< typename GridOrTreeType::Ptr > &segments)
Separates disjoint SDF surfaces into distinct grids or trees.
GridOrTreeType::template ValueConverter< bool >::Type::Ptr sdfInteriorMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue=lsutilGridZero< GridOrTreeType >())
Threaded method to construct a boolean mask that represents interior regions in a signed distance fie...
GLint GLint GLsizei GLint GLenum GLenum type
Definition: glcorearb.h:107
#define SIZE
Definition: simple.C:40
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:518
static TreeTypePtr construct(const TreeType &, TreeTypePtr &tree)
GridOrTreeType::template ValueConverter< bool >::Type::Ptr extractIsosurfaceMask(const GridOrTreeType &volume, typename GridOrTreeType::ValueType isovalue)
Return a mask of the voxels that intersect the implicit surface with the given isovalue.
void operator()(const tbb::blocked_range< size_t > &range) const
typename TreeType::template ValueConverter< bool >::Type BoolTreeType
typename hboost::mpl::at< NodeChainType, hboost::mpl::int_< 1 > >::type InternalNodeType
const std::enable_if<!VecTraits< T >::IsVec, T >::type & min(const T &a, const T &b)
Definition: Composite.h:129
void operator()(const tbb::blocked_range< size_t > &range) const
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:390
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:135
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:267
MaskInteriorVoxels(ValueType isovalue, const LeafNodeType **nodes, BoolLeafNodeType **maskNodes)
Definition: LevelSetUtil.h:218
PopulateTree(TreeType &tree, LeafNodeType **leafnodes, const size_t *nodexIndexMap, ValueType background)
Definition: LevelSetUtil.h:292
TreeType & tree() const
Return a reference to the tree associated with this accessor.
void operator()(const tbb::blocked_range< size_t > &range) const
Definition: LevelSetUtil.h:349