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