OpenVDB  1.2.0
VolumeToMesh.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 
31 #ifndef OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
32 #define OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/tree/ValueAccessor.h>
35 #include <openvdb/util/Util.h> // for COORD_OFFSETS
36 #include <openvdb/math/Operators.h> // for ISGradient
37 #include <openvdb/tools/Morphology.h> // for dilateVoxels()
38 #include <openvdb/tree/LeafManager.h>
39 
40 #include <boost/scoped_array.hpp>
41 #include <boost/scoped_ptr.hpp>
42 
43 #include <tbb/blocked_range.h>
44 #include <tbb/parallel_for.h>
45 #include <tbb/parallel_reduce.h>
46 
47 #include <vector>
48 
49 namespace openvdb {
51 namespace OPENVDB_VERSION_NAME {
52 namespace tools {
53 
55 
56 
57 // Wrapper functions for the VolumeToMesh converter
58 
59 
66 template<typename GridType>
67 void
69  const GridType& grid,
70  std::vector<Vec3s>& points,
71  std::vector<Vec4I>& quads,
72  double isovalue = 0.0);
73 
74 
83 template<typename GridType>
84 void
86  const GridType& grid,
87  std::vector<Vec3s>& points,
88  std::vector<Vec3I>& triangles,
89  std::vector<Vec4I>& quads,
90  double isovalue = 0.0,
91  double adaptivity = 0.0);
92 
93 
95 
96 
98 enum {
100 };
101 
102 
105 {
106 public:
108  : mNumQuads(0)
109  , mNumTriangles(0)
110  , mQuads(NULL)
111  , mTriangles(NULL)
112  , mQuadFlags(NULL)
113  , mTriangleFlags(NULL)
114  {
115  }
116 
117  void resetQuads(size_t size)
118  {
119  mNumQuads = size;
120  mQuads.reset(new openvdb::Vec4I[mNumQuads]);
121  mQuadFlags.reset(new char[mNumQuads]);
122  }
123 
124  void clearQuads()
125  {
126  mNumQuads = 0;
127  mQuads.reset(NULL);
128  mQuadFlags.reset(NULL);
129  }
130 
131  void resetTriangles(size_t size)
132  {
133  mNumTriangles = size;
134  mTriangles.reset(new openvdb::Vec3I[mNumTriangles]);
135  mTriangleFlags.reset(new char[mNumTriangles]);
136  }
137 
139  {
140  mNumTriangles = 0;
141  mTriangles.reset(NULL);
142  mTriangleFlags.reset(NULL);
143  }
144 
145  const size_t& numQuads() const { return mNumQuads; }
146  const size_t& numTriangles() const { return mNumTriangles; }
147 
148  // polygon accessor methods
149  openvdb::Vec4I& quad(size_t n) { return mQuads[n]; }
150  const openvdb::Vec4I& quad(size_t n) const { return mQuads[n]; }
151 
152  openvdb::Vec3I& triangle(size_t n) { return mTriangles[n]; }
153  const openvdb::Vec3I& triangle(size_t n) const { return mTriangles[n]; }
154 
155  // polygon flags accessor methods
156  char& quadFlags(size_t n) { return mQuadFlags[n]; }
157  const char& quadFlags(size_t n) const { return mQuadFlags[n]; }
158 
159  char& triangleFlags(size_t n) { return mTriangleFlags[n]; }
160  const char& triangleFlags(size_t n) const { return mTriangleFlags[n]; }
161 
162 private:
163  size_t mNumQuads, mNumTriangles;
164  boost::scoped_array<openvdb::Vec4I> mQuads;
165  boost::scoped_array<openvdb::Vec3I> mTriangles;
166  boost::scoped_array<char> mQuadFlags, mTriangleFlags;
167 };
168 
169 
172 typedef boost::scoped_array<openvdb::Vec3s> PointList;
173 typedef boost::scoped_array<PolygonPool> PolygonPoolList;
175 
176 
178 
179 
182 {
183 public:
186  VolumeToMesh(double isovalue = 0, double adaptivity = 0);
187 
188  PointList& pointList();
189  const size_t& pointListSize() const;
190 
191  PolygonPoolList& polygonPoolList();
192  const PolygonPoolList& polygonPoolList() const;
193  const size_t& polygonPoolListSize() const;
194 
197  template<typename GridT>
198  void operator()(const GridT&);
199 
200 
224  void setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity = 0,
225  bool smoothSeams = true);
226 
227 
231  void setSurfaceMask(const GridBase::ConstPtr& mask, bool invertMask = false);
232 
233 
236  void setSpatialAdaptivity(const GridBase::ConstPtr& grid);
237 
238 
241  void setAdaptivityMask(const TreeBase::ConstPtr& tree);
242 
243 
247  void partition(unsigned partitions = 1, unsigned activePart = 0);
248 
249 
250 private:
251 
252  PointList mPoints;
253  PolygonPoolList mPolygons;
254 
255  size_t mPointListSize, mPolygonPoolListSize;
256  double mIsovalue, mPrimAdaptivity, mSecAdaptivity;
257 
258  GridBase::ConstPtr mRefGrid, mSurfaceMaskGrid, mAdaptivityGrid;
259  TreeBase::ConstPtr mAdaptivityMaskTree;
260 
261  TreeBase::Ptr mRefEdgeTree, mRefTopologyMaskTree, mSeamPointTree;
262 
263  bool mSmoothSeams, mInvertSurfaceMask;
264  unsigned mPartitions, mActivePart;
265 };
266 
267 
269 
270 
278  const std::vector<Vec3d>& points,
279  const std::vector<Vec3d>& normals)
280 {
281  typedef math::Mat3d Mat3d;
282 
283  Vec3d avgPos(0.0);
284 
285  if (points.empty()) return avgPos;
286 
287  for (size_t n = 0, N = points.size(); n < N; ++n) {
288  avgPos += points[n];
289  }
290 
291  avgPos /= double(points.size());
292 
293  // Unique components of the 3x3 A^TA matrix, where A is
294  // the matrix of normals.
295  double m00=0,m01=0,m02=0,
296  m11=0,m12=0,
297  m22=0;
298 
299  // The rhs vector, A^Tb, where b = n dot p
300  Vec3d rhs(0.0);
301 
302  for (size_t n = 0, N = points.size(); n < N; ++n) {
303 
304  const Vec3d& n_ref = normals[n];
305 
306  // A^TA
307  m00 += n_ref[0] * n_ref[0]; // diagonal
308  m11 += n_ref[1] * n_ref[1];
309  m22 += n_ref[2] * n_ref[2];
310 
311  m01 += n_ref[0] * n_ref[1]; // Upper-tri
312  m02 += n_ref[0] * n_ref[2];
313  m12 += n_ref[1] * n_ref[2];
314 
315  // A^Tb (centered around the origin)
316  rhs += n_ref * n_ref.dot(points[n] - avgPos);
317  }
318 
319  Mat3d A(m00,m01,m02,
320  m01,m11,m12,
321  m02,m12,m22);
322 
323  /*
324  // Inverse
325  const double det = A.det();
326  if (det > 0.01) {
327  Mat3d A_inv = A.adjoint();
328  A_inv *= (1.0 / det);
329 
330  return avgPos + A_inv * rhs;
331  }
332  */
333 
334  // Compute the pseudo inverse
335 
336  math::Mat3d eigenVectors;
337  Vec3d eigenValues;
338 
339  diagonalizeSymmetricMatrix(A, eigenVectors, eigenValues, 300);
340 
341  Mat3d D = Mat3d::identity();
342 
343 
344  double tolerance = std::max(std::abs(eigenValues[0]), std::abs(eigenValues[1]));
345  tolerance = std::max(tolerance, std::abs(eigenValues[2]));
346  tolerance *= 0.01;
347 
348  int clamped = 0;
349  for (int i = 0; i < 3; ++i ) {
350  if (std::abs(eigenValues[i]) < tolerance) {
351  D[i][i] = 0.0;
352  ++clamped;
353  } else {
354  D[i][i] = 1.0 / eigenValues[i];
355  }
356  }
357 
358  // Assemble the pseudo inverse and calc. the intersection point
359  if (clamped < 3) {
360  Mat3d pseudoInv = eigenVectors * D * eigenVectors.transpose();
361  return avgPos + pseudoInv * rhs;
362  }
363 
364  return avgPos;
365 }
366 
367 
369 
370 
371 // Internal utility methods
372 
373 
374 namespace internal {
375 
376 
377 // Bit-flags
378 enum { INSIDE = 0x1, XEDGE = 0x2, YEDGE = 0x4, ZEDGE = 0x8 };
379 
380 const bool sAmbiguous[256] =
381  {0,0,0,0,0,1,0,0,
382  0,0,1,0,0,0,0,0,
383  0,0,1,0,1,1,1,0,
384  1,0,1,0,1,0,1,0,
385  0,1,0,0,1,1,0,0,
386  1,1,1,0,1,1,0,0,
387  0,0,0,0,1,1,0,0,
388  1,0,1,0,1,1,1,0,
389  0,1,1,1,0,1,0,0,
390  1,1,1,1,0,0,0,0,
391  1,1,1,1,1,1,1,1,
392  1,1,1,1,1,1,1,1,
393  0,1,0,0,0,1,0,0,
394  1,1,1,1,0,1,0,0,
395  0,0,0,0,0,1,0,0,
396  1,1,1,1,1,1,1,0,
397  0,1,1,1,1,1,1,1,
398  0,0,1,0,0,0,0,0,
399  0,0,1,0,1,1,1,1,
400  0,0,1,0,0,0,1,0,
401  1,1,1,1,1,1,1,1,
402  1,1,1,1,1,1,1,1,
403  0,0,0,0,1,1,1,1,
404  0,0,1,0,1,1,1,0,
405  0,1,1,1,0,1,0,1,
406  0,0,1,1,0,0,0,0,
407  0,0,1,1,0,1,1,1,
408  0,0,1,1,0,0,1,0,
409  0,1,0,1,0,1,0,1,
410  0,1,1,1,0,1,0,0,
411  0,0,0,0,0,1,0,0,
412  0,0,1,0,0,0,0,0};
413 
414 
415 template<class AccessorT>
416 inline bool isAmbiguous(const AccessorT& accessor, const Coord& ijk,
417  typename AccessorT::ValueType isovalue, int dim)
418 {
419  unsigned signs = 0;
420  Coord coord = ijk; // i, j, k
421  if (accessor.getValue(coord) < isovalue) signs |= 1u;
422  coord[0] += dim; // i+dim, j, k
423  if (accessor.getValue(coord) < isovalue) signs |= 2u;
424  coord[2] += dim; // i+dim, j, k+dim
425  if (accessor.getValue(coord) < isovalue) signs |= 4u;
426  coord[0] = ijk[0]; // i, j, k+dim
427  if (accessor.getValue(coord) < isovalue) signs |= 8u;
428  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
429  if (accessor.getValue(coord) < isovalue) signs |= 16u;
430  coord[0] += dim; // i+dim, j+dim, k
431  if (accessor.getValue(coord) < isovalue) signs |= 32u;
432  coord[2] += dim; // i+dim, j+dim, k+dim
433  if (accessor.getValue(coord) < isovalue) signs |= 64u;
434  coord[0] = ijk[0]; // i, j+dim, k+dim
435  if (accessor.getValue(coord) < isovalue) signs |= 128u;
436  return sAmbiguous[signs];
437 }
438 
439 
440 template<class AccessorT>
441 inline bool
442 isNonManifold(const AccessorT& accessor, const Coord& ijk,
443  typename AccessorT::ValueType isovalue, const int dim)
444 {
445  int hDim = dim >> 1;
446  bool m, p[8]; // Corner signs
447 
448  Coord coord = ijk; // i, j, k
449  p[0] = accessor.getValue(coord) < isovalue;
450  coord[0] += dim; // i+dim, j, k
451  p[1] = accessor.getValue(coord) < isovalue;
452  coord[2] += dim; // i+dim, j, k+dim
453  p[2] = accessor.getValue(coord) < isovalue;
454  coord[0] = ijk[0]; // i, j, k+dim
455  p[3] = accessor.getValue(coord) < isovalue;
456  coord[1] += dim; coord[2] = ijk[2]; // i, j+dim, k
457  p[4] = accessor.getValue(coord) < isovalue;
458  coord[0] += dim; // i+dim, j+dim, k
459  p[5] = accessor.getValue(coord) < isovalue;
460  coord[2] += dim; // i+dim, j+dim, k+dim
461  p[6] = accessor.getValue(coord) < isovalue;
462  coord[0] = ijk[0]; // i, j+dim, k+dim
463  p[7] = accessor.getValue(coord) < isovalue;
464 
465  // Check if the corner sign configuration is ambiguous
466  unsigned signs = 0;
467  if (p[0]) signs |= 1u;
468  if (p[1]) signs |= 2u;
469  if (p[2]) signs |= 4u;
470  if (p[3]) signs |= 8u;
471  if (p[4]) signs |= 16u;
472  if (p[5]) signs |= 32u;
473  if (p[6]) signs |= 64u;
474  if (p[7]) signs |= 128u;
475  if (sAmbiguous[signs]) return true;
476 
477  // Manifold check
478 
479  // Evaluate edges
480  int i = ijk[0], ip = ijk[0] + hDim, ipp = ijk[0] + dim;
481  int j = ijk[1], jp = ijk[1] + hDim, jpp = ijk[1] + dim;
482  int k = ijk[2], kp = ijk[2] + hDim, kpp = ijk[2] + dim;
483 
484  // edge 1
485  coord.reset(ip, j, k);
486  m = accessor.getValue(coord) < isovalue;
487  if (p[0] != m && p[1] != m) return true;
488 
489  // edge 2
490  coord.reset(ipp, j, kp);
491  m = accessor.getValue(coord) < isovalue;
492  if (p[1] != m && p[2] != m) return true;
493 
494  // edge 3
495  coord.reset(ip, j, kpp);
496  m = accessor.getValue(coord) < isovalue;
497  if (p[2] != m && p[3] != m) return true;
498 
499  // edge 4
500  coord.reset(i, j, kp);
501  m = accessor.getValue(coord) < isovalue;
502  if (p[0] != m && p[3] != m) return true;
503 
504  // edge 5
505  coord.reset(ip, jpp, k);
506  m = accessor.getValue(coord) < isovalue;
507  if (p[4] != m && p[5] != m) return true;
508 
509  // edge 6
510  coord.reset(ipp, jpp, kp);
511  m = accessor.getValue(coord) < isovalue;
512  if (p[5] != m && p[6] != m) return true;
513 
514  // edge 7
515  coord.reset(ip, jpp, kpp);
516  m = accessor.getValue(coord) < isovalue;
517  if (p[6] != m && p[7] != m) return true;
518 
519  // edge 8
520  coord.reset(i, jpp, kp);
521  m = accessor.getValue(coord) < isovalue;
522  if (p[7] != m && p[4] != m) return true;
523 
524  // edge 9
525  coord.reset(i, jp, k);
526  m = accessor.getValue(coord) < isovalue;
527  if (p[0] != m && p[4] != m) return true;
528 
529  // edge 10
530  coord.reset(ipp, jp, k);
531  m = accessor.getValue(coord) < isovalue;
532  if (p[1] != m && p[5] != m) return true;
533 
534  // edge 11
535  coord.reset(ipp, jp, kpp);
536  m = accessor.getValue(coord) < isovalue;
537  if (p[2] != m && p[6] != m) return true;
538 
539 
540  // edge 12
541  coord.reset(i, jp, kpp);
542  m = accessor.getValue(coord) < isovalue;
543  if (p[3] != m && p[7] != m) return true;
544 
545 
546  // Evaluate faces
547 
548  // face 1
549  coord.reset(ip, jp, k);
550  m = accessor.getValue(coord) < isovalue;
551  if (p[0] != m && p[1] != m && p[4] != m && p[5] != m) return true;
552 
553  // face 2
554  coord.reset(ipp, jp, kp);
555  m = accessor.getValue(coord) < isovalue;
556  if (p[1] != m && p[2] != m && p[5] != m && p[6] != m) return true;
557 
558  // face 3
559  coord.reset(ip, jp, kpp);
560  m = accessor.getValue(coord) < isovalue;
561  if (p[2] != m && p[3] != m && p[6] != m && p[7] != m) return true;
562 
563  // face 4
564  coord.reset(i, jp, kp);
565  m = accessor.getValue(coord) < isovalue;
566  if (p[0] != m && p[3] != m && p[4] != m && p[7] != m) return true;
567 
568  // face 5
569  coord.reset(ip, j, kp);
570  m = accessor.getValue(coord) < isovalue;
571  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m) return true;
572 
573  // face 6
574  coord.reset(ip, jpp, kp);
575  m = accessor.getValue(coord) < isovalue;
576  if (p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
577 
578  // test cube center
579  coord.reset(ip, jp, kp);
580  m = accessor.getValue(coord) < isovalue;
581  if (p[0] != m && p[1] != m && p[2] != m && p[3] != m &&
582  p[4] != m && p[5] != m && p[6] != m && p[7] != m) return true;
583 
584  return false;
585 }
586 
587 
589 
590 
591 template <class LeafType>
592 inline void
593 mergeVoxels(LeafType& leaf, const Coord& start, int dim, int regionId)
594 {
595  Coord ijk, end = start;
596  end[0] += dim;
597  end[1] += dim;
598  end[2] += dim;
599 
600  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
601  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
602  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
603  if(leaf.isValueOn(ijk)) leaf.setValue(ijk, regionId);
604  }
605  }
606  }
607 }
608 
609 
610 // Note that we must use ValueType::value_type or else Visual C++ gets confused
611 // thinking that it is a constructor.
612 template <class LeafType>
613 inline bool
614 isMergable(LeafType& leaf, const Coord& start, int dim,
615  typename LeafType::ValueType::value_type adaptivity)
616 {
617  if (adaptivity < 1e-6) return false;
618 
619  typedef typename LeafType::ValueType VecT;
620  Coord ijk, end = start;
621  end[0] += dim;
622  end[1] += dim;
623  end[2] += dim;
624 
625  std::vector<VecT> norms;
626  for (ijk[0] = start[0]; ijk[0] < end[0]; ++ijk[0]) {
627  for (ijk[1] = start[1]; ijk[1] < end[1]; ++ijk[1]) {
628  for (ijk[2] = start[2]; ijk[2] < end[2]; ++ijk[2]) {
629 
630  if(!leaf.isValueOn(ijk)) continue;
631  norms.push_back(leaf.getValue(ijk));
632  }
633  }
634  }
635 
636  size_t N = norms.size();
637  for (size_t ni = 0; ni < N; ++ni) {
638  VecT n_i = norms[ni];
639  for (size_t nj = 0; nj < N; ++nj) {
640  VecT n_j = norms[nj];
641  if ((1.0 - n_i.dot(n_j)) > adaptivity) return false;
642  }
643  }
644  return true;
645 }
646 
647 
649 
650 
651 template <class TreeT>
653 {
654 public:
655  typedef std::vector<const typename TreeT::LeafNodeType *> ListT;
656 
657  LeafCPtrList(const TreeT& tree)
658  {
659  mLeafNodes.reserve(tree.leafCount());
660  typename TreeT::LeafCIter iter = tree.cbeginLeaf();
661  for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
662  }
663 
664  size_t size() const { return mLeafNodes.size(); }
665 
666  const typename TreeT::LeafNodeType* operator[](size_t n) const
667  { return mLeafNodes[n]; }
668 
669  tbb::blocked_range<size_t> getRange() const
670  { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
671 
672  const ListT& getList() const { return mLeafNodes; }
673 
674 private:
675  ListT mLeafNodes;
676 };
677 
678 
679 template <class TreeT>
681 {
682 public:
683  typedef std::vector<typename TreeT::LeafNodeType *> ListT;
684 
685  LeafPtrList(TreeT& tree)
686  {
687  mLeafNodes.reserve(tree.leafCount());
688  typename TreeT::LeafIter iter = tree.beginLeaf();
689  for ( ; iter; ++iter) mLeafNodes.push_back(iter.getLeaf());
690  }
691 
692  size_t size() const { return mLeafNodes.size(); }
693 
694  typename TreeT::LeafNodeType* operator[](size_t n) const
695  { return mLeafNodes[n]; }
696 
697  tbb::blocked_range<size_t> getRange() const
698  { return tbb::blocked_range<size_t>(0, mLeafNodes.size()); }
699 
700  const ListT& getList() const { return mLeafNodes; }
701 
702 private:
703  ListT mLeafNodes;
704 };
705 
706 
708 
709 template<typename DistTreeT>
711 {
712  typedef typename DistTreeT::ValueType DistValueT;
713  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
714  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
715  typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
716 
718  : mDistTree(NULL)
719  , mEdgeTree(NULL)
720  , mTopologyMaskTree(NULL)
721  , mSeamPointTree(NULL)
722  , mSeamMaskTree(typename BoolTreeT::Ptr())
723  , mSmoothingMaskTree(typename BoolTreeT::Ptr())
724  , mInternalAdaptivity(DistValueT(0.0))
725  {
726  }
727 
728  bool isValid() const
729  {
730  return mDistTree && mEdgeTree && mTopologyMaskTree && mSeamMaskTree;
731  }
732 
733  const DistTreeT* mDistTree;
737  typename BoolTreeT::Ptr mSeamMaskTree, mSmoothingMaskTree;
739 };
740 
741 
743 
744 
745 template <class DistTreeT>
746 class Count
747 {
748 public:
749  Count(const LeafPtrList<DistTreeT>&, std::vector<size_t>&);
750  inline Count(const Count<DistTreeT>&);
751 
752  void runParallel();
753  void runSerial();
754 
755  inline void operator()(const tbb::blocked_range<size_t>&) const;
756 private:
757  const LeafPtrList<DistTreeT>& mLeafNodes;
758  std::vector<size_t>& mLeafRegionCount;
759 };
760 
761 
762 template <class DistTreeT>
764  const LeafPtrList<DistTreeT>& leafs,
765  std::vector<size_t>& leafRegionCount)
766  : mLeafNodes(leafs)
767  , mLeafRegionCount(leafRegionCount)
768 {
769 }
770 
771 
772 template <class DistTreeT>
773 inline
775  : mLeafNodes(rhs.mLeafNodes)
776  , mLeafRegionCount(rhs.mLeafRegionCount)
777 {
778 }
779 
780 
781 template <class DistTreeT>
782 void
784 {
785  tbb::parallel_for(mLeafNodes.getRange(), *this);
786 }
787 
788 
789 template <class DistTreeT>
790 void
792 {
793  (*this)(mLeafNodes.getRange());
794 }
795 
796 
797 template <class DistTreeT>
798 inline void
799 Count<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
800 {
801  for (size_t n = range.begin(); n != range.end(); ++n) {
802  mLeafRegionCount[n] = size_t(mLeafNodes[n]->onVoxelCount());
803  }
804 }
805 
806 
808 
809 
810 template <class DistTreeT>
811 class Merge
812 {
813 public:
814  typedef typename DistTreeT::ValueType DistValueT;
815  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
816  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
817  typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
819 
820  Merge(
821  const DistTreeT& distTree,
822  LeafPtrList<IntTreeT>& auxLeafs,
823  std::vector<size_t>& leafRegionCount,
824  const DistValueT iso,
825  const DistValueT adaptivity,
826  const BoolTreeT * pointMask = NULL);
827 
828  inline Merge(const Merge<DistTreeT>&);
829 
831 
833  const math::Transform& distGridXForm, const FloatGridT& adaptivityField);
834 
835  void run(bool threaded = true);
836 
837  void operator()(const tbb::blocked_range<size_t>&) const;
838 
839 private:
840  const DistTreeT& mDistTree;
841  LeafPtrList<IntTreeT>& mAuxLeafs;
842  std::vector<size_t>& mLeafRegionCount;
843  const DistValueT mIsovalue, mAdaptivity;
844  const ReferenceData<DistTreeT>* mRefData;
845 
846  const math::Transform * mTransform;
847  const FloatGridT * mAdaptivityGrid;
848  const BoolTreeT * mPointMask;
849 };
850 
851 
852 template <class DistTreeT>
854  const DistTreeT& distTree,
855  LeafPtrList<IntTreeT>& auxLeafs,
856  std::vector<size_t>& leafRegionCount,
857  const DistValueT iso,
858  const DistValueT adaptivity,
859  const BoolTreeT * pointMask)
860  : mDistTree(distTree)
861  , mAuxLeafs(auxLeafs)
862  , mLeafRegionCount(leafRegionCount)
863  , mIsovalue(iso)
864  , mAdaptivity(adaptivity)
865  , mRefData(NULL)
866  , mTransform(NULL)
867  , mAdaptivityGrid(NULL)
868  , mPointMask(pointMask)
869 {
870 }
871 
872 
873 template <class DistTreeT>
874 inline
876  : mDistTree(rhs.mDistTree)
877  , mAuxLeafs(rhs.mAuxLeafs)
878  , mLeafRegionCount(rhs.mLeafRegionCount)
879  , mIsovalue(rhs.mIsovalue)
880  , mAdaptivity(rhs.mAdaptivity)
881  , mRefData(rhs.mRefData)
882  , mTransform(rhs.mTransform)
883  , mAdaptivityGrid(rhs.mAdaptivityGrid)
884  , mPointMask(rhs.mPointMask)
885 {
886 }
887 
888 
889 template <class DistTreeT>
890 void
891 Merge<DistTreeT>::run(bool threaded)
892 {
893  if (threaded) {
894  tbb::parallel_for(mAuxLeafs.getRange(), *this);
895  } else {
896  (*this)(mAuxLeafs.getRange());
897  }
898 }
899 
900 
901 template <class DistTreeT>
902 void
904 {
905  mRefData = &refData;
906 }
907 
908 
909 template <class DistTreeT>
910 void
912  const math::Transform& distGridXForm, const FloatGridT& adaptivityField)
913 {
914  mTransform = &distGridXForm;
915  mAdaptivityGrid = &adaptivityField;
916 }
917 
918 
919 template <class DistTreeT>
920 void
921 Merge<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
922 {
923  typedef math::Vec3<DistValueT> Vec3T;
924  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
925  typedef typename IntTreeT::LeafNodeType IntLeafT;
926  typedef typename BoolLeafT::template ValueConverter<Vec3T>::Type Vec3LeafT;
927  typedef typename BoolLeafT::template ValueConverter<DistValueT>::Type DistLeafT;
928 
929  typedef typename IntLeafT::ValueOnIter IntIterT;
930  typedef typename BoolLeafT::ValueOnCIter BoolCIterT;
931 
932  typedef typename tree::ValueAccessor<BoolTreeT> BoolTreeAccessorT;
933  typedef typename tree::ValueAccessor<const BoolTreeT> BoolTreeCAccessorT;
934  typedef typename tree::ValueAccessor<const FloatTreeT> FloatTreeCAccessorT;
935 
936  boost::scoped_ptr<FloatTreeCAccessorT> adaptivityAcc;
937  if (mAdaptivityGrid) {
938  adaptivityAcc.reset(new FloatTreeCAccessorT(mAdaptivityGrid->tree()));
939  }
940 
941  boost::scoped_ptr<BoolTreeAccessorT> seamMaskAcc;
942  boost::scoped_ptr<BoolTreeCAccessorT> topologyMaskAcc;
943  if (mRefData && mRefData->isValid()) {
944  seamMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
945  topologyMaskAcc.reset(new BoolTreeCAccessorT(*mRefData->mTopologyMaskTree));
946  }
947  const bool hasRefData = seamMaskAcc && topologyMaskAcc;
948 
949  const int LeafDim = BoolLeafT::DIM;
950  tree::ValueAccessor<const DistTreeT> distAcc(mDistTree);
951 
952  // Allocate reusable leaf buffers
953  BoolLeafT mask;
954  Vec3LeafT gradientBuffer;
955  Coord ijk, coord, end;
956 
957  for (size_t n = range.begin(); n != range.end(); ++n) {
958 
959  DistValueT adaptivity = mAdaptivity;
960  IntLeafT& auxLeaf = *mAuxLeafs[n];
961 
962  const Coord& origin = auxLeaf.getOrigin();
963  end[0] = origin[0] + LeafDim;
964  end[1] = origin[1] + LeafDim;
965  end[2] = origin[2] + LeafDim;
966 
967  mask.setValuesOff();
968 
969  // Mask off seam line adjacent voxels
970  if (hasRefData) {
971  const BoolLeafT* seamMask = seamMaskAcc->probeConstLeaf(origin);
972  if (seamMask != NULL) {
973  for (BoolCIterT it = seamMask->cbeginValueOn(); it; ++it) {
974  ijk = it.getCoord();
975  coord[0] = ijk[0] - (ijk[0] % 2);
976  coord[1] = ijk[1] - (ijk[1] % 2);
977  coord[2] = ijk[2] - (ijk[2] % 2);
978  mask.setActiveState(coord, true);
979  }
980  }
981  if (topologyMaskAcc->probeConstLeaf(origin) == NULL) {
982  adaptivity = mRefData->mInternalAdaptivity;
983  }
984  }
985 
986 
987  DistLeafT adaptivityLeaf(origin, adaptivity);
988 
989  if (mAdaptivityGrid) {
990  for (Index offset = 0; offset < DistLeafT::NUM_VALUES; ++offset) {
991  ijk = adaptivityLeaf.offset2globalCoord(offset);
992  Vec3d xyz = mAdaptivityGrid->transform().worldToIndex(
993  mTransform->indexToWorld(ijk));
994  DistValueT tmpA = DistValueT(adaptivityAcc->getValue(util::nearestCoord(xyz)));
995  adaptivityLeaf.setValueOnly(offset, tmpA * adaptivity);
996  }
997  }
998 
999  // Mask off ambiguous voxels
1000  for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1001  ijk = it.getCoord();
1002  coord[0] = ijk[0] - (ijk[0] % 2);
1003  coord[1] = ijk[1] - (ijk[1] % 2);
1004  coord[2] = ijk[2] - (ijk[2] % 2);
1005  if(mask.isValueOn(coord)) continue;
1006  mask.setActiveState(coord, isAmbiguous(distAcc, ijk, mIsovalue, 1));
1007  }
1008 
1009  int dim = 2;
1010  // Mask off topologically ambiguous 2x2x2 voxel sub-blocks
1011  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
1012  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
1013  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
1014  if (isNonManifold(distAcc, ijk, mIsovalue, dim)) {
1015  mask.setActiveState(ijk, true);
1016  }
1017  }
1018  }
1019  }
1020 
1021  // Compute the gradient for the remaining voxels
1022  gradientBuffer.setValuesOff();
1023 
1024  for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1025 
1026  ijk = it.getCoord();
1027  coord[0] = ijk[0] - (ijk[0] % dim);
1028  coord[1] = ijk[1] - (ijk[1] % dim);
1029  coord[2] = ijk[2] - (ijk[2] % dim);
1030  if(mask.isValueOn(coord)) continue;
1031 
1032  Vec3T norm(math::ISGradient<math::CD_2ND>::result(distAcc, ijk));
1033  // Normalize (Vec3's normalize uses isApproxEqual, which uses abs and does more work)
1034  DistValueT length = norm.length();
1035  if (length > DistValueT(1.0e-7)) {
1036  norm *= DistValueT(1.0) / length;
1037  }
1038  gradientBuffer.setValue(ijk, norm);
1039  }
1040 
1041  int regionId = 1, next_dim = dim << 1;
1042 
1043  // Process the first adaptivity level.
1044  for (ijk[0] = 0; ijk[0] < LeafDim; ijk[0] += dim) {
1045  coord[0] = ijk[0] - (ijk[0] % next_dim);
1046  for (ijk[1] = 0; ijk[1] < LeafDim; ijk[1] += dim) {
1047  coord[1] = ijk[1] - (ijk[1] % next_dim);
1048  for (ijk[2] = 0; ijk[2] < LeafDim; ijk[2] += dim) {
1049  coord[2] = ijk[2] - (ijk[2] % next_dim);
1050  adaptivity = adaptivityLeaf.getValue(ijk);
1051  if(mask.isValueOn(ijk) || !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
1052  mask.setActiveState(coord, true);
1053  continue;
1054  }
1055  mergeVoxels(auxLeaf, ijk, dim, regionId++);
1056  }
1057  }
1058  }
1059 
1060 
1061  // Process remaining adaptivity levels
1062  for (dim = 4; dim < LeafDim; dim = dim << 1) {
1063  next_dim = dim << 1;
1064  coord[0] = ijk[0] - (ijk[0] % next_dim);
1065  for (ijk[0] = origin[0]; ijk[0] < end[0]; ijk[0] += dim) {
1066  coord[1] = ijk[1] - (ijk[1] % next_dim);
1067  for (ijk[1] = origin[1]; ijk[1] < end[1]; ijk[1] += dim) {
1068  coord[2] = ijk[2] - (ijk[2] % next_dim);
1069  for (ijk[2] = origin[2]; ijk[2] < end[2]; ijk[2] += dim) {
1070  adaptivity = adaptivityLeaf.getValue(ijk);
1071  if (mask.isValueOn(ijk) || isNonManifold(distAcc, ijk, mIsovalue, dim) ||
1072  !isMergable(gradientBuffer, ijk, dim, adaptivity)) {
1073  mask.setActiveState(coord, true);
1074  continue;
1075  }
1076  mergeVoxels(auxLeaf, ijk, dim, regionId++);
1077  }
1078  }
1079  }
1080  }
1081 
1082  adaptivity = adaptivityLeaf.getValue(origin);
1083  if (!(mask.isValueOn(origin) || isNonManifold(distAcc, origin, mIsovalue, LeafDim))
1084  && isMergable(gradientBuffer, origin, LeafDim, adaptivity)) {
1085  mergeVoxels(auxLeaf, origin, LeafDim, regionId++);
1086  }
1087 
1088  if (mPointMask) {
1089  const BoolLeafT * ptnMaskLeaf = mPointMask->probeConstLeaf(auxLeaf.getOrigin());
1090 
1091  if (ptnMaskLeaf) { // remove unused regions
1092 
1093  for (IntIterT it = auxLeaf.beginValueOn(); it; ++it) {
1094  if (it.getValue() == 0 && !ptnMaskLeaf->isValueOn(it.pos())) {
1095  it.setValueOff();
1096  }
1097  }
1098 
1099  IntLeafT tmpLeaf(auxLeaf);
1100 
1101  while (tmpLeaf.onVoxelCount() > 0) {
1102 
1103  IntIterT it = tmpLeaf.beginValueOn();
1104  regionId = it.getValue();
1105  bool removeRegion = true;
1106  for (; it; ++it) {
1107  if (it.getValue() == regionId) {
1108  if (ptnMaskLeaf->isValueOn(it.pos())) {
1109  removeRegion = false;
1110  }
1111  it.setValueOff();
1112  }
1113  }
1114 
1115  if (removeRegion) {
1116  it = auxLeaf.beginValueOn();
1117  for (; it; ++it) {
1118  if (it.getValue() == regionId) {
1119  it.setValueOff();
1120  }
1121  }
1122  }
1123  }
1124  } else {
1125  mLeafRegionCount[n] = 0;
1126  continue;
1127  }
1128  }
1129 
1130  // Count unique regions
1131  size_t numVoxels = 0;
1132  IntLeafT tmpLeaf(auxLeaf);
1133  for (IntIterT it = tmpLeaf.beginValueOn(); it; ++it) {
1134  if(it.getValue() == 0) {
1135  it.setValueOff();
1136  ++numVoxels;
1137  }
1138  }
1139 
1140  while (tmpLeaf.onVoxelCount() > 0) {
1141  ++numVoxels;
1142  IntIterT it = tmpLeaf.beginValueOn();
1143  regionId = it.getValue();
1144  for (; it; ++it) {
1145  if (it.getValue() == regionId) it.setValueOff();
1146  }
1147  }
1148 
1149  mLeafRegionCount[n] = numVoxels;
1150  }
1151 }
1152 
1153 
1155 
1156 
1157 template <class DistTreeT>
1159 {
1160 public:
1161  typedef typename DistTreeT::ValueType DistValueT;
1163  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
1164 
1165  PointGen(
1166  const DistTreeT& distTree,
1167  const LeafPtrList<IntTreeT>& auxLeafs,
1168  std::vector<size_t>& leafIndices,
1169  const openvdb::math::Transform& xform,
1170  PointList& points,
1171  double iso = 0.0);
1172 
1173  PointGen(const PointGen<DistTreeT>&);
1174 
1175  void setRefData(const ReferenceData<DistTreeT>&);
1176 
1177  void runParallel();
1178  void runSerial();
1179 
1180  void operator()(const tbb::blocked_range<size_t>&) const;
1181 
1182 private:
1183  const DistTreeT& mDistTree;
1184  const LeafPtrList<IntTreeT>& mAuxLeafs;
1185  const std::vector<size_t>& mLeafIndices;
1186  const openvdb::math::Transform& mTransform;
1187  const PointList& mPoints;
1188  const double mIsovalue;
1189 
1190  const ReferenceData<DistTreeT>* mRefData;
1191 
1192  double root(double v0, double v1) const { return (mIsovalue - v0) / (v1 - v0); }
1193  int calcAvgPoint(DistTreeAccessorT&, const Coord&, openvdb::Vec3d&) const;
1194 };
1195 
1196 
1197 template <class DistTreeT>
1199  const DistTreeT& distTree,
1200  const LeafPtrList<IntTreeT>& auxLeafs,
1201  std::vector<size_t>& leafIndices,
1202  const openvdb::math::Transform& xform,
1203  PointList& points,
1204  double iso)
1205  : mDistTree(distTree)
1206  , mAuxLeafs(auxLeafs)
1207  , mLeafIndices(leafIndices)
1208  , mTransform(xform)
1209  , mPoints(points)
1210  , mIsovalue(iso)
1211  , mRefData(NULL)
1212 {
1213 }
1214 
1215 
1216 template <class DistTreeT>
1218  : mDistTree(rhs.mDistTree)
1219  , mAuxLeafs(rhs.mAuxLeafs)
1220  , mLeafIndices(rhs.mLeafIndices)
1221  , mTransform(rhs.mTransform)
1222  , mPoints(rhs.mPoints)
1223  , mIsovalue(rhs.mIsovalue)
1224  , mRefData(rhs.mRefData)
1225 {
1226 }
1227 
1228 
1229 template <class DistTreeT>
1230 void
1232  const ReferenceData<DistTreeT>& refData)
1233 {
1234  mRefData = &refData;
1235 }
1236 
1237 template <class DistTreeT>
1238 void
1240 {
1241  tbb::parallel_for(mAuxLeafs.getRange(), *this);
1242 }
1243 
1244 
1245 template <class DistTreeT>
1246 void
1248 {
1249  (*this)(mAuxLeafs.getRange());
1250 }
1251 
1252 
1253 template <class DistTreeT>
1254 void
1255 PointGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
1256 {
1257  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1258  typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
1259 
1260  typedef tree::ValueAccessor<BoolTreeT> BoolTreeAccessorT;
1261  typedef tree::ValueAccessor<Vec3sTreeT> Vec3sTreeAccessorT;
1262 
1263  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1264  typedef typename IntTreeT::LeafNodeType IntLeafT;
1265  typedef typename Vec3sTreeT::LeafNodeType Vec3sLeafT;
1266 
1267  boost::scoped_ptr<DistTreeAccessorT> refDistAcc;
1268  boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc, refSmoothMaskAcc;
1269  boost::scoped_ptr<Vec3sTreeAccessorT> refPtnAcc;
1270 
1271  if (mRefData && mRefData->isValid()) {
1272  refDistAcc.reset(new DistTreeAccessorT(*mRefData->mDistTree));
1273  refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mTopologyMaskTree));
1274  refSmoothMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSmoothingMaskTree));
1275  refPtnAcc.reset(new Vec3sTreeAccessorT(*mRefData->mSeamPointTree));
1276  }
1277 
1278 
1279  const bool hasRefData = refDistAcc && refMaskAcc;
1280  typename IntTreeT::LeafNodeType::ValueOnIter auxIter;
1281  DistTreeAccessorT distAcc(mDistTree);
1282 
1283  Coord ijk;
1284  openvdb::Vec3d avg, tmp;
1285 
1286  for (size_t n = range.begin(); n != range.end(); ++n) {
1287 
1288  size_t idx = mLeafIndices[n];
1289  IntLeafT& auxLeaf = *mAuxLeafs[n];
1290 
1291  BoolLeafT* maskLeaf = NULL;
1292  BoolLeafT* smoothMaskLeaf = NULL;
1293  Vec3sLeafT* ptnLeaf = NULL;
1294  if (hasRefData) {
1295  maskLeaf = refMaskAcc->probeLeaf(auxLeaf.getOrigin());
1296  smoothMaskLeaf = refSmoothMaskAcc->probeLeaf(auxLeaf.getOrigin());
1297  ptnLeaf = refPtnAcc->probeLeaf(auxLeaf.getOrigin());
1298  }
1299 
1300  for (auxIter = auxLeaf.beginValueOn(); auxIter; ++auxIter) {
1301 
1302  if(auxIter.getValue() == 0) {
1303 
1304  auxIter.setValue(idx);
1305  auxIter.setValueOff();
1306  ijk = auxIter.getCoord();
1307 
1308  if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1309 
1310  if (ptnLeaf && ptnLeaf->isValueOn(ijk)) {
1311  avg = ptnLeaf->getValue(ijk);
1312  } else {
1313  int e1 = calcAvgPoint(*refDistAcc.get(), ijk, avg);
1314 
1315  if (e1 != (XEDGE|YEDGE|ZEDGE)) {
1316  int e2 = calcAvgPoint(distAcc, ijk, tmp);
1317  if((e2 & (~e1)) != 0) smoothMaskLeaf->setValueOn(ijk);
1318  }
1319  }
1320  } else {
1321  calcAvgPoint(distAcc, ijk, avg);
1322  }
1323 
1324  openvdb::Vec3s& ptn = mPoints[idx];
1325  ptn[0] = float(avg[0]);
1326  ptn[1] = float(avg[1]);
1327  ptn[2] = float(avg[2]);
1328 
1329  ++idx;
1330  }
1331  }
1332 
1333  while(auxLeaf.onVoxelCount() > 0) {
1334 
1335  avg[0] = 0;
1336  avg[1] = 0;
1337  avg[2] = 0;
1338 
1339  auxIter = auxLeaf.beginValueOn();
1340  int regionId = auxIter.getValue(), points = 0;
1341 
1342  for (; auxIter; ++auxIter) {
1343  if(auxIter.getValue() == regionId) {
1344 
1345  auxIter.setValue(idx);
1346  auxIter.setValueOff();
1347  ijk = auxIter.getCoord();
1348 
1349  if (hasRefData && maskLeaf && maskLeaf->isValueOn(ijk)) {
1350  calcAvgPoint(*refDistAcc.get(), ijk, tmp);
1351  } else {
1352  calcAvgPoint(distAcc, ijk, tmp);
1353  }
1354 
1355  avg += tmp;
1356  ++points;
1357  }
1358  }
1359 
1360 
1361  if (points > 1) {
1362  double w = 1.0 / double(points);
1363  avg[0] *= w;
1364  avg[1] *= w;
1365  avg[2] *= w;
1366  }
1367 
1368  openvdb::Vec3s& ptn = mPoints[idx];
1369  ptn[0] = float(avg[0]);
1370  ptn[1] = float(avg[1]);
1371  ptn[2] = float(avg[2]);
1372  ++idx;
1373  }
1374  }
1375 }
1376 
1377 template <class DistTreeT>
1378 int
1379 PointGen<DistTreeT>::calcAvgPoint(DistTreeAccessorT& acc,
1380  const Coord& ijk, openvdb::Vec3d& avg) const
1381 {
1382  double values[8];
1383  bool signMask[8];
1384  Coord coord;
1385 
1386  // Sample corner values
1387  coord = ijk;
1388  values[0] = double(acc.getValue(coord)); // i, j, k
1389 
1390  coord[0] += 1;
1391  values[1] = double(acc.getValue(coord)); // i+1, j, k
1392 
1393  coord[2] += 1;
1394  values[2] = double(acc.getValue(coord)); // i+i, j, k+1
1395 
1396  coord[0] = ijk[0];
1397  values[3] = double(acc.getValue(coord)); // i, j, k+1
1398 
1399  coord[1] += 1; coord[2] = ijk[2];
1400  values[4] = double(acc.getValue(coord)); // i, j+1, k
1401 
1402  coord[0] += 1;
1403  values[5] = double(acc.getValue(coord)); // i+1, j+1, k
1404 
1405  coord[2] += 1;
1406  values[6] = double(acc.getValue(coord)); // i+1, j+1, k+1
1407 
1408  coord[0] = ijk[0];
1409  values[7] = double(acc.getValue(coord)); // i, j+1, k+1
1410 
1411  // init sign mask
1412  for (int n = 0; n < 8; ++n) signMask[n] = (values[n] < mIsovalue);
1413 
1414  int samples = 0, edgeFlags = 0;
1415  avg[0] = 0.0;
1416  avg[1] = 0.0;
1417  avg[2] = 0.0;
1418 
1419  if (signMask[0] != signMask[1]) { // Edged: 0 - 1
1420  avg[0] += root(values[0], values[1]);
1421  ++samples;
1422  edgeFlags |= XEDGE;
1423  }
1424 
1425  if (signMask[1] != signMask[2]) { // Edged: 1 - 2
1426  avg[0] += 1.0;
1427  avg[2] += root(values[1], values[2]);
1428  ++samples;
1429  edgeFlags |= ZEDGE;
1430  }
1431 
1432  if (signMask[3] != signMask[2]) { // Edged: 3 - 2
1433  avg[0] += root(values[3], values[2]);
1434  avg[2] += 1.0;
1435  ++samples;
1436  edgeFlags |= XEDGE;
1437  }
1438 
1439  if (signMask[0] != signMask[3]) { // Edged: 0 - 3
1440  avg[2] += root(values[0], values[3]);
1441  ++samples;
1442  edgeFlags |= ZEDGE;
1443  }
1444 
1445  if (signMask[4] != signMask[5]) { // Edged: 4 - 5
1446  avg[0] += root(values[4], values[5]);
1447  avg[1] += 1.0;
1448  ++samples;
1449  edgeFlags |= XEDGE;
1450  }
1451 
1452  if (signMask[5] != signMask[6]) { // Edged: 5 - 6
1453  avg[0] += 1.0;
1454  avg[1] += 1.0;
1455  avg[2] += root(values[5], values[6]);
1456  ++samples;
1457  edgeFlags |= ZEDGE;
1458  }
1459 
1460  if (signMask[7] != signMask[6]) { // Edged: 7 - 6
1461  avg[0] += root(values[7], values[6]);
1462  avg[1] += 1.0;
1463  avg[2] += 1.0;
1464  ++samples;
1465  edgeFlags |= XEDGE;
1466  }
1467 
1468  if (signMask[4] != signMask[7]) { // Edged: 4 - 7
1469  avg[1] += 1.0;
1470  avg[2] += root(values[4], values[7]);
1471  ++samples;
1472  edgeFlags |= ZEDGE;
1473  }
1474 
1475  if (signMask[0] != signMask[4]) { // Edged: 0 - 4
1476  avg[1] += root(values[0], values[4]);
1477  ++samples;
1478  edgeFlags |= YEDGE;
1479  }
1480 
1481  if (signMask[1] != signMask[5]) { // Edged: 1 - 5
1482  avg[0] += 1.0;
1483  avg[1] += root(values[1], values[5]);
1484  ++samples;
1485  edgeFlags |= YEDGE;
1486  }
1487 
1488  if (signMask[2] != signMask[6]) { // Edged: 2 - 6
1489  avg[0] += 1.0;
1490  avg[1] += root(values[2], values[6]);
1491  avg[2] += 1.0;
1492  ++samples;
1493  edgeFlags |= YEDGE;
1494  }
1495 
1496  if (signMask[3] != signMask[7]) { // Edged: 3 - 7
1497  avg[1] += root(values[3], values[7]);
1498  avg[2] += 1.0;
1499  ++samples;
1500  edgeFlags |= YEDGE;
1501  }
1502 
1503  if (samples > 1) {
1504  double w = 1.0 / double(samples);
1505  avg[0] *= w;
1506  avg[1] *= w;
1507  avg[2] *= w;
1508  }
1509 
1510  // offset by cell-origin
1511  avg[0] += double(ijk[0]);
1512  avg[1] += double(ijk[1]);
1513  avg[2] += double(ijk[2]);
1514 
1515  avg = mTransform.indexToWorld(avg);
1516 
1517  return edgeFlags;
1518 }
1519 
1520 
1522 
1523 
1525 {
1526  QuadMeshOp(): mIdx(0), mPolygonPool(NULL) {}
1527 
1528  void init(const size_t upperBound, PolygonPool& quadPool)
1529  {
1530  mPolygonPool = &quadPool;
1531  mPolygonPool->resetQuads(upperBound);
1532  mIdx = 0;
1533  }
1534 
1535  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
1536  {
1537  if (!reverse) {
1538  mPolygonPool->quad(mIdx) = verts;
1539  } else {
1540  Vec4I& quad = mPolygonPool->quad(mIdx);
1541  quad[0] = verts[3];
1542  quad[1] = verts[2];
1543  quad[2] = verts[1];
1544  quad[3] = verts[0];
1545  }
1546  mPolygonPool->quadFlags(mIdx) = flags;
1547  ++mIdx;
1548  }
1549 
1550  void done() {}
1551 
1552 private:
1553  size_t mIdx;
1554  PolygonPool* mPolygonPool;
1555 };
1556 
1557 
1559 {
1560  AdaptiveMeshOp(): mQuadIdx(0), mTriangleIdx(0), mPolygonPool(NULL), mTmpPolygonPool() {}
1561 
1562  void init(const size_t upperBound, PolygonPool& polygonPool)
1563  {
1564  mPolygonPool = &polygonPool;
1565 
1566  mTmpPolygonPool.resetQuads(upperBound);
1567  mTmpPolygonPool.resetTriangles(upperBound);
1568 
1569  mQuadIdx = 0;
1570  mTriangleIdx = 0;
1571  }
1572 
1573  void addPrim(const Vec4I& verts, bool reverse, char flags = 0)
1574  {
1575  if (verts[0] != verts[1] && verts[0] != verts[2] && verts[0] != verts[3]
1576  && verts[1] != verts[2] && verts[1] != verts[3] && verts[2] != verts[3]) {
1577  mTmpPolygonPool.quadFlags(mQuadIdx) = flags;
1578  addQuad(verts, reverse);
1579  } else if (
1580  verts[0] == verts[3] &&
1581  verts[1] != verts[2] &&
1582  verts[1] != verts[0] &&
1583  verts[2] != verts[0]) {
1584  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1585  addTriangle(verts[0], verts[1], verts[2], reverse);
1586  } else if (
1587  verts[1] == verts[2] &&
1588  verts[0] != verts[3] &&
1589  verts[0] != verts[1] &&
1590  verts[3] != verts[1]) {
1591  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1592  addTriangle(verts[0], verts[1], verts[3], reverse);
1593  } else if (
1594  verts[0] == verts[1] &&
1595  verts[2] != verts[3] &&
1596  verts[2] != verts[0] &&
1597  verts[3] != verts[0]) {
1598  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1599  addTriangle(verts[0], verts[2], verts[3], reverse);
1600  } else if (
1601  verts[2] == verts[3] &&
1602  verts[0] != verts[1] &&
1603  verts[0] != verts[2] &&
1604  verts[1] != verts[2]) {
1605  mTmpPolygonPool.triangleFlags(mTriangleIdx) = flags;
1606  addTriangle(verts[0], verts[1], verts[2], reverse);
1607  }
1608  }
1609 
1610 
1611  void done()
1612  {
1613  mPolygonPool->resetQuads(mQuadIdx);
1614  for (size_t i = 0; i < mQuadIdx; ++i) {
1615  mPolygonPool->quad(i) = mTmpPolygonPool.quad(i);
1616  mPolygonPool->quadFlags(i) = mTmpPolygonPool.quadFlags(i);
1617  }
1618  mTmpPolygonPool.clearQuads();
1619 
1620  mPolygonPool->resetTriangles(mTriangleIdx);
1621  for (size_t i = 0; i < mTriangleIdx; ++i) {
1622  mPolygonPool->triangle(i) = mTmpPolygonPool.triangle(i);
1623  mPolygonPool->triangleFlags(i) = mTmpPolygonPool.triangleFlags(i);
1624  }
1625  mTmpPolygonPool.clearTriangles();
1626  }
1627 
1628 private:
1629 
1630  void addQuad(const Vec4I& verts, bool reverse)
1631  {
1632  if (!reverse) {
1633  mTmpPolygonPool.quad(mQuadIdx) = verts;
1634  } else {
1635  Vec4I& quad = mTmpPolygonPool.quad(mQuadIdx);
1636  quad[0] = verts[3];
1637  quad[1] = verts[2];
1638  quad[2] = verts[1];
1639  quad[3] = verts[0];
1640  }
1641  ++mQuadIdx;
1642  }
1643 
1644  void addTriangle(unsigned v0, unsigned v1, unsigned v2, bool reverse)
1645  {
1646  Vec3I& prim = mTmpPolygonPool.triangle(mTriangleIdx);
1647 
1648  prim[1] = v1;
1649 
1650  if (!reverse) {
1651  prim[0] = v0;
1652  prim[2] = v2;
1653  } else {
1654  prim[0] = v2;
1655  prim[2] = v0;
1656  }
1657  ++mTriangleIdx;
1658  }
1659 
1660  size_t mQuadIdx, mTriangleIdx;
1661  PolygonPool *mPolygonPool;
1662  PolygonPool mTmpPolygonPool;
1663 };
1664 
1665 
1667 
1668 
1669 template<class DistTreeT, class MeshingOp>
1670 class MeshGen
1671 {
1672 public:
1673  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
1674  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
1675 
1676  MeshGen(const LeafCPtrList<CharTreeT>& edgeLeafs, const IntTreeT& auxTree, PolygonPoolList&);
1678 
1679  void setRefData(const ReferenceData<DistTreeT>&);
1680 
1681  void runParallel();
1682  void runSerial();
1683 
1684  void operator()(const tbb::blocked_range<size_t>&) const;
1685 
1686 private:
1687  const LeafCPtrList<CharTreeT>& mEdgeLeafs;
1688  const IntTreeT& mAuxTree;
1689  const PolygonPoolList& mPolygonPoolList;
1690  size_t mID;
1691  const ReferenceData<DistTreeT>* mRefData;
1692 };
1693 
1694 
1695 template<class DistTreeT, class MeshingOp>
1697  const IntTreeT& auxTree, PolygonPoolList& polygonPoolList)
1698  : mEdgeLeafs(edgeLeafs)
1699  , mAuxTree(auxTree)
1700  , mPolygonPoolList(polygonPoolList)
1701  , mRefData(NULL)
1702 {
1703 }
1704 
1705 
1706 template<class DistTreeT, class MeshingOp>
1708  : mEdgeLeafs(rhs.mEdgeLeafs)
1709  , mAuxTree(rhs.mAuxTree)
1710  , mPolygonPoolList(rhs.mPolygonPoolList)
1711  , mRefData(rhs.mRefData)
1712 {
1713 }
1714 
1715 
1716 template<class DistTreeT, class MeshingOp>
1717 void
1719  const ReferenceData<DistTreeT>& refData)
1720 {
1721  mRefData = &refData;
1722 }
1723 
1724 
1725 template<class DistTreeT, class MeshingOp>
1726 void
1728 {
1729  tbb::parallel_for(mEdgeLeafs.getRange(), *this);
1730 }
1731 
1732 
1733 template<class DistTreeT, class MeshingOp>
1734 void
1736 {
1737  (*this)(mEdgeLeafs.getRange());
1738 }
1739 
1740 
1741 template<class DistTreeT, class MeshingOp>
1742 void
1744  const tbb::blocked_range<size_t>& range) const
1745 {
1746  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
1747  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1748  typedef typename CharTreeT::LeafNodeType CharLeafT;
1749 
1750  typedef openvdb::tree::ValueAccessor<const CharTreeT> CharTreeAccessorT;
1751  typedef openvdb::tree::ValueAccessor<const IntTreeT> IntTreeAccessorT;
1752  typedef openvdb::tree::ValueAccessor<const BoolTreeT> BoolTreeAccessorT;
1753 
1754  boost::scoped_ptr<CharTreeAccessorT> refEdgeAcc;
1755  boost::scoped_ptr<BoolTreeAccessorT> refMaskAcc;
1756  if (mRefData && mRefData->isValid()) {
1757  refEdgeAcc.reset(new CharTreeAccessorT(*mRefData->mEdgeTree));
1758  refMaskAcc.reset(new BoolTreeAccessorT(*mRefData->mSeamMaskTree.get()));
1759  }
1760  const bool hasRefData = refEdgeAcc && refMaskAcc;
1761 
1762 
1763  typename CharTreeT::LeafNodeType::ValueOnCIter iter;
1764  IntTreeAccessorT auxAcc(mAuxTree);
1765 
1766  Coord ijk, coord;
1767  char refEdgeFlags, isSemLinePoly;
1768  const char isExteriorPoly[2] = {0, char(POLYFLAG_EXTERIOR)};
1769  openvdb::Vec4I quad;
1770  size_t edgeCount;
1771 
1772  MeshingOp mesher;
1773 
1774  for (size_t n = range.begin(); n != range.end(); ++n) {
1775 
1776  const Coord origin = mEdgeLeafs[n]->getOrigin();
1777 
1778  // Get an upper bound on the number of primitives.
1779  edgeCount = 0;
1780  iter = mEdgeLeafs[n]->cbeginValueOn();
1781  for (; iter; ++iter) {
1782  char edgeFlags = iter.getValue() >> 1;
1783  edgeCount += edgeFlags & 0x1;
1784 
1785  edgeFlags = edgeFlags >> 1;
1786  edgeCount += edgeFlags & 0x1;
1787 
1788  edgeFlags = edgeFlags >> 1;
1789  edgeCount += edgeFlags & 0x1;
1790  }
1791 
1792  mesher.init(edgeCount, mPolygonPoolList[n]);
1793 
1794  const CharLeafT* refEdgeLeaf = NULL;
1795  const BoolLeafT* refMaskLeaf = NULL;
1796 
1797  if (hasRefData) {
1798  refEdgeLeaf = refEdgeAcc->probeConstLeaf(origin);
1799  refMaskLeaf = refMaskAcc->probeConstLeaf(origin);
1800  }
1801 
1802  iter = mEdgeLeafs[n]->cbeginValueOn();
1803  for (; iter; ++iter) {
1804  ijk = iter.getCoord();
1805  const char& edgeFlags = iter.getValue();
1806 
1807  const bool isInside = edgeFlags & INSIDE;
1808 
1809  refEdgeFlags = 0;
1810  isSemLinePoly = 0;
1811  if (hasRefData) {
1812  if(refEdgeLeaf) refEdgeFlags = refEdgeLeaf->getValue(ijk);
1813  if (refMaskLeaf && refMaskLeaf->isValueOn(ijk)) {
1814  isSemLinePoly = char(POLYFLAG_FRACTURE_SEAM);
1815  }
1816  }
1817 
1818 
1819  int v0 = auxAcc.getValue(ijk);
1820 
1821  if (edgeFlags & XEDGE) {
1822 
1823  quad[0] = v0;
1824  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
1825  quad[1] = auxAcc.getValue(coord);
1826  coord[2] -= 1; // i, j-1, k-1
1827  quad[2] = auxAcc.getValue(coord);
1828  coord[1] = ijk[1]; // i, j, k-1
1829  quad[3] = auxAcc.getValue(coord);
1830 
1831  mesher.addPrim(quad, isInside,
1832  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & XEDGE)]));
1833  }
1834 
1835 
1836  if (edgeFlags & YEDGE) {
1837 
1838  quad[0] = v0;
1839  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1; // i, j, k-1
1840  quad[1] = auxAcc.getValue(coord);
1841  coord[0] -= 1; // i-1, j, k-1
1842  quad[2] = auxAcc.getValue(coord);
1843  coord[2] = ijk[2]; // i-1, j, k
1844  quad[3] = auxAcc.getValue(coord);
1845 
1846  mesher.addPrim(quad, isInside,
1847  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & YEDGE)]));
1848  }
1849 
1850  if (edgeFlags & ZEDGE) {
1851 
1852  quad[0] = v0;
1853  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]; // i, j-1, k
1854  quad[1] = auxAcc.getValue(coord);
1855  coord[0] -= 1; // i-1, j-1, k
1856  quad[2] = auxAcc.getValue(coord);
1857  coord[1] = ijk[1]; // i, j, k
1858  quad[3] = auxAcc.getValue(coord);
1859 
1860  mesher.addPrim(quad, !isInside,
1861  (isSemLinePoly | isExteriorPoly[bool(refEdgeFlags & ZEDGE)]));
1862  }
1863  }
1864 
1865  mesher.done();
1866  }
1867 }
1868 
1869 
1871 
1872 // Masking and mesh partitioning
1873 
1874 struct PartOp
1875 {
1876 
1877  PartOp(size_t leafCount, size_t partitions, size_t activePart)
1878  {
1879  size_t leafSegments = leafCount / partitions;
1880  mStart = leafSegments * activePart;
1881  mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
1882  }
1883 
1884  template <typename LeafNodeType>
1885  void operator()(LeafNodeType &leaf, size_t leafIndex) const
1886  {
1887  if (leafIndex < mStart || leafIndex >= mEnd) leaf.setValuesOff();
1888  }
1889 
1890 private:
1891  size_t mStart, mEnd;
1892 };
1893 
1894 template<typename SrcTreeT>
1895 class PartGen
1896 {
1897 public:
1899  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
1901 
1903 
1904 
1905  PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart);
1906 
1907  void run(bool threaded = true);
1908 
1909  BoolTreeT& tree() { return mTree; }
1910 
1911 
1913 
1914  PartGen(PartGen&, tbb::split);
1915  void operator()(const tbb::blocked_range<size_t>&);
1916  void join(PartGen& rhs) { mTree.merge(rhs.mTree); }
1917 
1918 private:
1919  const LeafManagerT& mLeafManager;
1920  BoolTreeT mTree;
1921  size_t mStart, mEnd;
1922 };
1923 
1924 template<typename SrcTreeT>
1925 PartGen<SrcTreeT>::PartGen(const LeafManagerT& leafs, size_t partitions, size_t activePart)
1926  : mLeafManager(leafs)
1927  , mTree(false)
1928 {
1929  size_t leafCount = leafs.leafCount();
1930  size_t leafSegments = leafCount / partitions;
1931  mStart = leafSegments * activePart;
1932  mEnd = activePart >= (partitions - 1) ? leafCount : mStart + leafSegments;
1933 }
1934 
1935 template<typename SrcTreeT>
1937  : mLeafManager(rhs.mLeafManager)
1938  , mTree(false)
1939  , mStart(rhs.mStart)
1940  , mEnd(rhs.mEnd)
1941 {
1942 }
1943 
1944 
1945 template<typename SrcTreeT>
1946 void
1948 {
1949  if (threaded) {
1950  tbb::parallel_reduce(mLeafManager.getRange(), *this);
1951  } else {
1952  (*this)(mLeafManager.getRange());
1953  }
1954 }
1955 
1956 
1957 template<typename SrcTreeT>
1958 void
1959 PartGen<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
1960 {
1961  Coord ijk;
1962  BoolAccessorT acc(mTree);
1963 
1964  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
1965  typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
1966 
1967  for (size_t n = range.begin(); n != range.end(); ++n) {
1968  if (n < mStart || n >= mEnd) continue;
1969  BoolLeafT* leaf = acc.touchLeaf(mLeafManager.leaf(n).getOrigin());
1970  leaf->topologyUnion(mLeafManager.leaf(n));
1971  }
1972 }
1973 
1974 
1975 template<typename SrcTreeT>
1976 class MaskGen
1977 {
1978 public:
1980  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
1984 
1985 
1987 
1988 
1989  MaskGen(const BoolGridT& mask, const LeafManagerT& srcLeafs,
1990  const math::Transform& srcXForm, bool invertMask);
1991 
1992  void run(bool threaded = true);
1993 
1994  BoolTreeT& tree() { return mTree; }
1995 
1996 
1998 
1999  MaskGen(MaskGen&, tbb::split);
2000 
2001  void operator()(const tbb::blocked_range<size_t>&);
2002 
2003  void join(MaskGen& rhs) { mTree.merge(rhs.mTree); }
2004 
2005 private:
2006 
2007  const BoolGridT& mMask;
2008  const LeafManagerT& mLeafManager;
2009  const math::Transform& mSrcXForm;
2010  bool mInvertMask;
2011  BoolTreeT mTree;
2012 };
2013 
2014 template<typename SrcTreeT>
2016  const math::Transform& srcXForm, bool invertMask)
2017  : mMask(mask)
2018  , mLeafManager(srcLeafs)
2019  , mSrcXForm(srcXForm)
2020  , mInvertMask(invertMask)
2021  , mTree(false)
2022 {
2023 }
2024 
2025 template<typename SrcTreeT>
2027  : mMask(rhs.mMask)
2028  , mLeafManager(rhs.mLeafManager)
2029  , mSrcXForm(rhs.mSrcXForm)
2030  , mInvertMask(rhs.mInvertMask)
2031  , mTree(false)
2032 {
2033 }
2034 
2035 
2036 template<typename SrcTreeT>
2037 void
2039 {
2040  if (threaded) {
2041  tbb::parallel_reduce(mLeafManager.getRange(), *this);
2042  } else {
2043  (*this)(mLeafManager.getRange());
2044  }
2045 }
2046 
2047 
2048 template<typename SrcTreeT>
2049 void
2050 MaskGen<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
2051 {
2052  Coord ijk;
2053  Vec3d xyz;
2054  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2055  const math::Transform& maskXForm = mMask.transform();
2056  tree::ValueAccessor<const BoolTreeT> maskAcc(mMask.tree());
2057  tree::ValueAccessor<BoolTreeT> acc(mTree);
2058 
2059  typename SrcTreeT::LeafNodeType::ValueOnCIter iter;
2060  for (size_t n = range.begin(); n != range.end(); ++n) {
2061 
2062  ijk = mLeafManager.leaf(n).getOrigin();
2063  BoolLeafT* leaf = new BoolLeafT(ijk, false);
2064  bool addLeaf = false;
2065 
2066  if (maskXForm == mSrcXForm) {
2067 
2068  leaf->topologyUnion(mLeafManager.leaf(n));
2069  const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
2070  if (maskLeaf) {
2071  leaf->topologyDifference(*maskLeaf, false);
2072  addLeaf = !leaf->isEmpty();
2073  } else {
2074  addLeaf = maskAcc.isValueOn(ijk) != mInvertMask;
2075  }
2076 
2077  } else {
2078  for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2079  ijk = iter.getCoord();
2080  xyz = maskXForm.worldToIndex(mSrcXForm.indexToWorld(ijk));
2081  if(maskAcc.isValueOn(util::nearestCoord(xyz)) != mInvertMask) {
2082  leaf->setValueOn(iter.pos());
2083  addLeaf = true;
2084  }
2085  }
2086  }
2087 
2088  if (addLeaf) acc.addLeaf(leaf);
2089  else delete leaf;
2090  }
2091 }
2092 
2093 template<typename IntTreeT>
2095 {
2097  typedef typename IntTreeT::LeafNodeType IntLeafT;
2098 
2099  BoundaryMergeOp(IntTreeT& aux, IntTreeT& boundaryAux)
2100  : mAuxAcc(aux), mBoundaryAuxAcc(boundaryAux) {}
2101 
2102  template <typename LeafNodeType>
2103  void operator()(LeafNodeType &leaf, size_t/*leafIndex*/) const
2104  {
2105  Coord ijk = leaf.getOrigin();
2106 
2107  IntLeafT* auxLeaf = const_cast<IntLeafT*>(mAuxAcc.probeLeaf(ijk));
2108  const IntLeafT* boundaryLeaf = mBoundaryAuxAcc.probeConstLeaf(ijk);
2109 
2110  if (auxLeaf && boundaryLeaf) auxLeaf->merge(*boundaryLeaf);
2111  }
2112 
2113 private:
2114  IntAccessorT mAuxAcc, mBoundaryAuxAcc;
2115 };
2116 
2117 template<typename SrcTreeT>
2119 {
2120 public:
2121  typedef typename SrcTreeT::template ValueConverter<int>::Type IntTreeT;
2122  typedef typename SrcTreeT::template ValueConverter<bool>::Type BoolTreeT;
2125 
2127 
2129 
2130 
2131  BoundaryMaskGen(const LeafManagerT& auxLeafs, const SrcTreeT&, const BoolTreeT&);
2132 
2133  void run(bool threaded = true);
2134 
2135  BoolTreeT& tree() { return mTree; }
2136 
2138 
2139  BoundaryMaskGen(BoundaryMaskGen&, tbb::split);
2140 
2141  void operator()(const tbb::blocked_range<size_t>&);
2142 
2143  void join(BoundaryMaskGen& rhs) { mTree.merge(rhs.mTree); }
2144 
2145 private:
2146 
2147  const LeafManagerT& mLeafManager;
2148  const SrcTreeT& mSrcTree;
2149  const BoolTreeT& mMaskTree;
2150  BoolTreeT mTree;
2151 };
2152 
2153 template<typename SrcTreeT>
2155  const LeafManagerT& auxLeafs, const SrcTreeT& srcTree, const BoolTreeT& boolTree)
2156  : mLeafManager(auxLeafs)
2157  , mSrcTree(srcTree)
2158  , mMaskTree(boolTree)
2159  , mTree(false)
2160 {
2161 }
2162 
2163 template<typename SrcTreeT>
2165  : mLeafManager(rhs.mLeafManager)
2166  , mSrcTree(rhs.mSrcTree)
2167  , mMaskTree(rhs.mMaskTree)
2168  , mTree(false)
2169 {
2170 }
2171 
2172 template<typename SrcTreeT>
2173 void
2175 {
2176  if (threaded) {
2177  tbb::parallel_reduce(mLeafManager.getRange(), *this);
2178  } else {
2179  (*this)(mLeafManager.getRange());
2180  }
2181 }
2182 
2183 
2184 template<typename SrcTreeT>
2185 void
2186 BoundaryMaskGen<SrcTreeT>::operator()(const tbb::blocked_range<size_t>& range)
2187 {
2188  Coord ijk;
2189  typedef typename BoolTreeT::LeafNodeType BoolLeafT;
2190  typedef typename SrcTreeT::LeafNodeType SrcLeafT;
2191 
2192  tree::ValueAccessor<const SrcTreeT> srcAcc(mSrcTree);
2193  tree::ValueAccessor<const BoolTreeT> maskAcc(mMaskTree);
2194  tree::ValueAccessor<BoolTreeT> acc(mTree);
2195 
2196  typename IntTreeT::LeafNodeType::ValueOnCIter iter;
2197  for (size_t n = range.begin(); n != range.end(); ++n) {
2198 
2199  ijk = mLeafManager.leaf(n).getOrigin();
2200 
2201  const SrcLeafT* srcLeaf = srcAcc.probeConstLeaf(ijk);
2202  if (!srcLeaf) continue;
2203 
2204  BoolLeafT* leaf = new BoolLeafT(ijk, false);
2205  leaf->topologyUnion(*srcLeaf);
2206 
2207  const BoolLeafT* maskLeaf = maskAcc.probeConstLeaf(ijk);
2208  if (maskLeaf) {
2209 
2210  leaf->topologyDifference(*maskLeaf, false);
2211 
2212  if (!leaf->isEmpty()) acc.addLeaf(leaf);
2213  else delete leaf;
2214 
2215  } else {
2216  acc.addLeaf(leaf);
2217  }
2218  }
2219 }
2220 
2221 
2223 
2224 
2225 template<
2226  typename SrcTreeT,
2227  typename LeafManagerT = tree::LeafManager<const SrcTreeT>,
2228  typename AuxDataT = int>
2230 {
2231 public:
2232  typedef typename SrcTreeT::ValueType SrcValueT;
2233  typedef typename SrcTreeT::template ValueConverter<char>::Type CharTreeT;
2234  typedef typename SrcTreeT::template ValueConverter<AuxDataT>::Type AuxTreeT;
2238 
2240 
2241 
2242  AuxiliaryData(const SrcTreeT&, const LeafManagerT&, double iso, bool extraCheck = false);
2243 
2244  void run(bool threaded = true);
2245 
2246  typename CharTreeT::Ptr edgeTree() const { return mEdgeTree; }
2247  typename AuxTreeT::Ptr auxTree() const { return mAuxTree; }
2248 
2249 
2251 
2252  AuxiliaryData(AuxiliaryData&, tbb::split);
2253 
2254  void operator()(const tbb::blocked_range<size_t>&);
2255 
2256  void join(const AuxiliaryData& rhs)
2257  {
2258  mEdgeTree->merge(*rhs.mEdgeTree);
2259  mAuxTree->merge(*rhs.mAuxTree);
2260  }
2261 
2262 private:
2263 
2264  int edgeCheck(const Coord& ijk, const bool thisInside);
2265 
2266  const LeafManagerT& mLeafManager;
2267 
2268  const SrcTreeT& mSourceTree;
2269  SrcAccessorT mSourceAccessor;
2270 
2271  typename CharTreeT::Ptr mEdgeTree;
2272  CharAccessorT mEdgeAccessor;
2273 
2274  typename AuxTreeT::Ptr mAuxTree;
2275  AuxAccessorT mAuxAccessor;
2276 
2277  const double mIsovalue;
2278  const bool mExtraCheck;
2279 };
2280 
2281 template<typename SrcTreeT, typename LeafManagerT, typename AuxDataT>
2283  const LeafManagerT& leafs, double iso, bool extraCheck)
2284  : mLeafManager(leafs)
2285  , mSourceTree(tree)
2286  , mSourceAccessor(mSourceTree)
2287  , mEdgeTree(new CharTreeT(0))
2288  , mEdgeAccessor(*mEdgeTree)
2289  , mAuxTree(new AuxTreeT(AuxDataT(0)))
2290  , mAuxAccessor(*mAuxTree)
2291  , mIsovalue(iso)
2292  , mExtraCheck(extraCheck)
2293 {
2294 }
2295 
2296 template<typename SrcTreeT, typename LeafManagerT, typename AuxDataT>
2298  : mLeafManager(rhs.mLeafManager)
2299  , mSourceTree(rhs.mSourceTree)
2300  , mSourceAccessor(mSourceTree)
2301  , mEdgeTree(new CharTreeT(0))
2302  , mEdgeAccessor(*mEdgeTree)
2303  , mAuxTree(new AuxTreeT(AuxDataT(0)))
2304  , mAuxAccessor(*mAuxTree)
2305  , mIsovalue(rhs.mIsovalue)
2306  , mExtraCheck(rhs.mExtraCheck)
2307 {
2308 }
2309 
2310 
2311 template<typename SrcTreeT, typename LeafManagerT, typename AuxDataT>
2312 void
2314 {
2315  if (threaded) {
2316  tbb::parallel_reduce(mLeafManager.getRange(), *this);
2317  } else {
2318  (*this)(mLeafManager.getRange());
2319  }
2320 }
2321 
2322 
2323 template<typename SrcTreeT, typename LeafManagerT, typename AuxDataT>
2324 void
2325 AuxiliaryData<SrcTreeT, LeafManagerT, AuxDataT>::operator()(const tbb::blocked_range<size_t>& range)
2326 {
2327  typename LeafManagerT::TreeType::LeafNodeType::ValueOnCIter iter;
2328  Coord ijk;
2329  bool thisInside;
2330  int edgeFlags;
2331  SrcValueT val;
2332 
2333  if (!mExtraCheck) {
2334  for (size_t n = range.begin(); n != range.end(); ++n) {
2335  for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2336 
2337  ijk = iter.getCoord();
2338 
2339  thisInside = mSourceAccessor.getValue(ijk) < mIsovalue;
2340  edgeFlags = edgeCheck(ijk, thisInside);
2341 
2342  if (edgeFlags != 0) {
2343  edgeFlags |= int(thisInside);
2344  mEdgeAccessor.setValue(ijk, char(edgeFlags));
2345  }
2346  }
2347  }
2348  } else {
2349  for (size_t n = range.begin(); n != range.end(); ++n) {
2350  for (iter = mLeafManager.leaf(n).cbeginValueOn(); iter; ++iter) {
2351 
2352  ijk = iter.getCoord();
2353 
2354  thisInside = mSourceAccessor.getValue(ijk) < mIsovalue;
2355  edgeFlags = edgeCheck(ijk, thisInside);
2356 
2357  if (edgeFlags != 0) {
2358  edgeFlags |= int(thisInside);
2359  mEdgeAccessor.setValue(ijk, char(edgeFlags));
2360  }
2361 
2362  --ijk[0];
2363  if (!mSourceAccessor.probeValue(ijk, val)) {
2364  thisInside = val < mIsovalue;
2365  edgeFlags = edgeCheck(ijk, thisInside);
2366 
2367  if (edgeFlags != 0) {
2368  edgeFlags |= int(thisInside);
2369  mEdgeAccessor.setValue(ijk, char(edgeFlags));
2370  }
2371  }
2372 
2373  ++ijk[0];
2374  --ijk[1];
2375  if (!mSourceAccessor.probeValue(ijk, val)) {
2376  thisInside = val < mIsovalue;
2377  edgeFlags = edgeCheck(ijk, thisInside);
2378 
2379  if (edgeFlags != 0) {
2380  edgeFlags |= int(thisInside);
2381  mEdgeAccessor.setValue(ijk, char(edgeFlags));
2382  }
2383  }
2384 
2385  ++ijk[1];
2386  --ijk[2];
2387  if (!mSourceAccessor.probeValue(ijk, val)) {
2388  thisInside = val < mIsovalue;
2389  edgeFlags = edgeCheck(ijk, thisInside);
2390 
2391  if (edgeFlags != 0) {
2392  edgeFlags |= int(thisInside);
2393  mEdgeAccessor.setValue(ijk, char(edgeFlags));
2394  }
2395  }
2396  }
2397  }
2398  }
2399 }
2400 
2401 template<typename SrcTreeT, typename LeafManagerT, typename AuxDataT>
2402 inline int
2403 AuxiliaryData<SrcTreeT, LeafManagerT, AuxDataT>::edgeCheck(const Coord& ijk, const bool thisInside)
2404 {
2405  int edgeFlags = 0;
2406  Coord n_ijk, coord;
2407 
2408  // Eval upwind x-edge
2409  n_ijk = ijk; ++n_ijk[0];
2410  bool otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2411  if (otherInside != thisInside) {
2412 
2413  edgeFlags = XEDGE;
2414 
2415  mAuxAccessor.setActiveState(ijk, true);
2416 
2417  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2418  mAuxAccessor.setActiveState(coord, true);
2419 
2420  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
2421  mAuxAccessor.setActiveState(coord, true);
2422 
2423  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2424  mAuxAccessor.setActiveState(coord, true);
2425  }
2426 
2427  // Eval upwind y-edge
2428  n_ijk[0] = ijk[0]; ++n_ijk[1];
2429  otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2430  if (otherInside != thisInside) {
2431 
2432  edgeFlags |= YEDGE;
2433 
2434  mAuxAccessor.setActiveState(ijk, true);
2435 
2436  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2437  mAuxAccessor.setActiveState(coord, true);
2438 
2439  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2440  mAuxAccessor.setActiveState(coord, true);
2441 
2442  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2443  mAuxAccessor.setActiveState(coord, true);
2444  }
2445 
2446  // Eval upwind z-edge
2447  n_ijk[1] = ijk[1]; ++n_ijk[2];
2448  otherInside = (mSourceAccessor.getValue(n_ijk) < mIsovalue);
2449  if (otherInside != thisInside) {
2450 
2451  edgeFlags |= ZEDGE;
2452 
2453  mAuxAccessor.setActiveState(ijk, true);
2454 
2455  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2456  mAuxAccessor.setActiveState(coord, true);
2457 
2458  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2459  mAuxAccessor.setActiveState(coord, true);
2460 
2461  coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2462  mAuxAccessor.setActiveState(coord, true);
2463  }
2464  return edgeFlags;
2465 }
2466 
2467 
2468 
2469 
2471 
2472 
2473 template <class DistTreeT>
2475 {
2476 public:
2477  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
2478  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
2481 
2482  SeamMaskGen(LeafPtrList<BoolTreeT>& seamMaskLeafs,
2483  const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree);
2484 
2486 
2487  void runParallel();
2488  void runSerial();
2489 
2490  void operator()(const tbb::blocked_range<size_t>&) const;
2491 
2492 private:
2493  LeafPtrList<BoolTreeT>& mSeamMaskLeafs;
2494  const BoolTreeT& mTopologyMaskTree;
2495  BoolTreeAccessorT mTopologyMaskAcc;
2496  const IntTreeT& mAuxTree;
2497  IntTreeAccessorT mAuxAcc;
2498 };
2499 
2500 
2501 template <class DistTreeT>
2503  const BoolTreeT& topologyMaskTree, const IntTreeT& auxTree)
2504  : mSeamMaskLeafs(seamMaskLeafs)
2505  , mTopologyMaskTree(topologyMaskTree)
2506  , mTopologyMaskAcc(mTopologyMaskTree)
2507  , mAuxTree(auxTree)
2508  , mAuxAcc(mAuxTree)
2509 {
2510 }
2511 
2512 template <class DistTreeT>
2514  : mSeamMaskLeafs(rhs.mSeamMaskLeafs)
2515  , mTopologyMaskTree(rhs.mTopologyMaskTree)
2516  , mTopologyMaskAcc(mTopologyMaskTree)
2517  , mAuxTree(rhs.mAuxTree)
2518  , mAuxAcc(mAuxTree)
2519 {
2520 }
2521 
2522 template <class DistTreeT>
2523 void
2525 {
2526  tbb::parallel_for(mSeamMaskLeafs.getRange(), *this);
2527 }
2528 
2529 template <class DistTreeT>
2530 void
2532 {
2533  (*this)(mSeamMaskLeafs.getRange());
2534 }
2535 
2536 template <class DistTreeT>
2537 void
2538 SeamMaskGen<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
2539 {
2540  typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
2541  Coord ijk, n_ijk;
2542  for (size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
2543  ValueOnIterT it = mSeamMaskLeafs[leafIdx]->beginValueOn();
2544  for (; it; ++it) {
2545  ijk = it.getCoord();
2546  if (!mTopologyMaskAcc.isValueOn(ijk)) {
2547  it.setValueOff();
2548  } else {
2549  bool turnOff = true;
2550  for (size_t n = 0; n < 6; ++n) {
2551  n_ijk = ijk + util::COORD_OFFSETS[n];
2552  if (!mAuxTree.isValueOn(n_ijk) && mTopologyMaskAcc.isValueOn(n_ijk)) {
2553  turnOff = false;
2554  break;
2555  }
2556  }
2557  if (turnOff) it.setValueOff();
2558  }
2559  }
2560  }
2561 }
2562 
2564 
2565 
2566 template <class DistTreeT>
2568 {
2569 public:
2570  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
2571  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
2574 
2575  EdgeSmooth(
2576  LeafPtrList<BoolTreeT>& leafs,
2577  const BoolTreeT& edgeMaskTree,
2578  const IntTreeT& auxTree,
2579  PointList& points,
2580  const math::Transform& xform);
2581 
2583 
2584  void runParallel(const size_t iterations);
2585  void runSerial(const size_t iterations);
2586 
2587  void operator()(const tbb::blocked_range<size_t>&) const;
2588 
2589 private:
2590  LeafPtrList<BoolTreeT>& mLeafs;
2591  const BoolTreeT& mEdgeMaskTree;
2592  const IntTreeT& mAuxTree;
2593  PointList& mPoints;
2594  const math::Transform& mTransform;
2595 
2596  bool pointInAABB(const Vec3s& p, const Vec3s& bmin, const Vec3s& bmax) const
2597  {
2598  for (int i = 0; i < 3; ++i) {
2599  if (p[i] < bmin[i] || p[i] > bmax[i]) {
2600  return false;
2601  }
2602  }
2603  return true;
2604  }
2605 
2606 };
2607 
2608 
2609 template <class DistTreeT>
2611  LeafPtrList<BoolTreeT>& leafs,
2612  const BoolTreeT& edgeMaskTree,
2613  const IntTreeT& auxTree,
2614  PointList& points,
2615  const math::Transform& xform)
2616  : mLeafs(leafs)
2617  , mEdgeMaskTree(edgeMaskTree)
2618  , mAuxTree(auxTree)
2619  , mPoints(points)
2620  , mTransform(xform)
2621 {
2622 }
2623 
2624 template <class DistTreeT>
2626  : mLeafs(rhs.mLeafs)
2627  , mEdgeMaskTree(rhs.mEdgeMaskTree)
2628  , mAuxTree(rhs.mAuxTree)
2629  , mPoints(rhs.mPoints)
2630  , mTransform(rhs.mTransform)
2631 {
2632 }
2633 
2634 template <class DistTreeT>
2635 void
2636 EdgeSmooth<DistTreeT>::runParallel(const size_t iterations)
2637 {
2638  for (size_t i = 0; i < iterations; ++i) {
2639  tbb::parallel_for(mLeafs.getRange(), *this);
2640  }
2641 }
2642 
2643 template <class DistTreeT>
2644 void
2645 EdgeSmooth<DistTreeT>::runSerial(const size_t iterations)
2646 {
2647  for (size_t i = 0; i < iterations; ++i) {
2648  (*this)(mLeafs.getRange());
2649  }
2650 }
2651 
2652 template <class DistTreeT>
2653 void
2654 EdgeSmooth<DistTreeT>::operator()(const tbb::blocked_range<size_t>& range) const
2655 {
2656  typedef typename BoolTreeT::LeafNodeType::ValueOnIter ValueOnIterT;
2657  BoolTreeAccessorT maskAcc(mEdgeMaskTree);
2658  IntTreeAccessorT auxAcc(mAuxTree);
2659 
2660  float dx = float(mTransform.voxelSize()[0]);
2661  openvdb::Vec3s avg, bmin, bmax;
2662  Coord ijk, nijk;
2663  int count;
2664 
2665  for (size_t leafIdx = range.begin(); leafIdx != range.end(); ++leafIdx) {
2666  typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = mLeafs[leafIdx]->beginValueOn();
2667  for ( ; valueIt; ++valueIt) {
2668 
2669  ijk = valueIt.getCoord();
2670  openvdb::Vec3s& ptn = mPoints[auxAcc.getValue(ijk)];
2671 
2672  avg = ptn;
2673  count = 1;
2674  for (int n = 0; n < 26; ++n) {
2675  nijk = ijk + util::COORD_OFFSETS[n];
2676  if (maskAcc.isValueOn(nijk)) {
2677  avg += mPoints[auxAcc.getValue(nijk)];
2678  ++count;
2679  }
2680  }
2681 
2682  if (count > 2) {
2683  avg *= (1.0 / float(count));
2684 
2685  // Constrain to current cell
2686  bmin = openvdb::Vec3s(mTransform.indexToWorld(ijk));
2687  bmax = bmin + dx;
2688 
2689  bool inCell = false;
2690  for (int i = 0; i < 10; ++i) {
2691 
2692  inCell = pointInAABB(avg, bmin, bmax);
2693 
2694  if (inCell) {
2695  break;
2696  }
2697 
2698  avg += ptn;
2699  avg *= 0.5;
2700  }
2701 
2702  if (inCell) ptn = avg;
2703 
2704  /*Vec3d a = mTransform.worldToIndex(ptn);
2705  Vec3d b = mTransform.worldToIndex(avg);
2706 
2707  bmin[0] = float(ijk[0]);
2708  bmin[1] = float(ijk[1]);
2709  bmin[2] = float(ijk[2]);
2710  bmax = bmin + 1.0;
2711 
2712  if(!pointInAABB(a, bmin, bmax)) continue;
2713 
2714  while (true) {
2715  if(pointInAABB(b, bmin, bmax)) break;
2716 
2717  b += a;
2718  b *= 0.5;
2719  }
2720 
2721  ptn = mTransform.indexToWorld(b);*/
2722  }
2723  }
2724  }
2725 }
2726 
2727 
2729 
2730 
2731 template<class CharAccessorT, typename AuxAccessorT>
2733 {
2734 public:
2735  AuxDataGenerator(CharAccessorT& edgeAcc, AuxAccessorT& auxAcc)
2736  : mEdgeAcc(edgeAcc), mAuxAcc(auxAcc) {}
2737 
2738 
2739  void setXEdge(char edgeFlags, const Coord& ijk)
2740  {
2741  mEdgeAcc.setValue(ijk, edgeFlags | XEDGE);
2742 
2743  mAuxAcc.setActiveState(ijk, true);
2744 
2745  Coord coord = ijk;
2746  coord[1] = ijk[1]-1;
2747  mAuxAcc.setActiveState(coord, true);
2748 
2749  coord[0] = ijk[0]; coord[1] = ijk[1]-1; coord[2] = ijk[2]-1;
2750  mAuxAcc.setActiveState(coord, true);
2751 
2752  coord[0] = ijk[0]; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2753  mAuxAcc.setActiveState(coord, true);
2754  }
2755 
2756  void setYEdge(char edgeFlags, const Coord& ijk)
2757  {
2758  mEdgeAcc.setValue(ijk, edgeFlags | YEDGE);
2759 
2760  mAuxAcc.setActiveState(ijk, true);
2761 
2762  Coord coord = ijk;
2763  coord[2] = ijk[2]-1;
2764  mAuxAcc.setActiveState(coord, true);
2765 
2766  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2767  mAuxAcc.setActiveState(coord, true);
2768 
2769  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2]-1;
2770  mAuxAcc.setActiveState(coord, true);
2771  }
2772 
2773  void setZEdge(char edgeFlags, const Coord& ijk)
2774  {
2775  mEdgeAcc.setValue(ijk, edgeFlags | ZEDGE);
2776 
2777  mAuxAcc.setActiveState(ijk, true);
2778 
2779  Coord coord = ijk;
2780  coord[1] = ijk[1]-1;
2781  mAuxAcc.setActiveState(coord, true);
2782 
2783  coord[0] = ijk[0]-1; coord[1] = ijk[1]; coord[2] = ijk[2];
2784  mAuxAcc.setActiveState(coord, true);
2785 
2786  coord[0] = ijk[0]-1; coord[1] = ijk[1]-1; coord[2] = ijk[2];
2787  mAuxAcc.setActiveState(coord, true);
2788  }
2789 
2790 private:
2791  CharAccessorT& mEdgeAcc;
2792  AuxAccessorT& mAuxAcc;
2793 };
2794 
2795 
2797 
2798 
2799 template<class DistTreeT, class AuxTreeT, class CharTreeT>
2800 inline void
2802  const DistTreeT& distTree, CharTreeT& edgeTree, AuxTreeT& auxTree,
2803  double iso)
2804 {
2805  typedef tree::ValueAccessor<const DistTreeT> DistAccessorT;
2806  typedef tree::ValueAccessor<AuxTreeT> AuxTreeAccessorT;
2807  typedef tree::ValueAccessor<CharTreeT> CharTreeAccessorT;
2808 
2809  typename DistTreeT::ValueType isoValue = typename DistTreeT::ValueType(iso);
2810 
2811  DistAccessorT distAcc(distTree);
2812  CharTreeAccessorT edgeAcc(edgeTree);
2813  AuxTreeAccessorT auxAcc(auxTree);
2814 
2816 
2817  Coord ijk, nijk;
2818  typename DistTreeT::ValueType value;
2819  CoordBBox bbox;
2820  bool processTileFace;
2821 
2822  typename DistTreeT::ValueOnCIter tileIter(distTree);
2823  tileIter.setMaxDepth(DistTreeT::ValueOnCIter::LEAF_DEPTH - 1);
2824 
2825  for ( ; tileIter; ++tileIter) {
2826  tileIter.getBoundingBox(bbox);
2827 
2828  const bool thisInside = (distAcc.getValue(bbox.min()) < isoValue);
2829  const int thisDepth = distAcc.getValueDepth(bbox.min());
2830 
2831 
2832  // Eval x-edges
2833 
2834  ijk = bbox.max();
2835  nijk = ijk;
2836  ++nijk[0];
2837 
2838  processTileFace = true;
2839  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2840  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2841  }
2842 
2843  if (processTileFace) {
2844  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2845  nijk[1] = ijk[1];
2846  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2847  nijk[2] = ijk[2];
2848  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2849  auxData.setXEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2850  }
2851  }
2852  }
2853  }
2854 
2855  ijk = bbox.min();
2856  --ijk[0];
2857 
2858  processTileFace = true;
2859  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2860  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2861  }
2862 
2863  if (processTileFace) {
2864  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2865  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2866  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2867  auxData.setXEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2868  }
2869  }
2870  }
2871  }
2872 
2873 
2874  // Eval y-edges
2875 
2876  ijk = bbox.max();
2877  nijk = ijk;
2878  ++nijk[1];
2879 
2880  processTileFace = true;
2881  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2882  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2883  }
2884 
2885  if (processTileFace) {
2886  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2887  nijk[0] = ijk[0];
2888  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2889  nijk[2] = ijk[2];
2890 
2891  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2892  auxData.setYEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2893  }
2894  }
2895  }
2896  }
2897 
2898 
2899  ijk = bbox.min();
2900  --ijk[1];
2901 
2902  processTileFace = true;
2903  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2904  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2905  }
2906 
2907  if (processTileFace) {
2908  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2909  for (ijk[2] = bbox.min()[2]; ijk[2] <= bbox.max()[2]; ++ijk[2]) {
2910 
2911  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2912  auxData.setYEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2913  }
2914  }
2915  }
2916  }
2917 
2918 
2919  // Eval z-edges
2920 
2921  ijk = bbox.max();
2922  nijk = ijk;
2923  ++nijk[2];
2924 
2925  processTileFace = true;
2926  if (thisDepth >= distAcc.getValueDepth(nijk)) {
2927  processTileFace = thisInside != (distAcc.getValue(nijk) < isoValue);
2928  }
2929 
2930  if (processTileFace) {
2931  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2932  nijk[0] = ijk[0];
2933  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2934  nijk[1] = ijk[1];
2935 
2936  if ((distAcc.getValue(nijk) < isoValue) != thisInside) {
2937  auxData.setZEdge(edgeAcc.getValue(ijk) | char(thisInside), ijk);
2938  }
2939  }
2940  }
2941  }
2942 
2943  ijk = bbox.min();
2944  --ijk[2];
2945 
2946  processTileFace = true;
2947  if (thisDepth >= distAcc.getValueDepth(ijk)) {
2948  processTileFace = !distAcc.probeValue(ijk, value) && thisInside != (value < isoValue);
2949  }
2950 
2951  if (processTileFace) {
2952  for (ijk[0] = bbox.min()[0]; ijk[0] <= bbox.max()[0]; ++ijk[0]) {
2953  for (ijk[1] = bbox.min()[1]; ijk[1] <= bbox.max()[1]; ++ijk[1]) {
2954 
2955  if (!distAcc.probeValue(ijk, value) && (value < isoValue) != thisInside) {
2956  auxData.setZEdge(edgeAcc.getValue(ijk) | char(!thisInside), ijk);
2957  }
2958  }
2959  }
2960  }
2961  }
2962 }
2963 
2964 
2966 
2967 
2968 // Utility class for the volumeToMesh wrapper
2970 {
2971 public:
2972  PointListCopy(const PointList& pointsIn, std::vector<Vec3s>& pointsOut)
2973  : mPointsIn(pointsIn) , mPointsOut(pointsOut)
2974  {
2975  }
2976 
2977  void operator()(const tbb::blocked_range<size_t>& range) const
2978  {
2979  for (size_t n = range.begin(); n < range.end(); ++n) {
2980  mPointsOut[n] = mPointsIn[n];
2981  }
2982  }
2983 
2984 private:
2985  const PointList& mPointsIn;
2986  std::vector<Vec3s>& mPointsOut;
2987 };
2988 
2989 
2990 } // namespace internal
2991 
2992 
2994 
2995 
2996 inline VolumeToMesh::VolumeToMesh(double isovalue, double adaptivity)
2997  : mPointListSize(0)
2998  , mPolygonPoolListSize(0)
2999  , mIsovalue(isovalue)
3000  , mPrimAdaptivity(adaptivity)
3001  , mSecAdaptivity(0.0)
3002  , mRefGrid(GridBase::ConstPtr())
3003  , mSurfaceMaskGrid(GridBase::ConstPtr())
3004  , mAdaptivityGrid(GridBase::ConstPtr())
3005  , mAdaptivityMaskTree(TreeBase::ConstPtr())
3006  , mRefEdgeTree(TreeBase::Ptr())
3007  , mRefTopologyMaskTree(TreeBase::Ptr())
3008  , mSeamPointTree(TreeBase::Ptr())
3009  , mSmoothSeams(false)
3010  , mInvertSurfaceMask(false)
3011  , mPartitions(1)
3012  , mActivePart(0)
3013 {
3014 }
3015 
3016 inline PointList&
3018 {
3019  return mPoints;
3020 }
3021 
3022 inline const size_t&
3024 {
3025  return mPointListSize;
3026 }
3027 
3028 
3029 inline PolygonPoolList&
3031 {
3032  return mPolygons;
3033 }
3034 
3035 
3036 inline const PolygonPoolList&
3038 {
3039  return mPolygons;
3040 }
3041 
3042 
3043 inline const size_t&
3045 {
3046  return mPolygonPoolListSize;
3047 }
3048 
3049 
3050 inline void
3051 VolumeToMesh::setRefGrid(const GridBase::ConstPtr& grid, double secAdaptivity, bool smoothSeams)
3052 {
3053  mRefGrid = grid;
3054  mSecAdaptivity = secAdaptivity;
3055  // Clear out old auxiliary data
3056  mRefEdgeTree = TreeBase::Ptr();
3057  mRefTopologyMaskTree = TreeBase::Ptr();
3058  mSeamPointTree = TreeBase::Ptr();
3059  mSmoothSeams = smoothSeams;
3060 }
3061 
3062 inline void
3064 {
3065  mSurfaceMaskGrid = mask;
3066  mInvertSurfaceMask = invertMask;
3067 }
3068 
3069 inline void
3071 {
3072  mAdaptivityGrid = grid;
3073 }
3074 
3075 inline void
3077 {
3078  mAdaptivityMaskTree = tree;
3079 }
3080 
3081 inline void
3082 VolumeToMesh::partition(unsigned partitions, unsigned activePart)
3083 {
3084  mPartitions = std::max(partitions, unsigned(1));
3085  mActivePart = std::min(activePart, mPartitions-1);
3086 }
3087 
3088 
3089 
3090 
3091 template<typename GridT>
3092 inline void
3093 VolumeToMesh::operator()(const GridT& distGrid)
3094 {
3095  typedef typename GridT::TreeType DistTreeT;
3096  typedef typename DistTreeT::ValueType DistValueT;
3097 
3098  typedef typename DistTreeT::template ValueConverter<char>::Type CharTreeT;
3099  typedef typename DistTreeT::template ValueConverter<int>::Type IntTreeT;
3100  typedef typename DistTreeT::template ValueConverter<Vec3s>::Type Vec3sTreeT;
3101  typedef typename DistTreeT::template ValueConverter<bool>::Type BoolTreeT;
3102  typedef Grid<BoolTreeT> BoolGridT;
3103 
3104  typedef typename DistTreeT::template ValueConverter<float>::Type FloatTreeT;
3105  typedef Grid<FloatTreeT> FloatGridT;
3106 
3107  const openvdb::math::Transform& transform = distGrid.transform();
3108  const DistTreeT& distTree = distGrid.tree();
3109  typename CharTreeT::Ptr edgeTree; // edge flags
3110  typename IntTreeT::Ptr auxTree; // auxiliary data
3111  typename BoolTreeT::Ptr pointMask;
3112 
3113  const bool nonLevelSetGrid = distGrid.getGridClass() != GRID_LEVEL_SET;
3114  const bool noAdaptivity = mPrimAdaptivity < 1e-6 && mSecAdaptivity < 1e-6;
3115 
3116  const bool extraSignCheck = nonLevelSetGrid ||
3117  (std::abs(mIsovalue - double(distGrid.background())) < (1.5 * transform.voxelSize()[0]));
3118 
3119  const BoolGridT * surfaceMask = NULL;
3120  if (mSurfaceMaskGrid) {
3121  surfaceMask = static_cast<const BoolGridT*>(mSurfaceMaskGrid.get());
3122  }
3123 
3124  const FloatGridT * adaptivityField = NULL;
3125  if (mAdaptivityGrid) {
3126  adaptivityField = static_cast<const FloatGridT*>(mAdaptivityGrid.get());
3127  }
3128 
3129  // Collect auxiliary data
3130  {
3131  if (surfaceMask || mPartitions > 1) {
3132 
3133  BoolTreeT valueMask(false);
3134 
3135  if (surfaceMask) {
3136 
3137  // Mask
3138  {
3139  tree::LeafManager<const DistTreeT> leafs(distTree);
3141  *surfaceMask, leafs, transform, mInvertSurfaceMask);
3142  masking.run();
3143  valueMask.merge(masking.tree());
3144  }
3145 
3146  // Partition
3147  if (mPartitions > 1) {
3148  tree::LeafManager<BoolTreeT> leafs(valueMask);
3149  leafs.foreach(internal::PartOp(leafs.leafCount() , mPartitions, mActivePart));
3150  valueMask.pruneInactive();
3151  }
3152 
3153  } else {
3154 
3155  // Partition
3156  tree::LeafManager<const DistTreeT> leafs(distTree);
3157  internal::PartGen<DistTreeT> partitioner(leafs, mPartitions, mActivePart);
3158  partitioner.run();
3159  valueMask.merge(partitioner.tree());
3160  }
3161 
3162  {
3163  tree::LeafManager<BoolTreeT> leafs(valueMask);
3165  op(distTree, leafs, mIsovalue, extraSignCheck);
3166  op.run();
3167 
3168  edgeTree = op.edgeTree();
3169  auxTree = op.auxTree();
3170  }
3171 
3172 
3173  if (!noAdaptivity) {
3174  pointMask = typename BoolTreeT::Ptr(new BoolTreeT(false));
3175  pointMask->topologyUnion(*auxTree);
3176 
3177 // tree::LeafManager<const IntTreeT> auxLeafs(*auxTree);
3178 // internal::BoundaryMaskGen<DistTreeT> boundary(auxLeafs, distTree, valueMask);
3179 // boundary.run();
3180 
3181  {
3182 // tree::LeafManager<BoolTreeT> leafs(boundary.tree());
3183 // internal::AuxiliaryData<DistTreeT, tree::LeafManager<BoolTreeT> >
3184 // op(distTree, leafs, mIsovalue, extraSignCheck);
3185 
3186  tree::LeafManager<const DistTreeT> leafs(distTree);
3188  op(distTree, leafs, mIsovalue, extraSignCheck);
3189 
3190  tree::LeafManager<const IntTreeT> auxLeafs(*auxTree);
3191 
3192  op.run();
3193  auxLeafs.foreach(internal::BoundaryMergeOp<IntTreeT>(*auxTree, *op.auxTree()));
3194  }
3195  }
3196 
3197  } else {
3198 
3199  tree::LeafManager<const DistTreeT> leafs(distTree);
3200  internal::AuxiliaryData<DistTreeT> op(distTree, leafs, mIsovalue, extraSignCheck);
3201  op.run();
3202 
3203  edgeTree = op.edgeTree();
3204  auxTree = op.auxTree();
3205  }
3206 
3207  // Collect auxiliary data from active tiles
3208  if (nonLevelSetGrid) {
3209  internal::tileAuxiliaryData(distTree, *edgeTree, *auxTree, mIsovalue);
3210  }
3211  }
3212 
3213 
3214  // Optionally collect auxiliary data from a reference surface.
3216  if (mRefGrid) {
3217  const GridT* refGrid = static_cast<const GridT*>(mRefGrid.get());
3218  if (refGrid && refGrid->activeVoxelCount() > 0) {
3219 
3220  refData.mDistTree = &refGrid->tree();
3221  refData.mInternalAdaptivity = DistValueT(mSecAdaptivity);
3222 
3223  // Cache reference data for reuse.
3224  if (!mRefEdgeTree && !mRefTopologyMaskTree) {
3227  *refData.mDistTree, leafs, mIsovalue, extraSignCheck);
3228  op.run();
3229  mRefEdgeTree = op.edgeTree();
3230  mRefTopologyMaskTree = op.auxTree();
3231  mSeamPointTree = typename Vec3sTreeT::Ptr(new Vec3sTreeT(Vec3s(0.0)));
3232  }
3233 
3234  if (mRefEdgeTree && mRefTopologyMaskTree) {
3235  refData.mEdgeTree = static_cast<CharTreeT*>(mRefEdgeTree.get());
3236  refData.mTopologyMaskTree = static_cast<BoolTreeT*>(mRefTopologyMaskTree.get());
3237  refData.mSeamPointTree = static_cast<Vec3sTreeT*>(mSeamPointTree.get());
3238  refData.mSeamMaskTree = typename BoolTreeT::Ptr(new BoolTreeT(false));
3239  refData.mSmoothingMaskTree = typename BoolTreeT::Ptr(new BoolTreeT(false));
3240  }
3241  }
3242  }
3243 
3244  BoolTreeT edgeMaskTree(0.0);
3245 
3246  // Generate the seamline mask
3247  if (refData.mSeamMaskTree) {
3248  refData.mSeamMaskTree->topologyUnion(*auxTree.get());
3249 
3250  internal::LeafPtrList<BoolTreeT> leafs(*refData.mSeamMaskTree.get());
3251  internal::SeamMaskGen<DistTreeT> op(leafs, *refData.mTopologyMaskTree, *auxTree.get());
3252  op.runParallel();
3253 
3254  refData.mSeamMaskTree->pruneInactive();
3255 
3256  edgeMaskTree.topologyUnion(*refData.mSeamMaskTree);
3257 
3258  dilateVoxels(*refData.mSeamMaskTree);
3259  dilateVoxels(*refData.mSeamMaskTree);
3260  dilateVoxels(*refData.mSeamMaskTree);
3261 
3262 
3263  if (mAdaptivityMaskTree) {
3264 
3265  const BoolTreeT* refGrid = static_cast<const BoolTreeT*>(
3266  mAdaptivityMaskTree.get());
3267 
3268  refData.mSeamMaskTree->topologyUnion(*refGrid);
3269  }
3270  }
3271 
3272  // Process auxiliary data
3273  {
3274  internal::LeafPtrList<IntTreeT> auxLeafs(*auxTree);
3275  std::vector<size_t> regions(auxLeafs.size(), 0);
3276 
3277  {
3278  if (noAdaptivity) {
3279  internal::Count<IntTreeT> count(auxLeafs, regions);
3280  count.runParallel();
3281  } else {
3282  internal::Merge<DistTreeT> merge(distTree, auxLeafs,
3283  regions, DistValueT(mIsovalue), DistValueT(mPrimAdaptivity), pointMask.get());
3284 
3285  merge.setRefData(refData);
3286 
3287  if (adaptivityField) {
3288  merge.setSpatialAdaptivity(transform, *adaptivityField);
3289  }
3290 
3291  merge.run();
3292  }
3293 
3294  mPointListSize = 0;
3295  size_t tmp = 0;
3296  for (size_t n = 0, N = regions.size(); n < N; ++n) {
3297  tmp = regions[n];
3298  regions[n] = mPointListSize;
3299  mPointListSize += tmp;
3300  }
3301  }
3302 
3303  if (refData.isValid()) { // match leaf topology
3305  for (size_t n = 0, N = auxLeafs.size(); n < N; ++n) {
3306  acc.touchLeaf(auxLeafs[n]->getOrigin());
3307  }
3308  }
3309 
3310  // Generate the unique point list
3311  mPoints.reset(new openvdb::Vec3s[mPointListSize]);
3312 
3314  pointGen(distTree, auxLeafs, regions, transform, mPoints, mIsovalue);
3315  pointGen.setRefData(refData);
3316  pointGen.runParallel();
3317  }
3318 
3319  // Smooth seam line edges
3320  if (mSmoothSeams && refData.isValid()) {
3321 
3322  refData.mSmoothingMaskTree->pruneInactive();
3324  internal::EdgeSmooth<DistTreeT> op(leafs, edgeMaskTree, *auxTree, mPoints, transform);
3325  op.runParallel(3);
3326 
3327  // Cache shared points
3329  tree::ValueAccessor<IntTreeT> auxAcc(*auxTree);
3330  Coord ijk;
3331 
3332  typename BoolTreeT::LeafIter leafIt = refData.mSeamMaskTree->beginLeaf();
3333  for ( ; leafIt; ++leafIt) {
3334  typename BoolTreeT::LeafNodeType::ValueOnIter valueIt = leafIt->beginValueOn();
3335  for ( ; valueIt; ++valueIt) {
3336  ijk = valueIt.getCoord();
3337  const int idx = auxAcc.getValue(ijk);
3338  if (idx != 0 && !ptnAcc.isValueOn(ijk)) {
3339  ptnAcc.setValue(ijk, mPoints[idx]);
3340  }
3341  }
3342  }
3343  }
3344 
3345  // Generate mesh
3346  {
3347  internal::LeafCPtrList<CharTreeT> edgeLeafs(*edgeTree);
3348  mPolygonPoolListSize = edgeLeafs.size();
3349  mPolygons.reset(new PolygonPool[mPolygonPoolListSize]);
3350 
3351  if (noAdaptivity) {
3353  meshGen(edgeLeafs, *auxTree, mPolygons);
3354  meshGen.setRefData(refData);
3355  meshGen.runParallel();
3356  } else {
3358  meshGen(edgeLeafs, *auxTree, mPolygons);
3359  meshGen.setRefData(refData);
3360  meshGen.runParallel();
3361  }
3362  }
3363 }
3364 
3365 
3367 
3368 
3369 template<typename GridType>
3370 inline void
3372  const GridType& grid,
3373  std::vector<Vec3s>& points,
3374  std::vector<Vec3I>& triangles,
3375  std::vector<Vec4I>& quads,
3376  double isovalue,
3377  double adaptivity)
3378 {
3379  VolumeToMesh mesher(isovalue, adaptivity);
3380  mesher(grid);
3381 
3382  // Preallocate the point list
3383  points.clear();
3384  points.resize(mesher.pointListSize());
3385 
3386  { // Copy points
3387  internal::PointListCopy ptnCpy(mesher.pointList(), points);
3388  tbb::parallel_for(tbb::blocked_range<size_t>(0, points.size()), ptnCpy);
3389  mesher.pointList().reset(NULL);
3390  }
3391 
3392  PolygonPoolList& polygonPoolList = mesher.polygonPoolList();
3393 
3394  { // Preallocate primitive lists
3395  size_t numQuads = 0, numTriangles = 0;
3396  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
3397  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
3398  numTriangles += polygons.numTriangles();
3399  numQuads += polygons.numQuads();
3400  }
3401 
3402  triangles.clear();
3403  triangles.resize(numTriangles);
3404  quads.clear();
3405  quads.resize(numQuads);
3406  }
3407 
3408  // Copy primitives
3409  size_t qIdx = 0, tIdx = 0;
3410  for (size_t n = 0, N = mesher.polygonPoolListSize(); n < N; ++n) {
3411  openvdb::tools::PolygonPool& polygons = polygonPoolList[n];
3412 
3413  for (size_t i = 0, I = polygons.numQuads(); i < I; ++i) {
3414  quads[qIdx++] = polygons.quad(i);
3415  }
3416 
3417  for (size_t i = 0, I = polygons.numTriangles(); i < I; ++i) {
3418  triangles[tIdx++] = polygons.triangle(i);
3419  }
3420  }
3421 }
3422 
3423 
3424 template<typename GridType>
3425 void
3427  const GridType& grid,
3428  std::vector<Vec3s>& points,
3429  std::vector<Vec4I>& quads,
3430  double isovalue)
3431 {
3432  std::vector<Vec3I> triangles(0);
3433  volumeToMesh(grid,points, triangles, quads, isovalue, 0.0);
3434 }
3435 
3436 } // namespace tools
3437 } // namespace OPENVDB_VERSION_NAME
3438 } // namespace openvdb
3439 
3440 #endif // OPENVDB_TOOLS_VOLUME_TO_MESH_HAS_BEEN_INCLUDED
3441 
3442 // Copyright (c) 2012-2013 DreamWorks Animation LLC
3443 // All rights reserved. This software is distributed under the
3444 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )