OpenVDB 11.0.0
Loading...
Searching...
No Matches
MeshToVolume.h
Go to the documentation of this file.
1// Copyright Contributors to the OpenVDB Project
2// SPDX-License-Identifier: MPL-2.0
3
4/// @file MeshToVolume.h
5///
6/// @brief Convert polygonal meshes that consist of quads and/or triangles
7/// into signed or unsigned distance field volumes.
8///
9/// @note The signed distance field conversion requires a closed surface
10/// but not necessarily a manifold surface. Supports surfaces with
11/// self intersections and degenerate faces and is independent of
12/// mesh surface normals / polygon orientation.
13///
14/// @author Mihai Alden
15
16#ifndef OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
17#define OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
18
19#include <openvdb/Platform.h> // for OPENVDB_HAS_CXX11
20#include <openvdb/Types.h>
21#include <openvdb/math/FiniteDifference.h> // for GodunovsNormSqrd
22#include <openvdb/math/Proximity.h> // for closestPointOnTriangleToPoint
24#include <openvdb/util/Util.h>
25#include <openvdb/thread/Threading.h>
26#include <openvdb/openvdb.h>
27
28#include "ChangeBackground.h"
29#include "Prune.h" // for pruneInactive and pruneLevelSet
30#include "SignedFloodFill.h" // for signedFloodFillWithValues
31
32#include <tbb/blocked_range.h>
33#include <tbb/enumerable_thread_specific.h>
34#include <tbb/parallel_for.h>
35#include <tbb/parallel_reduce.h>
36#include <tbb/partitioner.h>
37#include <tbb/task_group.h>
38#include <tbb/task_arena.h>
39
40#include <algorithm> // for std::sort()
41#include <cmath> // for std::isfinite(), std::isnan()
42#include <deque>
43#include <limits>
44#include <memory>
45#include <sstream>
46#include <type_traits>
47#include <vector>
48
49namespace openvdb {
51namespace OPENVDB_VERSION_NAME {
52namespace tools {
53
54
55////////////////////////////////////////
56
57
58/// @brief Mesh to volume conversion flags
60
61 /// Switch from the default signed distance field conversion that classifies
62 /// regions as either inside or outside the mesh boundary to a unsigned distance
63 /// field conversion that only computes distance values. This conversion type
64 /// does not require a closed watertight mesh.
66
67 /// Disable the cleanup step that removes voxels created by self intersecting
68 /// portions of the mesh.
70
71 /// Disable the distance renormalization step that smooths out bumps caused
72 /// by self intersecting or overlapping portions of the mesh
74
75 /// Disable the cleanup step that removes active voxels that exceed the
76 /// narrow band limits. (Only relevant for small limits)
78};
79
80
81/// @brief Different staregies how to determine sign of an SDF when using
82/// interior test.
84
85 /// Evaluates interior test at every voxel. This is usefull when we rebuild already
86 /// existing SDF where evaluating previous grid is cheap
88
89 /// Evaluates interior test at least once per tile and flood fills within the tile.
91};
92
93
94/// @brief Convert polygonal meshes that consist of quads and/or triangles into
95/// signed or unsigned distance field volumes.
96///
97/// @note Requires a closed surface but not necessarily a manifold surface.
98/// Supports surfaces with self intersections and degenerate faces
99/// and is independent of mesh surface normals.
100///
101/// @interface MeshDataAdapter
102/// Expected interface for the MeshDataAdapter class
103/// @code
104/// struct MeshDataAdapter {
105/// size_t polygonCount() const; // Total number of polygons
106/// size_t pointCount() const; // Total number of points
107/// size_t vertexCount(size_t n) const; // Vertex count for polygon n
108///
109/// // Return position pos in local grid index space for polygon n and vertex v
110/// void getIndexSpacePoint(size_t n, size_t v, openvdb::Vec3d& pos) const;
111/// };
112/// @endcode
113///
114/// @param mesh mesh data access class that conforms to the MeshDataAdapter
115/// interface
116/// @param transform world-to-index-space transform
117/// @param exteriorBandWidth exterior narrow band width in voxel units
118/// @param interiorBandWidth interior narrow band width in voxel units
119/// (set to std::numeric_limits<float>::max() to fill object
120/// interior with distance values)
121/// @param flags optional conversion flags defined in @c MeshToVolumeFlags
122/// @param polygonIndexGrid optional grid output that will contain the closest-polygon
123/// index for each voxel in the narrow band region
124/// @param interiorTest function `Coord -> Bool` that evaluates to true inside of the
125/// mesh and false outside, for more see evaluateInteriortest
126/// @param interiorTestStrat determines how the interiorTest is used, see InteriorTestStrategy
127template <typename GridType, typename MeshDataAdapter, typename InteriorTest = std::nullptr_t>
128typename GridType::Ptr
129meshToVolume(
130 const MeshDataAdapter& mesh,
131 const math::Transform& transform,
132 float exteriorBandWidth = 3.0f,
133 float interiorBandWidth = 3.0f,
134 int flags = 0,
135 typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = nullptr,
136 InteriorTest interiorTest = nullptr,
137 InteriorTestStrategy interiorTestStrat = EVAL_EVERY_VOXEL);
138
139
140/// @brief Convert polygonal meshes that consist of quads and/or triangles into
141/// signed or unsigned distance field volumes.
142///
143/// @param interrupter a callback to interrupt the conversion process that conforms
144/// to the util::NullInterrupter interface
145/// @param mesh mesh data access class that conforms to the MeshDataAdapter
146/// interface
147/// @param transform world-to-index-space transform
148/// @param exteriorBandWidth exterior narrow band width in voxel units
149/// @param interiorBandWidth interior narrow band width in voxel units (set this value to
150/// std::numeric_limits<float>::max() to fill interior regions
151/// with distance values)
152/// @param flags optional conversion flags defined in @c MeshToVolumeFlags
153/// @param polygonIndexGrid optional grid output that will contain the closest-polygon
154/// index for each voxel in the active narrow band region
155/// @param interiorTest function `Coord -> Bool` that evaluates to true inside of the
156/// mesh and false outside, for more see evaluatInteriorTest
157/// @param interiorTestStrat determines how the interiorTest is used, see InteriorTestStrategy
158template <typename GridType, typename MeshDataAdapter, typename Interrupter, typename InteriorTest = std::nullptr_t>
159typename GridType::Ptr
160meshToVolume(
161 Interrupter& interrupter,
162 const MeshDataAdapter& mesh,
163 const math::Transform& transform,
164 float exteriorBandWidth = 3.0f,
165 float interiorBandWidth = 3.0f,
166 int flags = 0,
167 typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid = nullptr,
168 InteriorTest interiorTest = nullptr,
169 InteriorTestStrategy interiorTestStrategy = EVAL_EVERY_VOXEL);
170
171
172////////////////////////////////////////
173
174
175/// @brief Contiguous quad and triangle data adapter class
176///
177/// @details PointType and PolygonType must provide element access
178/// through the square brackets operator.
179/// @details Points are assumed to be in local grid index space.
180/// @details The PolygonType tuple can have either three or four components
181/// this property must be specified in a static member variable
182/// named @c size, similar to the math::Tuple class.
183/// @details A four component tuple can represent a quads or a triangle
184/// if the fourth component set to @c util::INVALID_INDEX
185template<typename PointType, typename PolygonType>
187
188 QuadAndTriangleDataAdapter(const std::vector<PointType>& points,
189 const std::vector<PolygonType>& polygons)
190 : mPointArray(points.empty() ? nullptr : &points[0])
191 , mPointArraySize(points.size())
192 , mPolygonArray(polygons.empty() ? nullptr : &polygons[0])
193 , mPolygonArraySize(polygons.size())
194 {
195 }
196
197 QuadAndTriangleDataAdapter(const PointType * pointArray, size_t pointArraySize,
198 const PolygonType* polygonArray, size_t polygonArraySize)
199 : mPointArray(pointArray)
200 , mPointArraySize(pointArraySize)
201 , mPolygonArray(polygonArray)
202 , mPolygonArraySize(polygonArraySize)
203 {
204 }
205
206 size_t polygonCount() const { return mPolygonArraySize; }
207 size_t pointCount() const { return mPointArraySize; }
208
209 /// @brief Vertex count for polygon @a n
210 size_t vertexCount(size_t n) const {
211 return (PolygonType::size == 3 || mPolygonArray[n][3] == util::INVALID_IDX) ? 3 : 4;
212 }
213
214 /// @brief Returns position @a pos in local grid index space
215 /// for polygon @a n and vertex @a v
216 void getIndexSpacePoint(size_t n, size_t v, Vec3d& pos) const {
217 const PointType& p = mPointArray[mPolygonArray[n][int(v)]];
218 pos[0] = double(p[0]);
219 pos[1] = double(p[1]);
220 pos[2] = double(p[2]);
221 }
222
223private:
224 PointType const * const mPointArray;
225 size_t const mPointArraySize;
226 PolygonType const * const mPolygonArray;
227 size_t const mPolygonArraySize;
228}; // struct QuadAndTriangleDataAdapter
229
230
231////////////////////////////////////////
232
233
234// Convenience functions for the mesh to volume converter that wrap stl containers.
235//
236// Note the meshToVolume() method declared above is more flexible and better suited
237// for arbitrary data structures.
238
239
240/// @brief Convert a triangle mesh to a level set volume.
241///
242/// @return a grid of type @c GridType containing a narrow-band level set
243/// representation of the input mesh.
244///
245/// @throw TypeError if @c GridType is not scalar or not floating-point
246///
247/// @note Requires a closed surface but not necessarily a manifold surface.
248/// Supports surfaces with self intersections and degenerate faces
249/// and is independent of mesh surface normals.
250///
251/// @param xform transform for the output grid
252/// @param points list of world space point positions
253/// @param triangles triangle index list
254/// @param halfWidth half the width of the narrow band, in voxel units
255template<typename GridType>
256typename GridType::Ptr
257meshToLevelSet(
258 const openvdb::math::Transform& xform,
259 const std::vector<Vec3s>& points,
260 const std::vector<Vec3I>& triangles,
261 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
262
263/// Adds support for a @a interrupter callback used to cancel the conversion.
264template<typename GridType, typename Interrupter>
265typename GridType::Ptr
266meshToLevelSet(
267 Interrupter& interrupter,
268 const openvdb::math::Transform& xform,
269 const std::vector<Vec3s>& points,
270 const std::vector<Vec3I>& triangles,
271 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
272
273
274/// @brief Convert a quad mesh to a level set volume.
275///
276/// @return a grid of type @c GridType containing a narrow-band level set
277/// representation of the input mesh.
278///
279/// @throw TypeError if @c GridType is not scalar or not floating-point
280///
281/// @note Requires a closed surface but not necessarily a manifold surface.
282/// Supports surfaces with self intersections and degenerate faces
283/// and is independent of mesh surface normals.
284///
285/// @param xform transform for the output grid
286/// @param points list of world space point positions
287/// @param quads quad index list
288/// @param halfWidth half the width of the narrow band, in voxel units
289template<typename GridType>
290typename GridType::Ptr
291meshToLevelSet(
292 const openvdb::math::Transform& xform,
293 const std::vector<Vec3s>& points,
294 const std::vector<Vec4I>& quads,
295 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
296
297/// Adds support for a @a interrupter callback used to cancel the conversion.
298template<typename GridType, typename Interrupter>
299typename GridType::Ptr
300meshToLevelSet(
301 Interrupter& interrupter,
302 const openvdb::math::Transform& xform,
303 const std::vector<Vec3s>& points,
304 const std::vector<Vec4I>& quads,
305 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
306
307
308/// @brief Convert a triangle and quad mesh to a level set volume.
309///
310/// @return a grid of type @c GridType containing a narrow-band level set
311/// representation of the input mesh.
312///
313/// @throw TypeError if @c GridType is not scalar or not floating-point
314///
315/// @note Requires a closed surface but not necessarily a manifold surface.
316/// Supports surfaces with self intersections and degenerate faces
317/// and is independent of mesh surface normals.
318///
319/// @param xform transform for the output grid
320/// @param points list of world space point positions
321/// @param triangles triangle index list
322/// @param quads quad index list
323/// @param halfWidth half the width of the narrow band, in voxel units
324template<typename GridType>
325typename GridType::Ptr
326meshToLevelSet(
327 const openvdb::math::Transform& xform,
328 const std::vector<Vec3s>& points,
329 const std::vector<Vec3I>& triangles,
330 const std::vector<Vec4I>& quads,
331 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
332
333/// Adds support for a @a interrupter callback used to cancel the conversion.
334template<typename GridType, typename Interrupter>
335typename GridType::Ptr
336meshToLevelSet(
337 Interrupter& interrupter,
338 const openvdb::math::Transform& xform,
339 const std::vector<Vec3s>& points,
340 const std::vector<Vec3I>& triangles,
341 const std::vector<Vec4I>& quads,
342 float halfWidth = float(LEVEL_SET_HALF_WIDTH));
343
344
345/// @brief Convert a triangle and quad mesh to a signed distance field
346/// with an asymmetrical narrow band.
347///
348/// @return a grid of type @c GridType containing a narrow-band signed
349/// distance field representation of the input mesh.
350///
351/// @throw TypeError if @c GridType is not scalar or not floating-point
352///
353/// @note Requires a closed surface but not necessarily a manifold surface.
354/// Supports surfaces with self intersections and degenerate faces
355/// and is independent of mesh surface normals.
356///
357/// @param xform transform for the output grid
358/// @param points list of world space point positions
359/// @param triangles triangle index list
360/// @param quads quad index list
361/// @param exBandWidth the exterior narrow-band width in voxel units
362/// @param inBandWidth the interior narrow-band width in voxel units
363template<typename GridType>
364typename GridType::Ptr
365meshToSignedDistanceField(
366 const openvdb::math::Transform& xform,
367 const std::vector<Vec3s>& points,
368 const std::vector<Vec3I>& triangles,
369 const std::vector<Vec4I>& quads,
370 float exBandWidth,
371 float inBandWidth);
372
373/// Adds support for a @a interrupter callback used to cancel the conversion.
374template<typename GridType, typename Interrupter>
375typename GridType::Ptr
376meshToSignedDistanceField(
377 Interrupter& interrupter,
378 const openvdb::math::Transform& xform,
379 const std::vector<Vec3s>& points,
380 const std::vector<Vec3I>& triangles,
381 const std::vector<Vec4I>& quads,
382 float exBandWidth,
383 float inBandWidth);
384
385
386/// @brief Convert a triangle and quad mesh to an unsigned distance field.
387///
388/// @return a grid of type @c GridType containing a narrow-band unsigned
389/// distance field representation of the input mesh.
390///
391/// @throw TypeError if @c GridType is not scalar or not floating-point
392///
393/// @note Does not requires a closed surface.
394///
395/// @param xform transform for the output grid
396/// @param points list of world space point positions
397/// @param triangles triangle index list
398/// @param quads quad index list
399/// @param bandWidth the width of the narrow band, in voxel units
400template<typename GridType>
401typename GridType::Ptr
402meshToUnsignedDistanceField(
403 const openvdb::math::Transform& xform,
404 const std::vector<Vec3s>& points,
405 const std::vector<Vec3I>& triangles,
406 const std::vector<Vec4I>& quads,
407 float bandWidth);
408
409/// Adds support for a @a interrupter callback used to cancel the conversion.
410template<typename GridType, typename Interrupter>
411typename GridType::Ptr
412meshToUnsignedDistanceField(
413 Interrupter& interrupter,
414 const openvdb::math::Transform& xform,
415 const std::vector<Vec3s>& points,
416 const std::vector<Vec3I>& triangles,
417 const std::vector<Vec4I>& quads,
418 float bandWidth);
419
420
421////////////////////////////////////////
422
423
424/// @brief Return a grid of type @c GridType containing a narrow-band level set
425/// representation of a box.
426///
427/// @param bbox a bounding box in world units
428/// @param xform world-to-index-space transform
429/// @param halfWidth half the width of the narrow band, in voxel units
430template<typename GridType, typename VecType>
431typename GridType::Ptr
432createLevelSetBox(const math::BBox<VecType>& bbox,
433 const openvdb::math::Transform& xform,
434 typename VecType::ValueType halfWidth = LEVEL_SET_HALF_WIDTH);
435
436
437////////////////////////////////////////
438
439
440/// @brief Traces the exterior voxel boundary of closed objects in the input
441/// volume @a tree. Exterior voxels are marked with a negative sign,
442/// voxels with a value below @c 0.75 are left unchanged and act as
443/// the boundary layer.
444///
445/// @note Does not propagate sign information into tile regions.
446template <typename FloatTreeT>
447void
448traceExteriorBoundaries(FloatTreeT& tree);
449
450
451////////////////////////////////////////
452
453
454/// @brief Extracts and stores voxel edge intersection data from a mesh.
456{
457public:
458
459 //////////
460
461 ///@brief Internal edge data type.
462 struct EdgeData {
463 EdgeData(float dist = 1.0)
464 : mXDist(dist), mYDist(dist), mZDist(dist)
465 , mXPrim(util::INVALID_IDX)
466 , mYPrim(util::INVALID_IDX)
467 , mZPrim(util::INVALID_IDX)
468 {
469 }
470
471 //@{
472 /// Required by several of the tree nodes
473 /// @note These methods don't perform meaningful operations.
474 bool operator< (const EdgeData&) const { return false; }
475 bool operator> (const EdgeData&) const { return false; }
476 template<class T> EdgeData operator+(const T&) const { return *this; }
477 template<class T> EdgeData operator-(const T&) const { return *this; }
478 EdgeData operator-() const { return *this; }
479 //@}
480
481 bool operator==(const EdgeData& rhs) const
482 {
483 return mXPrim == rhs.mXPrim && mYPrim == rhs.mYPrim && mZPrim == rhs.mZPrim;
484 }
485
486 float mXDist, mYDist, mZDist;
487 Index32 mXPrim, mYPrim, mZPrim;
488 };
489
492
493
494 //////////
495
496
498
499
500 /// @brief Threaded method to extract voxel edge data, the closest
501 /// intersection point and corresponding primitive index,
502 /// from the given mesh.
503 ///
504 /// @param pointList List of points in grid index space, preferably unique
505 /// and shared by different polygons.
506 /// @param polygonList List of triangles and/or quads.
507 void convert(const std::vector<Vec3s>& pointList, const std::vector<Vec4I>& polygonList);
508
509
510 /// @brief Returns intersection points with corresponding primitive
511 /// indices for the given @c ijk voxel.
512 void getEdgeData(Accessor& acc, const Coord& ijk,
513 std::vector<Vec3d>& points, std::vector<Index32>& primitives);
514
515 /// @return An accessor of @c MeshToVoxelEdgeData::Accessor type that
516 /// provides random read access to the internal tree.
517 Accessor getAccessor() { return Accessor(mTree); }
518
519private:
520 void operator=(const MeshToVoxelEdgeData&) {}
521 TreeType mTree;
522 class GenEdgeData;
523};
524
525
526////////////////////////////////////////////////////////////////////////////////
527////////////////////////////////////////////////////////////////////////////////
528
529
530// Internal utility objects and implementation details
531
532/// @cond OPENVDB_DOCS_INTERNAL
533
534namespace mesh_to_volume_internal {
535
536template<typename PointType>
537struct TransformPoints {
538
539 TransformPoints(const PointType* pointsIn, PointType* pointsOut,
540 const math::Transform& xform)
541 : mPointsIn(pointsIn), mPointsOut(pointsOut), mXform(&xform)
542 {
543 }
544
545 void operator()(const tbb::blocked_range<size_t>& range) const {
546
547 Vec3d pos;
548
549 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
550
551 const PointType& wsP = mPointsIn[n];
552 pos[0] = double(wsP[0]);
553 pos[1] = double(wsP[1]);
554 pos[2] = double(wsP[2]);
555
556 pos = mXform->worldToIndex(pos);
557
558 PointType& isP = mPointsOut[n];
559 isP[0] = typename PointType::value_type(pos[0]);
560 isP[1] = typename PointType::value_type(pos[1]);
561 isP[2] = typename PointType::value_type(pos[2]);
562 }
563 }
564
565 PointType const * const mPointsIn;
566 PointType * const mPointsOut;
567 math::Transform const * const mXform;
568}; // TransformPoints
569
570
571template<typename ValueType>
572struct Tolerance
573{
574 static ValueType epsilon() { return ValueType(1e-7); }
575 static ValueType minNarrowBandWidth() { return ValueType(1.0 + 1e-6); }
576};
577
578
579////////////////////////////////////////
580
581
582template<typename TreeType>
583class CombineLeafNodes
584{
585public:
586
587 using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
588
589 using LeafNodeType = typename TreeType::LeafNodeType;
590 using Int32LeafNodeType = typename Int32TreeType::LeafNodeType;
591
592 CombineLeafNodes(TreeType& lhsDistTree, Int32TreeType& lhsIdxTree,
593 LeafNodeType ** rhsDistNodes, Int32LeafNodeType ** rhsIdxNodes)
594 : mDistTree(&lhsDistTree)
595 , mIdxTree(&lhsIdxTree)
596 , mRhsDistNodes(rhsDistNodes)
597 , mRhsIdxNodes(rhsIdxNodes)
598 {
599 }
600
601 void operator()(const tbb::blocked_range<size_t>& range) const {
602
603 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
604 tree::ValueAccessor<Int32TreeType> idxAcc(*mIdxTree);
605
606 using DistValueType = typename LeafNodeType::ValueType;
607 using IndexValueType = typename Int32LeafNodeType::ValueType;
608
609 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
610
611 const Coord& origin = mRhsDistNodes[n]->origin();
612
613 LeafNodeType* lhsDistNode = distAcc.probeLeaf(origin);
614 Int32LeafNodeType* lhsIdxNode = idxAcc.probeLeaf(origin);
615
616 DistValueType* lhsDistData = lhsDistNode->buffer().data();
617 IndexValueType* lhsIdxData = lhsIdxNode->buffer().data();
618
619 const DistValueType* rhsDistData = mRhsDistNodes[n]->buffer().data();
620 const IndexValueType* rhsIdxData = mRhsIdxNodes[n]->buffer().data();
621
622
623 for (Index32 offset = 0; offset < LeafNodeType::SIZE; ++offset) {
624
625 if (rhsIdxData[offset] != Int32(util::INVALID_IDX)) {
626
627 const DistValueType& lhsValue = lhsDistData[offset];
628 const DistValueType& rhsValue = rhsDistData[offset];
629
630 if (rhsValue < lhsValue) {
631 lhsDistNode->setValueOn(offset, rhsValue);
632 lhsIdxNode->setValueOn(offset, rhsIdxData[offset]);
633 } else if (math::isExactlyEqual(rhsValue, lhsValue)) {
634 lhsIdxNode->setValueOn(offset,
635 std::min(lhsIdxData[offset], rhsIdxData[offset]));
636 }
637 }
638 }
639
640 delete mRhsDistNodes[n];
641 delete mRhsIdxNodes[n];
642 }
643 }
644
645private:
646
647 TreeType * const mDistTree;
648 Int32TreeType * const mIdxTree;
649
650 LeafNodeType ** const mRhsDistNodes;
651 Int32LeafNodeType ** const mRhsIdxNodes;
652}; // class CombineLeafNodes
653
654
655////////////////////////////////////////
656
657
658template<typename TreeType>
659struct StashOriginAndStoreOffset
660{
661 using LeafNodeType = typename TreeType::LeafNodeType;
662
663 StashOriginAndStoreOffset(std::vector<LeafNodeType*>& nodes, Coord* coordinates)
664 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
665 {
666 }
667
668 void operator()(const tbb::blocked_range<size_t>& range) const {
669 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
670 Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
671 mCoordinates[n] = origin;
672 origin[0] = static_cast<int>(n);
673 }
674 }
675
676 LeafNodeType ** const mNodes;
677 Coord * const mCoordinates;
678};
679
680
681template<typename TreeType>
682struct RestoreOrigin
683{
684 using LeafNodeType = typename TreeType::LeafNodeType;
685
686 RestoreOrigin(std::vector<LeafNodeType*>& nodes, const Coord* coordinates)
687 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mCoordinates(coordinates)
688 {
689 }
690
691 void operator()(const tbb::blocked_range<size_t>& range) const {
692 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
693 Coord& origin = const_cast<Coord&>(mNodes[n]->origin());
694 origin[0] = mCoordinates[n][0];
695 }
696 }
697
698 LeafNodeType ** const mNodes;
699 Coord const * const mCoordinates;
700};
701
702
703template<typename TreeType>
704class ComputeNodeConnectivity
705{
706public:
707 using LeafNodeType = typename TreeType::LeafNodeType;
708
709 ComputeNodeConnectivity(const TreeType& tree, const Coord* coordinates,
710 size_t* offsets, size_t numNodes, const CoordBBox& bbox)
711 : mTree(&tree)
712 , mCoordinates(coordinates)
713 , mOffsets(offsets)
714 , mNumNodes(numNodes)
715 , mBBox(bbox)
716 {
717 }
718
719 ComputeNodeConnectivity(const ComputeNodeConnectivity&) = default;
720
721 // Disallow assignment
722 ComputeNodeConnectivity& operator=(const ComputeNodeConnectivity&) = delete;
723
724 void operator()(const tbb::blocked_range<size_t>& range) const {
725
726 size_t* offsetsNextX = mOffsets;
727 size_t* offsetsPrevX = mOffsets + mNumNodes;
728 size_t* offsetsNextY = mOffsets + mNumNodes * 2;
729 size_t* offsetsPrevY = mOffsets + mNumNodes * 3;
730 size_t* offsetsNextZ = mOffsets + mNumNodes * 4;
731 size_t* offsetsPrevZ = mOffsets + mNumNodes * 5;
732
733 tree::ValueAccessor<const TreeType> acc(*mTree);
734 Coord ijk;
735 const Int32 DIM = static_cast<Int32>(LeafNodeType::DIM);
736
737 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
738 const Coord& origin = mCoordinates[n];
739 offsetsNextX[n] = findNeighbourNode(acc, origin, Coord(DIM, 0, 0));
740 offsetsPrevX[n] = findNeighbourNode(acc, origin, Coord(-DIM, 0, 0));
741 offsetsNextY[n] = findNeighbourNode(acc, origin, Coord(0, DIM, 0));
742 offsetsPrevY[n] = findNeighbourNode(acc, origin, Coord(0, -DIM, 0));
743 offsetsNextZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, DIM));
744 offsetsPrevZ[n] = findNeighbourNode(acc, origin, Coord(0, 0, -DIM));
745 }
746 }
747
748 size_t findNeighbourNode(tree::ValueAccessor<const TreeType>& acc,
749 const Coord& start, const Coord& step) const
750 {
751 Coord ijk = start + step;
752 CoordBBox bbox(mBBox);
753
754 while (bbox.isInside(ijk)) {
755 const LeafNodeType* node = acc.probeConstLeaf(ijk);
756 if (node) return static_cast<size_t>(node->origin()[0]);
757 ijk += step;
758 }
759
760 return std::numeric_limits<size_t>::max();
761 }
762
763
764private:
765 TreeType const * const mTree;
766 Coord const * const mCoordinates;
767 size_t * const mOffsets;
768
769 const size_t mNumNodes;
770 const CoordBBox mBBox;
771}; // class ComputeNodeConnectivity
772
773
774template<typename TreeType>
775struct LeafNodeConnectivityTable
776{
777 enum { INVALID_OFFSET = std::numeric_limits<size_t>::max() };
778
779 using LeafNodeType = typename TreeType::LeafNodeType;
780
781 LeafNodeConnectivityTable(TreeType& tree)
782 {
783 mLeafNodes.reserve(tree.leafCount());
784 tree.getNodes(mLeafNodes);
785
786 if (mLeafNodes.empty()) return;
787
788 CoordBBox bbox;
789 tree.evalLeafBoundingBox(bbox);
790
791 const tbb::blocked_range<size_t> range(0, mLeafNodes.size());
792
793 // stash the leafnode origin coordinate and temporarily store the
794 // linear offset in the origin.x variable.
795 std::unique_ptr<Coord[]> coordinates{new Coord[mLeafNodes.size()]};
796 tbb::parallel_for(range,
797 StashOriginAndStoreOffset<TreeType>(mLeafNodes, coordinates.get()));
798
799 // build the leafnode offset table
800 mOffsets.reset(new size_t[mLeafNodes.size() * 6]);
801
802
803 tbb::parallel_for(range, ComputeNodeConnectivity<TreeType>(
804 tree, coordinates.get(), mOffsets.get(), mLeafNodes.size(), bbox));
805
806 // restore the leafnode origin coordinate
807 tbb::parallel_for(range, RestoreOrigin<TreeType>(mLeafNodes, coordinates.get()));
808 }
809
810 size_t size() const { return mLeafNodes.size(); }
811
812 std::vector<LeafNodeType*>& nodes() { return mLeafNodes; }
813 const std::vector<LeafNodeType*>& nodes() const { return mLeafNodes; }
814
815
816 const size_t* offsetsNextX() const { return mOffsets.get(); }
817 const size_t* offsetsPrevX() const { return mOffsets.get() + mLeafNodes.size(); }
818
819 const size_t* offsetsNextY() const { return mOffsets.get() + mLeafNodes.size() * 2; }
820 const size_t* offsetsPrevY() const { return mOffsets.get() + mLeafNodes.size() * 3; }
821
822 const size_t* offsetsNextZ() const { return mOffsets.get() + mLeafNodes.size() * 4; }
823 const size_t* offsetsPrevZ() const { return mOffsets.get() + mLeafNodes.size() * 5; }
824
825private:
826 std::vector<LeafNodeType*> mLeafNodes;
827 std::unique_ptr<size_t[]> mOffsets;
828}; // struct LeafNodeConnectivityTable
829
830
831template<typename TreeType>
832class SweepExteriorSign
833{
834public:
835
836 enum Axis { X_AXIS = 0, Y_AXIS = 1, Z_AXIS = 2 };
837
838 using ValueType = typename TreeType::ValueType;
839 using LeafNodeType = typename TreeType::LeafNodeType;
840 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
841
842 SweepExteriorSign(Axis axis, const std::vector<size_t>& startNodeIndices,
843 ConnectivityTable& connectivity)
844 : mStartNodeIndices(startNodeIndices.empty() ? nullptr : &startNodeIndices[0])
845 , mConnectivity(&connectivity)
846 , mAxis(axis)
847 {
848 }
849
850 void operator()(const tbb::blocked_range<size_t>& range) const {
851
852 constexpr Int32 DIM = static_cast<Int32>(LeafNodeType::DIM);
853
854 std::vector<LeafNodeType*>& nodes = mConnectivity->nodes();
855
856 // Z Axis
857 size_t idxA = 0, idxB = 1;
858 Int32 step = 1;
859
860 const size_t* nextOffsets = mConnectivity->offsetsNextZ();
861 const size_t* prevOffsets = mConnectivity->offsetsPrevZ();
862
863 if (mAxis == Y_AXIS) {
864
865 idxA = 0;
866 idxB = 2;
867 step = DIM;
868
869 nextOffsets = mConnectivity->offsetsNextY();
870 prevOffsets = mConnectivity->offsetsPrevY();
871
872 } else if (mAxis == X_AXIS) {
873
874 idxA = 1;
875 idxB = 2;
876 step = DIM*DIM;
877
878 nextOffsets = mConnectivity->offsetsNextX();
879 prevOffsets = mConnectivity->offsetsPrevX();
880 }
881
882 Coord ijk(0, 0, 0);
883
884 Int32& a = ijk[idxA];
885 Int32& b = ijk[idxB];
886
887 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
888
889 size_t startOffset = mStartNodeIndices[n];
890 size_t lastOffset = startOffset;
891
892 Int32 pos(0);
893
894 for (a = 0; a < DIM; ++a) {
895 for (b = 0; b < DIM; ++b) {
896
897 pos = static_cast<Int32>(LeafNodeType::coordToOffset(ijk));
898 size_t offset = startOffset;
899
900 // sweep in +axis direction until a boundary voxel is hit.
901 while ( offset != ConnectivityTable::INVALID_OFFSET &&
902 traceVoxelLine(*nodes[offset], pos, step) ) {
903
904 lastOffset = offset;
905 offset = nextOffsets[offset];
906 }
907
908 // find last leafnode in +axis direction
909 offset = lastOffset;
910 while (offset != ConnectivityTable::INVALID_OFFSET) {
911 lastOffset = offset;
912 offset = nextOffsets[offset];
913 }
914
915 // sweep in -axis direction until a boundary voxel is hit.
916 offset = lastOffset;
917 pos += step * (DIM - 1);
918 while ( offset != ConnectivityTable::INVALID_OFFSET &&
919 traceVoxelLine(*nodes[offset], pos, -step)) {
920 offset = prevOffsets[offset];
921 }
922 }
923 }
924 }
925 }
926
927
928 bool traceVoxelLine(LeafNodeType& node, Int32 pos, const Int32 step) const {
929
930 ValueType* data = node.buffer().data();
931
932 bool isOutside = true;
933
934 for (Index i = 0; i < LeafNodeType::DIM; ++i) {
935
936 assert(pos >= 0);
937 ValueType& dist = data[pos];
938
939 if (dist < ValueType(0.0)) {
940 isOutside = true;
941 } else {
942 // Boundary voxel check. (Voxel that intersects the surface)
943 if (!(dist > ValueType(0.75))) isOutside = false;
944
945 if (isOutside) dist = ValueType(-dist);
946 }
947
948 pos += step;
949 }
950
951 return isOutside;
952 }
953
954
955private:
956 size_t const * const mStartNodeIndices;
957 ConnectivityTable * const mConnectivity;
958
959 const Axis mAxis;
960}; // class SweepExteriorSign
961
962
963template<typename LeafNodeType>
964inline void
965seedFill(LeafNodeType& node)
966{
967 using ValueType = typename LeafNodeType::ValueType;
968 using Queue = std::deque<Index>;
969
970
971 ValueType* data = node.buffer().data();
972
973 // find seed points
974 Queue seedPoints;
975 for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
976 if (data[pos] < 0.0) seedPoints.push_back(pos);
977 }
978
979 if (seedPoints.empty()) return;
980
981 // clear sign information
982 for (Queue::iterator it = seedPoints.begin(); it != seedPoints.end(); ++it) {
983 ValueType& dist = data[*it];
984 dist = -dist;
985 }
986
987 // flood fill
988
989 Coord ijk(0, 0, 0);
990 Index pos(0), nextPos(0);
991
992 while (!seedPoints.empty()) {
993
994 pos = seedPoints.back();
995 seedPoints.pop_back();
996
997 ValueType& dist = data[pos];
998
999 if (!(dist < ValueType(0.0))) {
1000
1001 dist = -dist; // flip sign
1002
1003 ijk = LeafNodeType::offsetToLocalCoord(pos);
1004
1005 if (ijk[0] != 0) { // i - 1, j, k
1006 nextPos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
1007 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1008 }
1009
1010 if (ijk[0] != (LeafNodeType::DIM - 1)) { // i + 1, j, k
1011 nextPos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
1012 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1013 }
1014
1015 if (ijk[1] != 0) { // i, j - 1, k
1016 nextPos = pos - LeafNodeType::DIM;
1017 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1018 }
1019
1020 if (ijk[1] != (LeafNodeType::DIM - 1)) { // i, j + 1, k
1021 nextPos = pos + LeafNodeType::DIM;
1022 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1023 }
1024
1025 if (ijk[2] != 0) { // i, j, k - 1
1026 nextPos = pos - 1;
1027 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1028 }
1029
1030 if (ijk[2] != (LeafNodeType::DIM - 1)) { // i, j, k + 1
1031 nextPos = pos + 1;
1032 if (data[nextPos] > ValueType(0.75)) seedPoints.push_back(nextPos);
1033 }
1034 }
1035 }
1036} // seedFill()
1037
1038
1039template<typename LeafNodeType>
1040inline bool
1041scanFill(LeafNodeType& node)
1042{
1043 bool updatedNode = false;
1044
1045 using ValueType = typename LeafNodeType::ValueType;
1046 ValueType* data = node.buffer().data();
1047
1048 Coord ijk(0, 0, 0);
1049
1050 bool updatedSign = true;
1051 while (updatedSign) {
1052
1053 updatedSign = false;
1054
1055 for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1056
1057 ValueType& dist = data[pos];
1058
1059 if (!(dist < ValueType(0.0)) && dist > ValueType(0.75)) {
1060
1061 ijk = LeafNodeType::offsetToLocalCoord(pos);
1062
1063 // i, j, k - 1
1064 if (ijk[2] != 0 && data[pos - 1] < ValueType(0.0)) {
1065 updatedSign = true;
1066 dist = ValueType(-dist);
1067
1068 // i, j, k + 1
1069 } else if (ijk[2] != (LeafNodeType::DIM - 1) && data[pos + 1] < ValueType(0.0)) {
1070 updatedSign = true;
1071 dist = ValueType(-dist);
1072
1073 // i, j - 1, k
1074 } else if (ijk[1] != 0 && data[pos - LeafNodeType::DIM] < ValueType(0.0)) {
1075 updatedSign = true;
1076 dist = ValueType(-dist);
1077
1078 // i, j + 1, k
1079 } else if (ijk[1] != (LeafNodeType::DIM - 1)
1080 && data[pos + LeafNodeType::DIM] < ValueType(0.0))
1081 {
1082 updatedSign = true;
1083 dist = ValueType(-dist);
1084
1085 // i - 1, j, k
1086 } else if (ijk[0] != 0
1087 && data[pos - LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1088 {
1089 updatedSign = true;
1090 dist = ValueType(-dist);
1091
1092 // i + 1, j, k
1093 } else if (ijk[0] != (LeafNodeType::DIM - 1)
1094 && data[pos + LeafNodeType::DIM * LeafNodeType::DIM] < ValueType(0.0))
1095 {
1096 updatedSign = true;
1097 dist = ValueType(-dist);
1098 }
1099 }
1100 } // end value loop
1101
1102 updatedNode |= updatedSign;
1103 } // end update loop
1104
1105 return updatedNode;
1106} // scanFill()
1107
1108
1109template<typename TreeType>
1110class SeedFillExteriorSign
1111{
1112public:
1113 using ValueType = typename TreeType::ValueType;
1114 using LeafNodeType = typename TreeType::LeafNodeType;
1115
1116 SeedFillExteriorSign(std::vector<LeafNodeType*>& nodes, const bool* changedNodeMask)
1117 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1118 , mChangedNodeMask(changedNodeMask)
1119 {
1120 }
1121
1122 void operator()(const tbb::blocked_range<size_t>& range) const {
1123 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1124 if (mChangedNodeMask[n]) {
1125 //seedFill(*mNodes[n]);
1126 // Do not update the flag in mChangedNodeMask even if scanFill
1127 // returns false. mChangedNodeMask is queried by neighboring
1128 // accesses in ::SeedPoints which needs to know that this
1129 // node has values propagated on a previous iteration.
1130 scanFill(*mNodes[n]);
1131 }
1132 }
1133 }
1134
1135 LeafNodeType ** const mNodes;
1136 const bool * const mChangedNodeMask;
1137};
1138
1139
1140template<typename ValueType>
1141struct FillArray
1142{
1143 FillArray(ValueType* array, const ValueType v) : mArray(array), mValue(v) { }
1144
1145 void operator()(const tbb::blocked_range<size_t>& range) const {
1146 const ValueType v = mValue;
1147 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1148 mArray[n] = v;
1149 }
1150 }
1151
1152 ValueType * const mArray;
1153 const ValueType mValue;
1154};
1155
1156
1157template<typename ValueType>
1158inline void
1159fillArray(ValueType* array, const ValueType val, const size_t length)
1160{
1161 const auto grainSize = std::max<size_t>(
1162 length / tbb::this_task_arena::max_concurrency(), 1024);
1163 const tbb::blocked_range<size_t> range(0, length, grainSize);
1164 tbb::parallel_for(range, FillArray<ValueType>(array, val), tbb::simple_partitioner());
1165}
1166
1167
1168template<typename TreeType>
1169class SyncVoxelMask
1170{
1171public:
1172 using ValueType = typename TreeType::ValueType;
1173 using LeafNodeType = typename TreeType::LeafNodeType;
1174
1175 SyncVoxelMask(std::vector<LeafNodeType*>& nodes,
1176 const bool* changedNodeMask, bool* changedVoxelMask)
1177 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1178 , mChangedNodeMask(changedNodeMask)
1179 , mChangedVoxelMask(changedVoxelMask)
1180 {
1181 }
1182
1183 void operator()(const tbb::blocked_range<size_t>& range) const {
1184 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1185
1186 if (mChangedNodeMask[n]) {
1187 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1188
1189 ValueType* data = mNodes[n]->buffer().data();
1190
1191 for (Index pos = 0; pos < LeafNodeType::SIZE; ++pos) {
1192 if (mask[pos]) {
1193 data[pos] = ValueType(-data[pos]);
1194 mask[pos] = false;
1195 }
1196 }
1197 }
1198 }
1199 }
1200
1201 LeafNodeType ** const mNodes;
1202 bool const * const mChangedNodeMask;
1203 bool * const mChangedVoxelMask;
1204};
1205
1206
1207template<typename TreeType>
1208class SeedPoints
1209{
1210public:
1211 using ValueType = typename TreeType::ValueType;
1212 using LeafNodeType = typename TreeType::LeafNodeType;
1213 using ConnectivityTable = LeafNodeConnectivityTable<TreeType>;
1214
1215 SeedPoints(ConnectivityTable& connectivity,
1216 bool* changedNodeMask, bool* nodeMask, bool* changedVoxelMask)
1217 : mConnectivity(&connectivity)
1218 , mChangedNodeMask(changedNodeMask)
1219 , mNodeMask(nodeMask)
1220 , mChangedVoxelMask(changedVoxelMask)
1221 {
1222 }
1223
1224 void operator()(const tbb::blocked_range<size_t>& range) const {
1225
1226 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1227 bool changedValue = false;
1228
1229 changedValue |= processZ(n, /*firstFace=*/true);
1230 changedValue |= processZ(n, /*firstFace=*/false);
1231
1232 changedValue |= processY(n, /*firstFace=*/true);
1233 changedValue |= processY(n, /*firstFace=*/false);
1234
1235 changedValue |= processX(n, /*firstFace=*/true);
1236 changedValue |= processX(n, /*firstFace=*/false);
1237
1238 mNodeMask[n] = changedValue;
1239 }
1240 }
1241
1242
1243 bool processZ(const size_t n, bool firstFace) const
1244 {
1245 const size_t offset =
1246 firstFace ? mConnectivity->offsetsPrevZ()[n] : mConnectivity->offsetsNextZ()[n];
1247 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1248
1249 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1250
1251 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1252 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1253
1254 const Index lastOffset = LeafNodeType::DIM - 1;
1255 const Index lhsOffset =
1256 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1257
1258 Index tmpPos(0), pos(0);
1259 bool changedValue = false;
1260
1261 for (Index x = 0; x < LeafNodeType::DIM; ++x) {
1262 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1263 for (Index y = 0; y < LeafNodeType::DIM; ++y) {
1264 pos = tmpPos + (y << LeafNodeType::LOG2DIM);
1265
1266 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1267 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1268 changedValue = true;
1269 mask[pos + lhsOffset] = true;
1270 }
1271 }
1272 }
1273 }
1274
1275 return changedValue;
1276 }
1277
1278 return false;
1279 }
1280
1281 bool processY(const size_t n, bool firstFace) const
1282 {
1283 const size_t offset =
1284 firstFace ? mConnectivity->offsetsPrevY()[n] : mConnectivity->offsetsNextY()[n];
1285 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1286
1287 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1288
1289 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1290 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1291
1292 const Index lastOffset = LeafNodeType::DIM * (LeafNodeType::DIM - 1);
1293 const Index lhsOffset =
1294 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1295
1296 Index tmpPos(0), pos(0);
1297 bool changedValue = false;
1298
1299 for (Index x = 0; x < LeafNodeType::DIM; ++x) {
1300 tmpPos = x << (2 * LeafNodeType::LOG2DIM);
1301 for (Index z = 0; z < LeafNodeType::DIM; ++z) {
1302 pos = tmpPos + z;
1303
1304 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1305 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1306 changedValue = true;
1307 mask[pos + lhsOffset] = true;
1308 }
1309 }
1310 }
1311 }
1312
1313 return changedValue;
1314 }
1315
1316 return false;
1317 }
1318
1319 bool processX(const size_t n, bool firstFace) const
1320 {
1321 const size_t offset =
1322 firstFace ? mConnectivity->offsetsPrevX()[n] : mConnectivity->offsetsNextX()[n];
1323 if (offset != ConnectivityTable::INVALID_OFFSET && mChangedNodeMask[offset]) {
1324
1325 bool* mask = &mChangedVoxelMask[n * LeafNodeType::SIZE];
1326
1327 const ValueType* lhsData = mConnectivity->nodes()[n]->buffer().data();
1328 const ValueType* rhsData = mConnectivity->nodes()[offset]->buffer().data();
1329
1330 const Index lastOffset = LeafNodeType::DIM * LeafNodeType::DIM * (LeafNodeType::DIM-1);
1331 const Index lhsOffset =
1332 firstFace ? 0 : lastOffset, rhsOffset = firstFace ? lastOffset : 0;
1333
1334 Index tmpPos(0), pos(0);
1335 bool changedValue = false;
1336
1337 for (Index y = 0; y < LeafNodeType::DIM; ++y) {
1338 tmpPos = y << LeafNodeType::LOG2DIM;
1339 for (Index z = 0; z < LeafNodeType::DIM; ++z) {
1340 pos = tmpPos + z;
1341
1342 if (lhsData[pos + lhsOffset] > ValueType(0.75)) {
1343 if (rhsData[pos + rhsOffset] < ValueType(0.0)) {
1344 changedValue = true;
1345 mask[pos + lhsOffset] = true;
1346 }
1347 }
1348 }
1349 }
1350
1351 return changedValue;
1352 }
1353
1354 return false;
1355 }
1356
1357 ConnectivityTable * const mConnectivity;
1358 bool * const mChangedNodeMask;
1359 bool * const mNodeMask;
1360 bool * const mChangedVoxelMask;
1361};
1362
1363
1364////////////////////////////////////////
1365
1366template<typename TreeType, typename MeshDataAdapter>
1367struct ComputeIntersectingVoxelSign
1368{
1369 using ValueType = typename TreeType::ValueType;
1370 using LeafNodeType = typename TreeType::LeafNodeType;
1371 using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
1372 using Int32LeafNodeType = typename Int32TreeType::LeafNodeType;
1373
1374 using PointArray = std::unique_ptr<Vec3d[]>;
1375 using MaskArray = std::unique_ptr<bool[]>;
1376 using LocalData = std::pair<PointArray, MaskArray>;
1377 using LocalDataTable = tbb::enumerable_thread_specific<LocalData>;
1378
1379 ComputeIntersectingVoxelSign(
1380 std::vector<LeafNodeType*>& distNodes,
1381 const TreeType& distTree,
1382 const Int32TreeType& indexTree,
1383 const MeshDataAdapter& mesh)
1384 : mDistNodes(distNodes.empty() ? nullptr : &distNodes[0])
1385 , mDistTree(&distTree)
1386 , mIndexTree(&indexTree)
1387 , mMesh(&mesh)
1388 , mLocalDataTable(new LocalDataTable())
1389 {
1390 }
1391
1392
1393 void operator()(const tbb::blocked_range<size_t>& range) const {
1394
1395 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1396 tree::ValueAccessor<const Int32TreeType> idxAcc(*mIndexTree);
1397
1398 ValueType nval;
1399 CoordBBox bbox;
1400 Index xPos(0), yPos(0);
1401 Coord ijk, nijk, nodeMin, nodeMax;
1402 Vec3d cp, xyz, nxyz, dir1, dir2;
1403
1404 LocalData& localData = mLocalDataTable->local();
1405
1406 PointArray& points = localData.first;
1407 if (!points) points.reset(new Vec3d[LeafNodeType::SIZE * 2]);
1408
1409 MaskArray& mask = localData.second;
1410 if (!mask) mask.reset(new bool[LeafNodeType::SIZE]);
1411
1412
1413 typename LeafNodeType::ValueOnCIter it;
1414
1415 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1416
1417 LeafNodeType& node = *mDistNodes[n];
1418 ValueType* data = node.buffer().data();
1419
1420 const Int32LeafNodeType* idxNode = idxAcc.probeConstLeaf(node.origin());
1421 const Int32* idxData = idxNode->buffer().data();
1422
1423 nodeMin = node.origin();
1424 nodeMax = nodeMin.offsetBy(LeafNodeType::DIM - 1);
1425
1426 // reset computed voxel mask.
1427 memset(mask.get(), 0, sizeof(bool) * LeafNodeType::SIZE);
1428
1429 for (it = node.cbeginValueOn(); it; ++it) {
1430 Index pos = it.pos();
1431
1432 ValueType& dist = data[pos];
1433 if (dist < 0.0 || dist > 0.75) continue;
1434
1435 ijk = node.offsetToGlobalCoord(pos);
1436
1437 xyz[0] = double(ijk[0]);
1438 xyz[1] = double(ijk[1]);
1439 xyz[2] = double(ijk[2]);
1440
1441
1442 bbox.min() = Coord::maxComponent(ijk.offsetBy(-1), nodeMin);
1443 bbox.max() = Coord::minComponent(ijk.offsetBy(1), nodeMax);
1444
1445 bool flipSign = false;
1446
1447 for (nijk[0] = bbox.min()[0]; nijk[0] <= bbox.max()[0] && !flipSign; ++nijk[0]) {
1448 xPos = (nijk[0] & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
1449 for (nijk[1]=bbox.min()[1]; nijk[1] <= bbox.max()[1] && !flipSign; ++nijk[1]) {
1450 yPos = xPos + ((nijk[1] & (LeafNodeType::DIM-1u)) << LeafNodeType::LOG2DIM);
1451 for (nijk[2] = bbox.min()[2]; nijk[2] <= bbox.max()[2]; ++nijk[2]) {
1452 pos = yPos + (nijk[2] & (LeafNodeType::DIM - 1u));
1453
1454 const Int32& polyIdx = idxData[pos];
1455
1456 if (polyIdx == Int32(util::INVALID_IDX) || !(data[pos] < -0.75))
1457 continue;
1458
1459 const Index pointIndex = pos * 2;
1460
1461 if (!mask[pos]) {
1462
1463 mask[pos] = true;
1464
1465 nxyz[0] = double(nijk[0]);
1466 nxyz[1] = double(nijk[1]);
1467 nxyz[2] = double(nijk[2]);
1468
1469 Vec3d& point = points[pointIndex];
1470
1471 point = closestPoint(nxyz, polyIdx);
1472
1473 Vec3d& direction = points[pointIndex + 1];
1474 direction = nxyz - point;
1475 direction.normalize();
1476 }
1477
1478 dir1 = xyz - points[pointIndex];
1479 dir1.normalize();
1480
1481 if (points[pointIndex + 1].dot(dir1) > 0.0) {
1482 flipSign = true;
1483 break;
1484 }
1485 }
1486 }
1487 }
1488
1489 if (flipSign) {
1490 dist = -dist;
1491 } else {
1492 for (Int32 m = 0; m < 26; ++m) {
1493 nijk = ijk + util::COORD_OFFSETS[m];
1494
1495 if (!bbox.isInside(nijk) && distAcc.probeValue(nijk, nval) && nval<-0.75) {
1496 nxyz[0] = double(nijk[0]);
1497 nxyz[1] = double(nijk[1]);
1498 nxyz[2] = double(nijk[2]);
1499
1500 cp = closestPoint(nxyz, idxAcc.getValue(nijk));
1501
1502 dir1 = xyz - cp;
1503 dir1.normalize();
1504
1505 dir2 = nxyz - cp;
1506 dir2.normalize();
1507
1508 if (dir2.dot(dir1) > 0.0) {
1509 dist = -dist;
1510 break;
1511 }
1512 }
1513 }
1514 }
1515
1516 } // active voxel loop
1517 } // leaf node loop
1518 }
1519
1520private:
1521
1522 Vec3d closestPoint(const Vec3d& center, Int32 polyIdx) const
1523 {
1524 Vec3d a, b, c, cp, uvw;
1525
1526 const size_t polygon = size_t(polyIdx);
1527 mMesh->getIndexSpacePoint(polygon, 0, a);
1528 mMesh->getIndexSpacePoint(polygon, 1, b);
1529 mMesh->getIndexSpacePoint(polygon, 2, c);
1530
1531 cp = closestPointOnTriangleToPoint(a, c, b, center, uvw);
1532
1533 if (4 == mMesh->vertexCount(polygon)) {
1534
1535 mMesh->getIndexSpacePoint(polygon, 3, b);
1536
1537 c = closestPointOnTriangleToPoint(a, b, c, center, uvw);
1538
1539 if ((center - c).lengthSqr() < (center - cp).lengthSqr()) {
1540 cp = c;
1541 }
1542 }
1543
1544 return cp;
1545 }
1546
1547
1548 LeafNodeType ** const mDistNodes;
1549 TreeType const * const mDistTree;
1550 Int32TreeType const * const mIndexTree;
1551 MeshDataAdapter const * const mMesh;
1552
1553 SharedPtr<LocalDataTable> mLocalDataTable;
1554}; // ComputeIntersectingVoxelSign
1555
1556
1557////////////////////////////////////////
1558
1559
1560template<typename LeafNodeType>
1561inline void
1562maskNodeInternalNeighbours(const Index pos, bool (&mask)[26])
1563{
1564 using NodeT = LeafNodeType;
1565
1566 const Coord ijk = NodeT::offsetToLocalCoord(pos);
1567
1568 // Face adjacent neighbours
1569 // i+1, j, k
1570 mask[0] = ijk[0] != (NodeT::DIM - 1);
1571 // i-1, j, k
1572 mask[1] = ijk[0] != 0;
1573 // i, j+1, k
1574 mask[2] = ijk[1] != (NodeT::DIM - 1);
1575 // i, j-1, k
1576 mask[3] = ijk[1] != 0;
1577 // i, j, k+1
1578 mask[4] = ijk[2] != (NodeT::DIM - 1);
1579 // i, j, k-1
1580 mask[5] = ijk[2] != 0;
1581
1582 // Edge adjacent neighbour
1583 // i+1, j, k-1
1584 mask[6] = mask[0] && mask[5];
1585 // i-1, j, k-1
1586 mask[7] = mask[1] && mask[5];
1587 // i+1, j, k+1
1588 mask[8] = mask[0] && mask[4];
1589 // i-1, j, k+1
1590 mask[9] = mask[1] && mask[4];
1591 // i+1, j+1, k
1592 mask[10] = mask[0] && mask[2];
1593 // i-1, j+1, k
1594 mask[11] = mask[1] && mask[2];
1595 // i+1, j-1, k
1596 mask[12] = mask[0] && mask[3];
1597 // i-1, j-1, k
1598 mask[13] = mask[1] && mask[3];
1599 // i, j-1, k+1
1600 mask[14] = mask[3] && mask[4];
1601 // i, j-1, k-1
1602 mask[15] = mask[3] && mask[5];
1603 // i, j+1, k+1
1604 mask[16] = mask[2] && mask[4];
1605 // i, j+1, k-1
1606 mask[17] = mask[2] && mask[5];
1607
1608 // Corner adjacent neighbours
1609 // i-1, j-1, k-1
1610 mask[18] = mask[1] && mask[3] && mask[5];
1611 // i-1, j-1, k+1
1612 mask[19] = mask[1] && mask[3] && mask[4];
1613 // i+1, j-1, k+1
1614 mask[20] = mask[0] && mask[3] && mask[4];
1615 // i+1, j-1, k-1
1616 mask[21] = mask[0] && mask[3] && mask[5];
1617 // i-1, j+1, k-1
1618 mask[22] = mask[1] && mask[2] && mask[5];
1619 // i-1, j+1, k+1
1620 mask[23] = mask[1] && mask[2] && mask[4];
1621 // i+1, j+1, k+1
1622 mask[24] = mask[0] && mask[2] && mask[4];
1623 // i+1, j+1, k-1
1624 mask[25] = mask[0] && mask[2] && mask[5];
1625}
1626
1627
1628template<typename Compare, typename LeafNodeType>
1629inline bool
1630checkNeighbours(const Index pos, const typename LeafNodeType::ValueType * data, bool (&mask)[26])
1631{
1632 using NodeT = LeafNodeType;
1633
1634 // i, j, k - 1
1635 if (mask[5] && Compare::check(data[pos - 1])) return true;
1636 // i, j, k + 1
1637 if (mask[4] && Compare::check(data[pos + 1])) return true;
1638 // i, j - 1, k
1639 if (mask[3] && Compare::check(data[pos - NodeT::DIM])) return true;
1640 // i, j + 1, k
1641 if (mask[2] && Compare::check(data[pos + NodeT::DIM])) return true;
1642 // i - 1, j, k
1643 if (mask[1] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM])) return true;
1644 // i + 1, j, k
1645 if (mask[0] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM])) return true;
1646 // i+1, j, k-1
1647 if (mask[6] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM])) return true;
1648 // i-1, j, k-1
1649 if (mask[7] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - 1])) return true;
1650 // i+1, j, k+1
1651 if (mask[8] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + 1])) return true;
1652 // i-1, j, k+1
1653 if (mask[9] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + 1])) return true;
1654 // i+1, j+1, k
1655 if (mask[10] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM])) return true;
1656 // i-1, j+1, k
1657 if (mask[11] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM])) return true;
1658 // i+1, j-1, k
1659 if (mask[12] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM])) return true;
1660 // i-1, j-1, k
1661 if (mask[13] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM])) return true;
1662 // i, j-1, k+1
1663 if (mask[14] && Compare::check(data[pos - NodeT::DIM + 1])) return true;
1664 // i, j-1, k-1
1665 if (mask[15] && Compare::check(data[pos - NodeT::DIM - 1])) return true;
1666 // i, j+1, k+1
1667 if (mask[16] && Compare::check(data[pos + NodeT::DIM + 1])) return true;
1668 // i, j+1, k-1
1669 if (mask[17] && Compare::check(data[pos + NodeT::DIM - 1])) return true;
1670 // i-1, j-1, k-1
1671 if (mask[18] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM - 1])) return true;
1672 // i-1, j-1, k+1
1673 if (mask[19] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM - NodeT::DIM + 1])) return true;
1674 // i+1, j-1, k+1
1675 if (mask[20] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM + 1])) return true;
1676 // i+1, j-1, k-1
1677 if (mask[21] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM - NodeT::DIM - 1])) return true;
1678 // i-1, j+1, k-1
1679 if (mask[22] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM - 1])) return true;
1680 // i-1, j+1, k+1
1681 if (mask[23] && Compare::check(data[pos - NodeT::DIM * NodeT::DIM + NodeT::DIM + 1])) return true;
1682 // i+1, j+1, k+1
1683 if (mask[24] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM + 1])) return true;
1684 // i+1, j+1, k-1
1685 if (mask[25] && Compare::check(data[pos + NodeT::DIM * NodeT::DIM + NodeT::DIM - 1])) return true;
1686
1687 return false;
1688}
1689
1690
1691template<typename Compare, typename AccessorType>
1692inline bool
1693checkNeighbours(const Coord& ijk, AccessorType& acc, bool (&mask)[26])
1694{
1695 for (Int32 m = 0; m < 26; ++m) {
1696 if (!mask[m] && Compare::check(acc.getValue(ijk + util::COORD_OFFSETS[m]))) {
1697 return true;
1698 }
1699 }
1700
1701 return false;
1702}
1703
1704
1705template<typename TreeType>
1706struct ValidateIntersectingVoxels
1707{
1708 using ValueType = typename TreeType::ValueType;
1709 using LeafNodeType = typename TreeType::LeafNodeType;
1710
1711 struct IsNegative { static bool check(const ValueType v) { return v < ValueType(0.0); } };
1712
1713 ValidateIntersectingVoxels(TreeType& tree, std::vector<LeafNodeType*>& nodes)
1714 : mTree(&tree)
1715 , mNodes(nodes.empty() ? nullptr : &nodes[0])
1716 {
1717 }
1718
1719 void operator()(const tbb::blocked_range<size_t>& range) const
1720 {
1721 tree::ValueAccessor<const TreeType> acc(*mTree);
1722 bool neighbourMask[26];
1723
1724 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1725
1726 LeafNodeType& node = *mNodes[n];
1727 ValueType* data = node.buffer().data();
1728
1729 typename LeafNodeType::ValueOnCIter it;
1730 for (it = node.cbeginValueOn(); it; ++it) {
1731
1732 const Index pos = it.pos();
1733
1734 ValueType& dist = data[pos];
1735 if (dist < 0.0 || dist > 0.75) continue;
1736
1737 // Mask node internal neighbours
1738 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1739
1740 const bool hasNegativeNeighbour =
1741 checkNeighbours<IsNegative, LeafNodeType>(pos, data, neighbourMask) ||
1742 checkNeighbours<IsNegative>(node.offsetToGlobalCoord(pos), acc, neighbourMask);
1743
1744 if (!hasNegativeNeighbour) {
1745 // push over boundary voxel distance
1746 dist = ValueType(0.75) + Tolerance<ValueType>::epsilon();
1747 }
1748 }
1749 }
1750 }
1751
1752 TreeType * const mTree;
1753 LeafNodeType ** const mNodes;
1754}; // ValidateIntersectingVoxels
1755
1756
1757template<typename TreeType>
1758struct RemoveSelfIntersectingSurface
1759{
1760 using ValueType = typename TreeType::ValueType;
1761 using LeafNodeType = typename TreeType::LeafNodeType;
1762 using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
1763
1764 struct Comp { static bool check(const ValueType v) { return !(v > ValueType(0.75)); } };
1765
1766 RemoveSelfIntersectingSurface(std::vector<LeafNodeType*>& nodes,
1767 TreeType& distTree, Int32TreeType& indexTree)
1768 : mNodes(nodes.empty() ? nullptr : &nodes[0])
1769 , mDistTree(&distTree)
1770 , mIndexTree(&indexTree)
1771 {
1772 }
1773
1774 void operator()(const tbb::blocked_range<size_t>& range) const
1775 {
1776 tree::ValueAccessor<const TreeType> distAcc(*mDistTree);
1777 tree::ValueAccessor<Int32TreeType> idxAcc(*mIndexTree);
1778 bool neighbourMask[26];
1779
1780 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1781
1782 LeafNodeType& distNode = *mNodes[n];
1783 ValueType* data = distNode.buffer().data();
1784
1785 typename Int32TreeType::LeafNodeType* idxNode =
1786 idxAcc.probeLeaf(distNode.origin());
1787
1788 typename LeafNodeType::ValueOnCIter it;
1789 for (it = distNode.cbeginValueOn(); it; ++it) {
1790
1791 const Index pos = it.pos();
1792
1793 if (!(data[pos] > 0.75)) continue;
1794
1795 // Mask node internal neighbours
1796 maskNodeInternalNeighbours<LeafNodeType>(pos, neighbourMask);
1797
1798 const bool hasBoundaryNeighbour =
1799 checkNeighbours<Comp, LeafNodeType>(pos, data, neighbourMask) ||
1800 checkNeighbours<Comp>(distNode.offsetToGlobalCoord(pos),distAcc,neighbourMask);
1801
1802 if (!hasBoundaryNeighbour) {
1803 distNode.setValueOff(pos);
1804 idxNode->setValueOff(pos);
1805 }
1806 }
1807 }
1808 }
1809
1810 LeafNodeType * * const mNodes;
1811 TreeType * const mDistTree;
1812 Int32TreeType * const mIndexTree;
1813}; // RemoveSelfIntersectingSurface
1814
1815
1816////////////////////////////////////////
1817
1818
1819template<typename NodeType>
1820struct ReleaseChildNodes
1821{
1822 ReleaseChildNodes(NodeType ** nodes) : mNodes(nodes) {}
1823
1824 void operator()(const tbb::blocked_range<size_t>& range) const {
1825
1826 using NodeMaskType = typename NodeType::NodeMaskType;
1827
1828 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
1829 const_cast<NodeMaskType&>(mNodes[n]->getChildMask()).setOff();
1830 }
1831 }
1832
1833 NodeType ** const mNodes;
1834};
1835
1836
1837template<typename TreeType>
1838inline void
1839releaseLeafNodes(TreeType& tree)
1840{
1841 using RootNodeType = typename TreeType::RootNodeType;
1842 using NodeChainType = typename RootNodeType::NodeChainType;
1843 using InternalNodeType = typename NodeChainType::template Get<1>;
1844
1845 std::vector<InternalNodeType*> nodes;
1846 tree.getNodes(nodes);
1847
1848 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
1849 ReleaseChildNodes<InternalNodeType>(nodes.empty() ? nullptr : &nodes[0]));
1850}
1851
1852
1853template<typename TreeType>
1854struct StealUniqueLeafNodes
1855{
1856 using LeafNodeType = typename TreeType::LeafNodeType;
1857
1858 StealUniqueLeafNodes(TreeType& lhsTree, TreeType& rhsTree,
1859 std::vector<LeafNodeType*>& overlappingNodes)
1860 : mLhsTree(&lhsTree)
1861 , mRhsTree(&rhsTree)
1862 , mNodes(&overlappingNodes)
1863 {
1864 }
1865
1866 void operator()() const {
1867
1868 std::vector<LeafNodeType*> rhsLeafNodes;
1869
1870 rhsLeafNodes.reserve(mRhsTree->leafCount());
1871 //mRhsTree->getNodes(rhsLeafNodes);
1872 //releaseLeafNodes(*mRhsTree);
1873 mRhsTree->stealNodes(rhsLeafNodes);
1874
1875 tree::ValueAccessor<TreeType> acc(*mLhsTree);
1876
1877 for (size_t n = 0, N = rhsLeafNodes.size(); n < N; ++n) {
1878 if (!acc.probeLeaf(rhsLeafNodes[n]->origin())) {
1879 acc.addLeaf(rhsLeafNodes[n]);
1880 } else {
1881 mNodes->push_back(rhsLeafNodes[n]);
1882 }
1883 }
1884 }
1885
1886private:
1887 TreeType * const mLhsTree;
1888 TreeType * const mRhsTree;
1889 std::vector<LeafNodeType*> * const mNodes;
1890};
1891
1892
1893template<typename DistTreeType, typename IndexTreeType>
1894inline void
1895combineData(DistTreeType& lhsDist, IndexTreeType& lhsIdx,
1896 DistTreeType& rhsDist, IndexTreeType& rhsIdx)
1897{
1898 using DistLeafNodeType = typename DistTreeType::LeafNodeType;
1899 using IndexLeafNodeType = typename IndexTreeType::LeafNodeType;
1900
1901 std::vector<DistLeafNodeType*> overlappingDistNodes;
1902 std::vector<IndexLeafNodeType*> overlappingIdxNodes;
1903
1904 // Steal unique leafnodes
1905 tbb::task_group tasks;
1906 tasks.run(StealUniqueLeafNodes<DistTreeType>(lhsDist, rhsDist, overlappingDistNodes));
1907 tasks.run(StealUniqueLeafNodes<IndexTreeType>(lhsIdx, rhsIdx, overlappingIdxNodes));
1908 tasks.wait();
1909
1910 // Combine overlapping leaf nodes
1911 if (!overlappingDistNodes.empty() && !overlappingIdxNodes.empty()) {
1912 tbb::parallel_for(tbb::blocked_range<size_t>(0, overlappingDistNodes.size()),
1913 CombineLeafNodes<DistTreeType>(lhsDist, lhsIdx,
1914 &overlappingDistNodes[0], &overlappingIdxNodes[0]));
1915 }
1916}
1917
1918/// @brief TBB body object to voxelize a mesh of triangles and/or quads into a collection
1919/// of VDB grids, namely a squared distance grid, a closest primitive grid and an
1920/// intersecting voxels grid (masks the mesh intersecting voxels)
1921/// @note Only the leaf nodes that intersect the mesh are allocated, and only voxels in
1922/// a narrow band (of two to three voxels in proximity to the mesh's surface) are activated.
1923/// They are populated with distance values and primitive indices.
1924template<typename TreeType>
1925struct VoxelizationData {
1926
1927 using Ptr = std::unique_ptr<VoxelizationData>;
1928 using ValueType = typename TreeType::ValueType;
1929
1930 using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
1931 using UCharTreeType = typename TreeType::template ValueConverter<unsigned char>::Type;
1932
1933 using FloatTreeAcc = tree::ValueAccessor<TreeType>;
1934 using Int32TreeAcc = tree::ValueAccessor<Int32TreeType>;
1935 using UCharTreeAcc = tree::ValueAccessor<UCharTreeType>;
1936
1937
1938 VoxelizationData()
1939 : distTree(std::numeric_limits<ValueType>::max())
1940 , distAcc(distTree)
1941 , indexTree(Int32(util::INVALID_IDX))
1942 , indexAcc(indexTree)
1943 , primIdTree(MaxPrimId)
1944 , primIdAcc(primIdTree)
1945 , mPrimCount(0)
1946 {
1947 }
1948
1949 TreeType distTree;
1950 FloatTreeAcc distAcc;
1951
1952 Int32TreeType indexTree;
1953 Int32TreeAcc indexAcc;
1954
1955 UCharTreeType primIdTree;
1956 UCharTreeAcc primIdAcc;
1957
1958 unsigned char getNewPrimId() {
1959
1960 /// @warning Don't use parallel methods here!
1961 /// The primIdTree is used as a "scratch" pad to mark visits for a given polygon
1962 /// into voxels which it may contribute to. The tree is kept as lightweight as
1963 /// possible and is reset when a maximum count or size is reached. A previous
1964 /// bug here occurred due to the calling of tree methods with multi-threaded
1965 /// implementations, resulting in nested parallelization and re-use of the TLS
1966 /// from the initial task. This consequently resulted in non deterministic values
1967 /// of mPrimCount on the return of the initial task, and could potentially end up
1968 /// with a mPrimCount equal to that of the MaxPrimId. This is used as the background
1969 /// value of the scratch tree.
1970 /// @see jira.aswf.io/browse/OVDB-117, PR #564
1971 /// @todo Consider profiling this operator with tree.clear() and Investigate the
1972 /// chosen value of MaxPrimId
1973
1974 if (mPrimCount == MaxPrimId || primIdTree.leafCount() > 1000) {
1975 mPrimCount = 0;
1976 primIdTree.root().clear();
1977 primIdTree.clearAllAccessors();
1978 assert(mPrimCount == 0);
1979 }
1980
1981 return mPrimCount++;
1982 }
1983
1984private:
1985
1986 enum { MaxPrimId = 100 };
1987
1988 unsigned char mPrimCount;
1989};
1990
1991
1992template<typename TreeType, typename MeshDataAdapter, typename Interrupter = util::NullInterrupter>
1993class VoxelizePolygons
1994{
1995public:
1996
1997 using VoxelizationDataType = VoxelizationData<TreeType>;
1998 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
1999
2000 VoxelizePolygons(DataTable& dataTable,
2001 const MeshDataAdapter& mesh,
2002 Interrupter* interrupter = nullptr)
2003 : mDataTable(&dataTable)
2004 , mMesh(&mesh)
2005 , mInterrupter(interrupter)
2006 {
2007 }
2008
2009 void operator()(const tbb::blocked_range<size_t>& range) const {
2010
2011 typename VoxelizationDataType::Ptr& dataPtr = mDataTable->local();
2012 if (!dataPtr) dataPtr.reset(new VoxelizationDataType());
2013
2014 Triangle prim;
2015
2016 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2017
2018 if (this->wasInterrupted()) {
2019 thread::cancelGroupExecution();
2020 break;
2021 }
2022
2023 const size_t numVerts = mMesh->vertexCount(n);
2024
2025 // rasterize triangles and quads.
2026 if (numVerts == 3 || numVerts == 4) {
2027
2028 prim.index = Int32(n);
2029
2030 mMesh->getIndexSpacePoint(n, 0, prim.a);
2031 mMesh->getIndexSpacePoint(n, 1, prim.b);
2032 mMesh->getIndexSpacePoint(n, 2, prim.c);
2033
2034 evalTriangle(prim, *dataPtr);
2035
2036 if (numVerts == 4) {
2037 mMesh->getIndexSpacePoint(n, 3, prim.b);
2038 evalTriangle(prim, *dataPtr);
2039 }
2040 }
2041 }
2042 }
2043
2044private:
2045
2046 bool wasInterrupted() const { return mInterrupter && mInterrupter->wasInterrupted(); }
2047
2048 struct Triangle { Vec3d a, b, c; Int32 index; };
2049
2050 struct SubTask
2051 {
2052 enum { POLYGON_LIMIT = 1000 };
2053
2054 SubTask(const Triangle& prim, DataTable& dataTable,
2055 int subdivisionCount, size_t polygonCount,
2056 Interrupter* interrupter = nullptr)
2057 : mLocalDataTable(&dataTable)
2058 , mPrim(prim)
2059 , mSubdivisionCount(subdivisionCount)
2060 , mPolygonCount(polygonCount)
2061 , mInterrupter(interrupter)
2062 {
2063 }
2064
2065 void operator()() const
2066 {
2067 if (mSubdivisionCount <= 0 || mPolygonCount >= POLYGON_LIMIT) {
2068
2069 typename VoxelizationDataType::Ptr& dataPtr = mLocalDataTable->local();
2070 if (!dataPtr) dataPtr.reset(new VoxelizationDataType());
2071
2072 voxelizeTriangle(mPrim, *dataPtr, mInterrupter);
2073
2074 } else if (!(mInterrupter && mInterrupter->wasInterrupted())) {
2075 spawnTasks(mPrim, *mLocalDataTable, mSubdivisionCount, mPolygonCount, mInterrupter);
2076 }
2077 }
2078
2079 DataTable * const mLocalDataTable;
2080 Triangle const mPrim;
2081 int const mSubdivisionCount;
2082 size_t const mPolygonCount;
2083 Interrupter * const mInterrupter;
2084 }; // struct SubTask
2085
2086 inline static int evalSubdivisionCount(const Triangle& prim)
2087 {
2088 const double ax = prim.a[0], bx = prim.b[0], cx = prim.c[0];
2089 const double dx = std::max(ax, std::max(bx, cx)) - std::min(ax, std::min(bx, cx));
2090
2091 const double ay = prim.a[1], by = prim.b[1], cy = prim.c[1];
2092 const double dy = std::max(ay, std::max(by, cy)) - std::min(ay, std::min(by, cy));
2093
2094 const double az = prim.a[2], bz = prim.b[2], cz = prim.c[2];
2095 const double dz = std::max(az, std::max(bz, cz)) - std::min(az, std::min(bz, cz));
2096
2097 return int(std::max(dx, std::max(dy, dz)) / double(TreeType::LeafNodeType::DIM * 2));
2098 }
2099
2100 void evalTriangle(const Triangle& prim, VoxelizationDataType& data) const
2101 {
2102 const size_t polygonCount = mMesh->polygonCount();
2103 const int subdivisionCount =
2104 polygonCount < SubTask::POLYGON_LIMIT ? evalSubdivisionCount(prim) : 0;
2105
2106 if (subdivisionCount <= 0) {
2107 voxelizeTriangle(prim, data, mInterrupter);
2108 } else {
2109 spawnTasks(prim, *mDataTable, subdivisionCount, polygonCount, mInterrupter);
2110 }
2111 }
2112
2113 static void spawnTasks(
2114 const Triangle& mainPrim,
2115 DataTable& dataTable,
2116 int subdivisionCount,
2117 size_t polygonCount,
2118 Interrupter* const interrupter)
2119 {
2120 subdivisionCount -= 1;
2121 polygonCount *= 4;
2122
2123 tbb::task_group tasks;
2124
2125 const Vec3d ac = (mainPrim.a + mainPrim.c) * 0.5;
2126 const Vec3d bc = (mainPrim.b + mainPrim.c) * 0.5;
2127 const Vec3d ab = (mainPrim.a + mainPrim.b) * 0.5;
2128
2129 Triangle prim;
2130 prim.index = mainPrim.index;
2131
2132 prim.a = mainPrim.a;
2133 prim.b = ab;
2134 prim.c = ac;
2135 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2136
2137 prim.a = ab;
2138 prim.b = bc;
2139 prim.c = ac;
2140 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2141
2142 prim.a = ab;
2143 prim.b = mainPrim.b;
2144 prim.c = bc;
2145 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2146
2147 prim.a = ac;
2148 prim.b = bc;
2149 prim.c = mainPrim.c;
2150 tasks.run(SubTask(prim, dataTable, subdivisionCount, polygonCount, interrupter));
2151
2152 tasks.wait();
2153 }
2154
2155 static void voxelizeTriangle(const Triangle& prim, VoxelizationDataType& data, Interrupter* const interrupter)
2156 {
2157 std::deque<Coord> coordList;
2158 Coord ijk, nijk;
2159
2160 ijk = Coord::floor(prim.a);
2161 coordList.push_back(ijk);
2162
2163 // The first point may not be quite in bounds, and rely
2164 // on one of the neighbours to have the first valid seed,
2165 // so we cannot early-exit here.
2166 updateDistance(ijk, prim, data);
2167
2168 unsigned char primId = data.getNewPrimId();
2169 data.primIdAcc.setValueOnly(ijk, primId);
2170
2171 while (!coordList.empty()) {
2172 if (interrupter && interrupter->wasInterrupted()) {
2173 thread::cancelGroupExecution();
2174 break;
2175 }
2176 for (Int32 pass = 0; pass < 1048576 && !coordList.empty(); ++pass) {
2177 ijk = coordList.back();
2178 coordList.pop_back();
2179
2180 for (Int32 i = 0; i < 26; ++i) {
2181 nijk = ijk + util::COORD_OFFSETS[i];
2182 if (primId != data.primIdAcc.getValue(nijk)) {
2183 data.primIdAcc.setValueOnly(nijk, primId);
2184 if(updateDistance(nijk, prim, data)) coordList.push_back(nijk);
2185 }
2186 }
2187 }
2188 }
2189 }
2190
2191 static bool updateDistance(const Coord& ijk, const Triangle& prim, VoxelizationDataType& data)
2192 {
2193 Vec3d uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2194
2195 using ValueType = typename TreeType::ValueType;
2196
2197 const ValueType dist = ValueType((voxelCenter -
2198 closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, voxelCenter, uvw)).lengthSqr());
2199
2200 // Either the points may be NAN, or they could be far enough from
2201 // the origin that computing distance fails.
2202 if (std::isnan(dist))
2203 return false;
2204
2205 const ValueType oldDist = data.distAcc.getValue(ijk);
2206
2207 if (dist < oldDist) {
2208 data.distAcc.setValue(ijk, dist);
2209 data.indexAcc.setValue(ijk, prim.index);
2210 } else if (math::isExactlyEqual(dist, oldDist)) {
2211 // makes reduction deterministic when different polygons
2212 // produce the same distance value.
2213 data.indexAcc.setValueOnly(ijk, std::min(prim.index, data.indexAcc.getValue(ijk)));
2214 }
2215
2216 return !(dist > 0.75); // true if the primitive intersects the voxel.
2217 }
2218
2219 DataTable * const mDataTable;
2220 MeshDataAdapter const * const mMesh;
2221 Interrupter * const mInterrupter;
2222}; // VoxelizePolygons
2223
2224
2225////////////////////////////////////////
2226
2227
2228template<typename TreeType>
2229struct DiffLeafNodeMask
2230{
2231 using AccessorType = typename tree::ValueAccessor<TreeType>;
2232 using LeafNodeType = typename TreeType::LeafNodeType;
2233
2234 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2235 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2236
2237 DiffLeafNodeMask(const TreeType& rhsTree,
2238 std::vector<BoolLeafNodeType*>& lhsNodes)
2239 : mRhsTree(&rhsTree), mLhsNodes(lhsNodes.empty() ? nullptr : &lhsNodes[0])
2240 {
2241 }
2242
2243 void operator()(const tbb::blocked_range<size_t>& range) const {
2244
2245 tree::ValueAccessor<const TreeType> acc(*mRhsTree);
2246
2247 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2248
2249 BoolLeafNodeType* lhsNode = mLhsNodes[n];
2250 const LeafNodeType* rhsNode = acc.probeConstLeaf(lhsNode->origin());
2251
2252 if (rhsNode) lhsNode->topologyDifference(*rhsNode, false);
2253 }
2254 }
2255
2256private:
2257 TreeType const * const mRhsTree;
2258 BoolLeafNodeType ** const mLhsNodes;
2259};
2260
2261
2262template<typename LeafNodeTypeA, typename LeafNodeTypeB>
2263struct UnionValueMasks
2264{
2265 UnionValueMasks(std::vector<LeafNodeTypeA*>& nodesA, std::vector<LeafNodeTypeB*>& nodesB)
2266 : mNodesA(nodesA.empty() ? nullptr : &nodesA[0])
2267 , mNodesB(nodesB.empty() ? nullptr : &nodesB[0])
2268 {
2269 }
2270
2271 void operator()(const tbb::blocked_range<size_t>& range) const {
2272 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2273 mNodesA[n]->topologyUnion(*mNodesB[n]);
2274 }
2275 }
2276
2277private:
2278 LeafNodeTypeA ** const mNodesA;
2279 LeafNodeTypeB ** const mNodesB;
2280};
2281
2282
2283template<typename TreeType>
2284struct ConstructVoxelMask
2285{
2286 using LeafNodeType = typename TreeType::LeafNodeType;
2287
2288 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2289 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2290
2291 ConstructVoxelMask(BoolTreeType& maskTree, const TreeType& tree,
2292 std::vector<LeafNodeType*>& nodes)
2293 : mTree(&tree)
2294 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2295 , mLocalMaskTree(false)
2296 , mMaskTree(&maskTree)
2297 {
2298 }
2299
2300 ConstructVoxelMask(ConstructVoxelMask& rhs, tbb::split)
2301 : mTree(rhs.mTree)
2302 , mNodes(rhs.mNodes)
2303 , mLocalMaskTree(false)
2304 , mMaskTree(&mLocalMaskTree)
2305 {
2306 }
2307
2308 void operator()(const tbb::blocked_range<size_t>& range)
2309 {
2310 using Iterator = typename LeafNodeType::ValueOnCIter;
2311
2312 tree::ValueAccessor<const TreeType> acc(*mTree);
2313 tree::ValueAccessor<BoolTreeType> maskAcc(*mMaskTree);
2314
2315 Coord ijk, nijk, localCorod;
2316 Index pos, npos;
2317
2318 for (size_t n = range.begin(); n != range.end(); ++n) {
2319
2320 LeafNodeType& node = *mNodes[n];
2321
2322 CoordBBox bbox = node.getNodeBoundingBox();
2323 bbox.expand(-1);
2324
2325 BoolLeafNodeType& maskNode = *maskAcc.touchLeaf(node.origin());
2326
2327 for (Iterator it = node.cbeginValueOn(); it; ++it) {
2328 ijk = it.getCoord();
2329 pos = it.pos();
2330
2331 localCorod = LeafNodeType::offsetToLocalCoord(pos);
2332
2333 if (localCorod[2] < int(LeafNodeType::DIM - 1)) {
2334 npos = pos + 1;
2335 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2336 } else {
2337 nijk = ijk.offsetBy(0, 0, 1);
2338 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2339 }
2340
2341 if (localCorod[2] > 0) {
2342 npos = pos - 1;
2343 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2344 } else {
2345 nijk = ijk.offsetBy(0, 0, -1);
2346 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2347 }
2348
2349 if (localCorod[1] < int(LeafNodeType::DIM - 1)) {
2350 npos = pos + LeafNodeType::DIM;
2351 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2352 } else {
2353 nijk = ijk.offsetBy(0, 1, 0);
2354 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2355 }
2356
2357 if (localCorod[1] > 0) {
2358 npos = pos - LeafNodeType::DIM;
2359 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2360 } else {
2361 nijk = ijk.offsetBy(0, -1, 0);
2362 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2363 }
2364
2365 if (localCorod[0] < int(LeafNodeType::DIM - 1)) {
2366 npos = pos + LeafNodeType::DIM * LeafNodeType::DIM;
2367 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2368 } else {
2369 nijk = ijk.offsetBy(1, 0, 0);
2370 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2371 }
2372
2373 if (localCorod[0] > 0) {
2374 npos = pos - LeafNodeType::DIM * LeafNodeType::DIM;
2375 if (!node.isValueOn(npos)) maskNode.setValueOn(npos);
2376 } else {
2377 nijk = ijk.offsetBy(-1, 0, 0);
2378 if (!acc.isValueOn(nijk)) maskAcc.setValueOn(nijk);
2379 }
2380 }
2381 }
2382 }
2383
2384 void join(ConstructVoxelMask& rhs) { mMaskTree->merge(*rhs.mMaskTree); }
2385
2386private:
2387 TreeType const * const mTree;
2388 LeafNodeType ** const mNodes;
2389
2390 BoolTreeType mLocalMaskTree;
2391 BoolTreeType * const mMaskTree;
2392};
2393
2394
2395/// @note The interior and exterior widths should be in world space units and squared.
2396template<typename TreeType, typename MeshDataAdapter>
2397struct ExpandNarrowband
2398{
2399 using ValueType = typename TreeType::ValueType;
2400 using LeafNodeType = typename TreeType::LeafNodeType;
2401 using NodeMaskType = typename LeafNodeType::NodeMaskType;
2402 using Int32TreeType = typename TreeType::template ValueConverter<Int32>::Type;
2403 using Int32LeafNodeType = typename Int32TreeType::LeafNodeType;
2404 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
2405 using BoolLeafNodeType = typename BoolTreeType::LeafNodeType;
2406
2407 struct Fragment
2408 {
2409 Int32 idx, x, y, z;
2410 ValueType dist;
2411
2412 Fragment() : idx(0), x(0), y(0), z(0), dist(0.0) {}
2413
2414 Fragment(Int32 idx_, Int32 x_, Int32 y_, Int32 z_, ValueType dist_)
2415 : idx(idx_), x(x_), y(y_), z(z_), dist(dist_)
2416 {
2417 }
2418
2419 bool operator<(const Fragment& rhs) const { return idx < rhs.idx; }
2420 }; // struct Fragment
2421
2422 ////////////////////
2423
2424 ExpandNarrowband(
2425 std::vector<BoolLeafNodeType*>& maskNodes,
2426 BoolTreeType& maskTree,
2427 TreeType& distTree,
2428 Int32TreeType& indexTree,
2429 const MeshDataAdapter& mesh,
2430 ValueType exteriorBandWidth,
2431 ValueType interiorBandWidth,
2432 ValueType voxelSize)
2433 : mMaskNodes(maskNodes.empty() ? nullptr : &maskNodes[0])
2434 , mMaskTree(&maskTree)
2435 , mDistTree(&distTree)
2436 , mIndexTree(&indexTree)
2437 , mMesh(&mesh)
2438 , mNewMaskTree(false)
2439 , mDistNodes()
2440 , mUpdatedDistNodes()
2441 , mIndexNodes()
2442 , mUpdatedIndexNodes()
2443 , mExteriorBandWidth(exteriorBandWidth)
2444 , mInteriorBandWidth(interiorBandWidth)
2445 , mVoxelSize(voxelSize)
2446 {
2447 }
2448
2449 ExpandNarrowband(const ExpandNarrowband& rhs, tbb::split)
2450 : mMaskNodes(rhs.mMaskNodes)
2451 , mMaskTree(rhs.mMaskTree)
2452 , mDistTree(rhs.mDistTree)
2453 , mIndexTree(rhs.mIndexTree)
2454 , mMesh(rhs.mMesh)
2455 , mNewMaskTree(false)
2456 , mDistNodes()
2457 , mUpdatedDistNodes()
2458 , mIndexNodes()
2459 , mUpdatedIndexNodes()
2460 , mExteriorBandWidth(rhs.mExteriorBandWidth)
2461 , mInteriorBandWidth(rhs.mInteriorBandWidth)
2462 , mVoxelSize(rhs.mVoxelSize)
2463 {
2464 }
2465
2466 void join(ExpandNarrowband& rhs)
2467 {
2468 mDistNodes.insert(mDistNodes.end(), rhs.mDistNodes.begin(), rhs.mDistNodes.end());
2469 mIndexNodes.insert(mIndexNodes.end(), rhs.mIndexNodes.begin(), rhs.mIndexNodes.end());
2470
2471 mUpdatedDistNodes.insert(mUpdatedDistNodes.end(),
2472 rhs.mUpdatedDistNodes.begin(), rhs.mUpdatedDistNodes.end());
2473
2474 mUpdatedIndexNodes.insert(mUpdatedIndexNodes.end(),
2475 rhs.mUpdatedIndexNodes.begin(), rhs.mUpdatedIndexNodes.end());
2476
2477 mNewMaskTree.merge(rhs.mNewMaskTree);
2478 }
2479
2480 void operator()(const tbb::blocked_range<size_t>& range)
2481 {
2482 tree::ValueAccessor<BoolTreeType> newMaskAcc(mNewMaskTree);
2483 tree::ValueAccessor<TreeType> distAcc(*mDistTree);
2484 tree::ValueAccessor<Int32TreeType> indexAcc(*mIndexTree);
2485
2486 std::vector<Fragment> fragments;
2487 fragments.reserve(256);
2488
2489 std::unique_ptr<LeafNodeType> newDistNodePt;
2490 std::unique_ptr<Int32LeafNodeType> newIndexNodePt;
2491
2492 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2493
2494 BoolLeafNodeType& maskNode = *mMaskNodes[n];
2495 if (maskNode.isEmpty()) continue;
2496
2497 // Setup local caches
2498
2499 const Coord& origin = maskNode.origin();
2500
2501 LeafNodeType * distNodePt = distAcc.probeLeaf(origin);
2502 Int32LeafNodeType * indexNodePt = indexAcc.probeLeaf(origin);
2503
2504 assert(!distNodePt == !indexNodePt);
2505
2506 bool usingNewNodes = false;
2507
2508 if (!distNodePt && !indexNodePt) {
2509
2510 const ValueType backgroundDist = distAcc.getValue(origin);
2511
2512 if (!newDistNodePt.get() && !newIndexNodePt.get()) {
2513 newDistNodePt.reset(new LeafNodeType(origin, backgroundDist));
2514 newIndexNodePt.reset(new Int32LeafNodeType(origin, indexAcc.getValue(origin)));
2515 } else {
2516
2517 if ((backgroundDist < ValueType(0.0)) !=
2518 (newDistNodePt->getValue(0) < ValueType(0.0))) {
2519 newDistNodePt->buffer().fill(backgroundDist);
2520 }
2521
2522 newDistNodePt->setOrigin(origin);
2523 newIndexNodePt->setOrigin(origin);
2524 }
2525
2526 distNodePt = newDistNodePt.get();
2527 indexNodePt = newIndexNodePt.get();
2528
2529 usingNewNodes = true;
2530 }
2531
2532
2533 // Gather neighbour information
2534
2535 CoordBBox bbox(Coord::max(), Coord::min());
2536 for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2537 bbox.expand(it.getCoord());
2538 }
2539
2540 bbox.expand(1);
2541
2542 gatherFragments(fragments, bbox, distAcc, indexAcc);
2543
2544
2545 // Compute first voxel layer
2546
2547 bbox = maskNode.getNodeBoundingBox();
2548 NodeMaskType mask;
2549 bool updatedLeafNodes = false;
2550
2551 for (typename BoolLeafNodeType::ValueOnIter it = maskNode.beginValueOn(); it; ++it) {
2552
2553 const Coord ijk = it.getCoord();
2554
2555 if (updateVoxel(ijk, 5, fragments, *distNodePt, *indexNodePt, &updatedLeafNodes)) {
2556
2557 for (Int32 i = 0; i < 6; ++i) {
2558 const Coord nijk = ijk + util::COORD_OFFSETS[i];
2559 if (bbox.isInside(nijk)) {
2560 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2561 } else {
2562 newMaskAcc.setValueOn(nijk);
2563 }
2564 }
2565
2566 for (Int32 i = 6; i < 26; ++i) {
2567 const Coord nijk = ijk + util::COORD_OFFSETS[i];
2568 if (bbox.isInside(nijk)) {
2569 mask.setOn(BoolLeafNodeType::coordToOffset(nijk));
2570 }
2571 }
2572 }
2573 }
2574
2575 if (updatedLeafNodes) {
2576
2577 // Compute second voxel layer
2578 mask -= indexNodePt->getValueMask();
2579
2580 for (typename NodeMaskType::OnIterator it = mask.beginOn(); it; ++it) {
2581
2582 const Index pos = it.pos();
2583 const Coord ijk = maskNode.origin() + LeafNodeType::offsetToLocalCoord(pos);
2584
2585 if (updateVoxel(ijk, 6, fragments, *distNodePt, *indexNodePt)) {
2586 for (Int32 i = 0; i < 6; ++i) {
2587 newMaskAcc.setValueOn(ijk + util::COORD_OFFSETS[i]);
2588 }
2589 }
2590 }
2591
2592 // Export new distance values
2593 if (usingNewNodes) {
2594 newDistNodePt->topologyUnion(*newIndexNodePt);
2595 mDistNodes.push_back(newDistNodePt.release());
2596 mIndexNodes.push_back(newIndexNodePt.release());
2597 } else {
2598 mUpdatedDistNodes.push_back(distNodePt);
2599 mUpdatedIndexNodes.push_back(indexNodePt);
2600 }
2601 }
2602 } // end leafnode loop
2603 }
2604
2605 //////////
2606
2607 BoolTreeType& newMaskTree() { return mNewMaskTree; }
2608
2609 std::vector<LeafNodeType*>& newDistNodes() { return mDistNodes; }
2610 std::vector<LeafNodeType*>& updatedDistNodes() { return mUpdatedDistNodes; }
2611
2612 std::vector<Int32LeafNodeType*>& newIndexNodes() { return mIndexNodes; }
2613 std::vector<Int32LeafNodeType*>& updatedIndexNodes() { return mUpdatedIndexNodes; }
2614
2615private:
2616
2617 /// @note The output fragment list is ordered by the primitive index
2618 void
2619 gatherFragments(std::vector<Fragment>& fragments, const CoordBBox& bbox,
2620 tree::ValueAccessor<TreeType>& distAcc, tree::ValueAccessor<Int32TreeType>& indexAcc)
2621 {
2622 fragments.clear();
2623 const Coord nodeMin = bbox.min() & ~(LeafNodeType::DIM - 1);
2624 const Coord nodeMax = bbox.max() & ~(LeafNodeType::DIM - 1);
2625
2626 CoordBBox region;
2627 Coord ijk;
2628
2629 for (ijk[0] = nodeMin[0]; ijk[0] <= nodeMax[0]; ijk[0] += LeafNodeType::DIM) {
2630 for (ijk[1] = nodeMin[1]; ijk[1] <= nodeMax[1]; ijk[1] += LeafNodeType::DIM) {
2631 for (ijk[2] = nodeMin[2]; ijk[2] <= nodeMax[2]; ijk[2] += LeafNodeType::DIM) {
2632 if (LeafNodeType* distleaf = distAcc.probeLeaf(ijk)) {
2633 region.min() = Coord::maxComponent(bbox.min(), ijk);
2634 region.max() = Coord::minComponent(bbox.max(),
2635 ijk.offsetBy(LeafNodeType::DIM - 1));
2636 gatherFragments(fragments, region, *distleaf, *indexAcc.probeLeaf(ijk));
2637 }
2638 }
2639 }
2640 }
2641
2642 std::sort(fragments.begin(), fragments.end());
2643 }
2644
2645 void
2646 gatherFragments(std::vector<Fragment>& fragments, const CoordBBox& bbox,
2647 const LeafNodeType& distLeaf, const Int32LeafNodeType& idxLeaf) const
2648 {
2649 const typename LeafNodeType::NodeMaskType& mask = distLeaf.getValueMask();
2650 const ValueType* distData = distLeaf.buffer().data();
2651 const Int32* idxData = idxLeaf.buffer().data();
2652
2653 for (int x = bbox.min()[0]; x <= bbox.max()[0]; ++x) {
2654 const Index xPos = (x & (LeafNodeType::DIM - 1u)) << (2 * LeafNodeType::LOG2DIM);
2655 for (int y = bbox.min()[1]; y <= bbox.max()[1]; ++y) {
2656 const Index yPos = xPos + ((y & (LeafNodeType::DIM - 1u)) << LeafNodeType::LOG2DIM);
2657 for (int z = bbox.min()[2]; z <= bbox.max()[2]; ++z) {
2658 const Index pos = yPos + (z & (LeafNodeType::DIM - 1u));
2659 if (mask.isOn(pos)) {
2660 fragments.push_back(Fragment(idxData[pos],x,y,z, std::abs(distData[pos])));
2661 }
2662 }
2663 }
2664 }
2665 }
2666
2667 /// @note This method expects the fragment list to be ordered by the primitive index
2668 /// to avoid redundant distance computations.
2669 ValueType
2670 computeDistance(const Coord& ijk, const Int32 manhattanLimit,
2671 const std::vector<Fragment>& fragments, Int32& closestPrimIdx) const
2672 {
2673 Vec3d a, b, c, uvw, voxelCenter(ijk[0], ijk[1], ijk[2]);
2674 double primDist, tmpDist, dist = std::numeric_limits<double>::max();
2675 Int32 lastIdx = Int32(util::INVALID_IDX);
2676
2677 for (size_t n = 0, N = fragments.size(); n < N; ++n) {
2678
2679 const Fragment& fragment = fragments[n];
2680 if (lastIdx == fragment.idx) continue;
2681
2682 const Int32 dx = std::abs(fragment.x - ijk[0]);
2683 const Int32 dy = std::abs(fragment.y - ijk[1]);
2684 const Int32 dz = std::abs(fragment.z - ijk[2]);
2685
2686 const Int32 manhattan = dx + dy + dz;
2687 if (manhattan > manhattanLimit) continue;
2688
2689 lastIdx = fragment.idx;
2690
2691 const size_t polygon = size_t(lastIdx);
2692
2693 mMesh->getIndexSpacePoint(polygon, 0, a);
2694 mMesh->getIndexSpacePoint(polygon, 1, b);
2695 mMesh->getIndexSpacePoint(polygon, 2, c);
2696
2697 primDist = (voxelCenter -
2698 closestPointOnTriangleToPoint(a, c, b, voxelCenter, uvw)).lengthSqr();
2699
2700 // Split quad into a second triangle
2701 if (4 == mMesh->vertexCount(polygon)) {
2702
2703 mMesh->getIndexSpacePoint(polygon, 3, b);
2704
2705 tmpDist = (voxelCenter - closestPointOnTriangleToPoint(
2706 a, b, c, voxelCenter, uvw)).lengthSqr();
2707
2708 if (tmpDist < primDist) primDist = tmpDist;
2709 }
2710
2711 if (primDist < dist) {
2712 dist = primDist;
2713 closestPrimIdx = lastIdx;
2714 }
2715 }
2716
2717 return ValueType(std::sqrt(dist)) * mVoxelSize;
2718 }
2719
2720 /// @note Returns true if the current voxel was updated and neighbouring
2721 /// voxels need to be evaluated.
2722 bool
2723 updateVoxel(const Coord& ijk, const Int32 manhattanLimit,
2724 const std::vector<Fragment>& fragments,
2725 LeafNodeType& distLeaf, Int32LeafNodeType& idxLeaf, bool* updatedLeafNodes = nullptr)
2726 {
2727 Int32 closestPrimIdx = 0;
2728 const ValueType distance = computeDistance(ijk, manhattanLimit, fragments, closestPrimIdx);
2729
2730 const Index pos = LeafNodeType::coordToOffset(ijk);
2731 const bool inside = distLeaf.getValue(pos) < ValueType(0.0);
2732
2733 bool activateNeighbourVoxels = false;
2734
2735 if (!inside && distance < mExteriorBandWidth) {
2736 if (updatedLeafNodes) *updatedLeafNodes = true;
2737 activateNeighbourVoxels = (distance + mVoxelSize) < mExteriorBandWidth;
2738 distLeaf.setValueOnly(pos, distance);
2739 idxLeaf.setValueOn(pos, closestPrimIdx);
2740 } else if (inside && distance < mInteriorBandWidth) {
2741 if (updatedLeafNodes) *updatedLeafNodes = true;
2742 activateNeighbourVoxels = (distance + mVoxelSize) < mInteriorBandWidth;
2743 distLeaf.setValueOnly(pos, -distance);
2744 idxLeaf.setValueOn(pos, closestPrimIdx);
2745 }
2746
2747 return activateNeighbourVoxels;
2748 }
2749
2750 //////////
2751
2752 BoolLeafNodeType ** const mMaskNodes;
2753 BoolTreeType * const mMaskTree;
2754 TreeType * const mDistTree;
2755 Int32TreeType * const mIndexTree;
2756
2757 MeshDataAdapter const * const mMesh;
2758
2759 BoolTreeType mNewMaskTree;
2760
2761 std::vector<LeafNodeType*> mDistNodes, mUpdatedDistNodes;
2762 std::vector<Int32LeafNodeType*> mIndexNodes, mUpdatedIndexNodes;
2763
2764 const ValueType mExteriorBandWidth, mInteriorBandWidth, mVoxelSize;
2765}; // struct ExpandNarrowband
2766
2767
2768template<typename TreeType>
2769struct AddNodes {
2770 using LeafNodeType = typename TreeType::LeafNodeType;
2771
2772 AddNodes(TreeType& tree, std::vector<LeafNodeType*>& nodes)
2773 : mTree(&tree) , mNodes(&nodes)
2774 {
2775 }
2776
2777 void operator()() const {
2778 tree::ValueAccessor<TreeType> acc(*mTree);
2779 std::vector<LeafNodeType*>& nodes = *mNodes;
2780 for (size_t n = 0, N = nodes.size(); n < N; ++n) {
2781 acc.addLeaf(nodes[n]);
2782 }
2783 }
2784
2785 TreeType * const mTree;
2786 std::vector<LeafNodeType*> * const mNodes;
2787}; // AddNodes
2788
2789
2790template<typename TreeType, typename Int32TreeType, typename BoolTreeType, typename MeshDataAdapter>
2791inline void
2792expandNarrowband(
2793 TreeType& distTree,
2794 Int32TreeType& indexTree,
2795 BoolTreeType& maskTree,
2796 std::vector<typename BoolTreeType::LeafNodeType*>& maskNodes,
2797 const MeshDataAdapter& mesh,
2798 typename TreeType::ValueType exteriorBandWidth,
2799 typename TreeType::ValueType interiorBandWidth,
2800 typename TreeType::ValueType voxelSize)
2801{
2802 ExpandNarrowband<TreeType, MeshDataAdapter> expandOp(maskNodes, maskTree,
2803 distTree, indexTree, mesh, exteriorBandWidth, interiorBandWidth, voxelSize);
2804
2805 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, maskNodes.size()), expandOp);
2806
2807 tbb::parallel_for(tbb::blocked_range<size_t>(0, expandOp.updatedIndexNodes().size()),
2808 UnionValueMasks<typename TreeType::LeafNodeType, typename Int32TreeType::LeafNodeType>(
2809 expandOp.updatedDistNodes(), expandOp.updatedIndexNodes()));
2810
2811 tbb::task_group tasks;
2812 tasks.run(AddNodes<TreeType>(distTree, expandOp.newDistNodes()));
2813 tasks.run(AddNodes<Int32TreeType>(indexTree, expandOp.newIndexNodes()));
2814 tasks.wait();
2815
2816 maskTree.clear();
2817 maskTree.merge(expandOp.newMaskTree());
2818}
2819
2820
2821////////////////////////////////////////
2822
2823
2824// Transform values (sqrt, world space scaling and sign flip if sdf)
2825template<typename TreeType>
2826struct TransformValues
2827{
2828 using LeafNodeType = typename TreeType::LeafNodeType;
2829 using ValueType = typename TreeType::ValueType;
2830
2831 TransformValues(std::vector<LeafNodeType*>& nodes,
2832 ValueType voxelSize, bool unsignedDist)
2833 : mNodes(&nodes[0])
2834 , mVoxelSize(voxelSize)
2835 , mUnsigned(unsignedDist)
2836 {
2837 }
2838
2839 void operator()(const tbb::blocked_range<size_t>& range) const {
2840
2841 typename LeafNodeType::ValueOnIter iter;
2842
2843 const bool udf = mUnsigned;
2844 const ValueType w[2] = { -mVoxelSize, mVoxelSize };
2845
2846 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2847
2848 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2849 ValueType& val = const_cast<ValueType&>(iter.getValue());
2850 val = w[udf || (val < ValueType(0.0))] * std::sqrt(std::abs(val));
2851 }
2852 }
2853 }
2854
2855private:
2856 LeafNodeType * * const mNodes;
2857 const ValueType mVoxelSize;
2858 const bool mUnsigned;
2859};
2860
2861
2862// Inactivate values outside the (exBandWidth, inBandWidth) range.
2863template<typename TreeType>
2864struct InactivateValues
2865{
2866 using LeafNodeType = typename TreeType::LeafNodeType;
2867 using ValueType = typename TreeType::ValueType;
2868
2869 InactivateValues(std::vector<LeafNodeType*>& nodes,
2870 ValueType exBandWidth, ValueType inBandWidth)
2871 : mNodes(nodes.empty() ? nullptr : &nodes[0])
2872 , mExBandWidth(exBandWidth)
2873 , mInBandWidth(inBandWidth)
2874 {
2875 }
2876
2877 void operator()(const tbb::blocked_range<size_t>& range) const {
2878
2879 typename LeafNodeType::ValueOnIter iter;
2880 const ValueType exVal = mExBandWidth;
2881 const ValueType inVal = -mInBandWidth;
2882
2883 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2884
2885 for (iter = mNodes[n]->beginValueOn(); iter; ++iter) {
2886
2887 ValueType& val = const_cast<ValueType&>(iter.getValue());
2888
2889 const bool inside = val < ValueType(0.0);
2890
2891 if (inside && !(val > inVal)) {
2892 val = inVal;
2893 iter.setValueOff();
2894 } else if (!inside && !(val < exVal)) {
2895 val = exVal;
2896 iter.setValueOff();
2897 }
2898 }
2899 }
2900 }
2901
2902private:
2903 LeafNodeType * * const mNodes;
2904 const ValueType mExBandWidth, mInBandWidth;
2905};
2906
2907
2908template<typename TreeType>
2909struct OffsetValues
2910{
2911 using LeafNodeType = typename TreeType::LeafNodeType;
2912 using ValueType = typename TreeType::ValueType;
2913
2914 OffsetValues(std::vector<LeafNodeType*>& nodes, ValueType offset)
2915 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mOffset(offset)
2916 {
2917 }
2918
2919 void operator()(const tbb::blocked_range<size_t>& range) const {
2920
2921 const ValueType offset = mOffset;
2922
2923 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2924
2925 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
2926
2927 for (; iter; ++iter) {
2928 ValueType& val = const_cast<ValueType&>(iter.getValue());
2929 val += offset;
2930 }
2931 }
2932 }
2933
2934private:
2935 LeafNodeType * * const mNodes;
2936 const ValueType mOffset;
2937};
2938
2939
2940template<typename TreeType>
2941struct Renormalize
2942{
2943 using LeafNodeType = typename TreeType::LeafNodeType;
2944 using ValueType = typename TreeType::ValueType;
2945
2946 Renormalize(const TreeType& tree, const std::vector<LeafNodeType*>& nodes,
2947 ValueType* buffer, ValueType voxelSize)
2948 : mTree(&tree)
2949 , mNodes(nodes.empty() ? nullptr : &nodes[0])
2950 , mBuffer(buffer)
2951 , mVoxelSize(voxelSize)
2952 {
2953 }
2954
2955 void operator()(const tbb::blocked_range<size_t>& range) const
2956 {
2957 using Vec3Type = math::Vec3<ValueType>;
2958
2959 tree::ValueAccessor<const TreeType> acc(*mTree);
2960
2961 Coord ijk;
2962 Vec3Type up, down;
2963
2964 const ValueType dx = mVoxelSize, invDx = ValueType(1.0) / mVoxelSize;
2965
2966 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
2967
2968 ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
2969
2970 typename LeafNodeType::ValueOnCIter iter = mNodes[n]->cbeginValueOn();
2971 for (; iter; ++iter) {
2972
2973 const ValueType phi0 = *iter;
2974
2975 ijk = iter.getCoord();
2976
2977 up[0] = acc.getValue(ijk.offsetBy(1, 0, 0)) - phi0;
2978 up[1] = acc.getValue(ijk.offsetBy(0, 1, 0)) - phi0;
2979 up[2] = acc.getValue(ijk.offsetBy(0, 0, 1)) - phi0;
2980
2981 down[0] = phi0 - acc.getValue(ijk.offsetBy(-1, 0, 0));
2982 down[1] = phi0 - acc.getValue(ijk.offsetBy(0, -1, 0));
2983 down[2] = phi0 - acc.getValue(ijk.offsetBy(0, 0, -1));
2984
2985 const ValueType normSqGradPhi = math::GodunovsNormSqrd(phi0 > 0.0, down, up);
2986
2987 const ValueType diff = math::Sqrt(normSqGradPhi) * invDx - ValueType(1.0);
2988 const ValueType S = phi0 / (math::Sqrt(math::Pow2(phi0) + normSqGradPhi));
2989
2990 bufferData[iter.pos()] = phi0 - dx * S * diff;
2991 }
2992 }
2993 }
2994
2995private:
2996 TreeType const * const mTree;
2997 LeafNodeType const * const * const mNodes;
2998 ValueType * const mBuffer;
2999
3000 const ValueType mVoxelSize;
3001};
3002
3003
3004template<typename TreeType>
3005struct MinCombine
3006{
3007 using LeafNodeType = typename TreeType::LeafNodeType;
3008 using ValueType = typename TreeType::ValueType;
3009
3010 MinCombine(std::vector<LeafNodeType*>& nodes, const ValueType* buffer)
3011 : mNodes(nodes.empty() ? nullptr : &nodes[0]), mBuffer(buffer)
3012 {
3013 }
3014
3015 void operator()(const tbb::blocked_range<size_t>& range) const {
3016
3017 for (size_t n = range.begin(), N = range.end(); n < N; ++n) {
3018
3019 const ValueType* bufferData = &mBuffer[n * LeafNodeType::SIZE];
3020
3021 typename LeafNodeType::ValueOnIter iter = mNodes[n]->beginValueOn();
3022
3023 for (; iter; ++iter) {
3024 ValueType& val = const_cast<ValueType&>(iter.getValue());
3025 val = std::min(val, bufferData[iter.pos()]);
3026 }
3027 }
3028 }
3029
3030private:
3031 LeafNodeType * * const mNodes;
3032 ValueType const * const mBuffer;
3033};
3034
3035
3036} // mesh_to_volume_internal namespace
3037
3038/// @endcond
3039
3040
3041////////////////////////////////////////
3042
3043// Utility method implementation
3044
3045
3046template <typename FloatTreeT>
3047void
3049{
3050 using ConnectivityTable = mesh_to_volume_internal::LeafNodeConnectivityTable<FloatTreeT>;
3051
3052 // Build a node connectivity table where each leaf node has an offset into a
3053 // linearized list of nodes, and each leaf stores its six axis aligned neighbor
3054 // offsets
3055 ConnectivityTable nodeConnectivity(tree);
3056
3057 std::vector<size_t> zStartNodes, yStartNodes, xStartNodes;
3058
3059 // Store all nodes which do not have negative neighbors i.e. the nodes furthest
3060 // in -X, -Y, -Z. We sweep from lowest coordinate positions +axis and then
3061 // from the furthest positive coordinate positions -axis
3062 for (size_t n = 0; n < nodeConnectivity.size(); ++n) {
3063 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevX()[n]) {
3064 xStartNodes.push_back(n);
3065 }
3066
3067 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevY()[n]) {
3068 yStartNodes.push_back(n);
3069 }
3070
3071 if (ConnectivityTable::INVALID_OFFSET == nodeConnectivity.offsetsPrevZ()[n]) {
3072 zStartNodes.push_back(n);
3073 }
3074 }
3075
3076 using SweepingOp = mesh_to_volume_internal::SweepExteriorSign<FloatTreeT>;
3077
3078 // Sweep the exterior value signs (make them negative) up until the voxel intersection
3079 // with the isosurface. Do this in both lowest -> + and largest -> - directions
3080
3081 tbb::parallel_for(tbb::blocked_range<size_t>(0, zStartNodes.size()),
3082 SweepingOp(SweepingOp::Z_AXIS, zStartNodes, nodeConnectivity));
3083
3084 tbb::parallel_for(tbb::blocked_range<size_t>(0, yStartNodes.size()),
3085 SweepingOp(SweepingOp::Y_AXIS, yStartNodes, nodeConnectivity));
3086
3087 tbb::parallel_for(tbb::blocked_range<size_t>(0, xStartNodes.size()),
3088 SweepingOp(SweepingOp::X_AXIS, xStartNodes, nodeConnectivity));
3089
3090 const size_t numLeafNodes = nodeConnectivity.size();
3091 const size_t numVoxels = numLeafNodes * FloatTreeT::LeafNodeType::SIZE;
3092
3093 std::unique_ptr<bool[]> changedNodeMaskA{new bool[numLeafNodes]};
3094 std::unique_ptr<bool[]> changedNodeMaskB{new bool[numLeafNodes]};
3095 std::unique_ptr<bool[]> changedVoxelMask{new bool[numVoxels]};
3096
3097 mesh_to_volume_internal::fillArray(changedNodeMaskA.get(), true, numLeafNodes);
3098 mesh_to_volume_internal::fillArray(changedNodeMaskB.get(), false, numLeafNodes);
3099 mesh_to_volume_internal::fillArray(changedVoxelMask.get(), false, numVoxels);
3100
3101 const tbb::blocked_range<size_t> nodeRange(0, numLeafNodes);
3102
3103 bool nodesUpdated = false;
3104 do {
3105 // Perform per leaf node localized propagation of signs by looping over
3106 // all voxels and checking to see if any of their neighbors (within the
3107 // same leaf) are negative
3108 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedFillExteriorSign<FloatTreeT>(
3109 nodeConnectivity.nodes(), changedNodeMaskA.get()));
3110
3111 // For each leaf, check its axis aligned neighbors and propagate any changes
3112 // which occurred previously (in SeedFillExteriorSign OR in SyncVoxelMask) to
3113 // the leaf faces. Note that this operation stores the propagated face results
3114 // in a separate buffer (changedVoxelMask) to avoid writing to nodes being read
3115 // from other threads. Additionally mark any leaf nodes which will absorb any
3116 // changes from its neighbors in changedNodeMaskB
3117 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SeedPoints<FloatTreeT>(
3118 nodeConnectivity, changedNodeMaskA.get(), changedNodeMaskB.get(),
3119 changedVoxelMask.get()));
3120
3121 // Only nodes where a value was influenced by an adjacent node need to be
3122 // processed on the next pass.
3123 changedNodeMaskA.swap(changedNodeMaskB);
3124
3125 nodesUpdated = false;
3126 for (size_t n = 0; n < numLeafNodes; ++n) {
3127 nodesUpdated |= changedNodeMaskA[n];
3128 if (nodesUpdated) break;
3129 }
3130
3131 // Use the voxel mask updates in ::SeedPoints to actually assign the new values
3132 // across leaf node faces
3133 if (nodesUpdated) {
3134 tbb::parallel_for(nodeRange, mesh_to_volume_internal::SyncVoxelMask<FloatTreeT>(
3135 nodeConnectivity.nodes(), changedNodeMaskA.get(), changedVoxelMask.get()));
3136 }
3137 } while (nodesUpdated);
3138
3139} // void traceExteriorBoundaries()
3140
3141
3142////////////////////////////////////////
3143
3144template <typename T, Index Log2Dim, typename InteriorTest>
3145void
3146floodFillLeafNode(tree::LeafNode<T,Log2Dim>& leafNode, const InteriorTest& interiorTest) {
3147
3148 // Floods fills a single leaf node.
3149 // Starts with all voxels in NOT_VISITED.
3150 // Final result is voxels in either POSITIVE, NEGATIVE, or NOT_ASSIGNED.
3151 // Voxels that were categorized as NEGATIVE are negated.
3152 // The NOT_ASSIGNED is all voxels within 0.75 of the zero-crossing.
3153 //
3154 // NOT_VISITED voxels, if outside the 0.75 band, will query the oracle
3155 // to get a POSITIVE Or NEGATIVE sign (with interior being POSITIVE!)
3156 //
3157 // After setting a NOT_VISITED to either POSITIVE or NEGATIVE, an 8-way
3158 // depth-first floodfill is done, stopping at either the 0.75 boundary
3159 // or visited voxels.
3160 enum VoxelState {
3161 NOT_VISITED = 0,
3162 POSITIVE = 1,
3163 NEGATIVE = 2,
3164 NOT_ASSIGNED = 3
3165 };
3166
3167 const auto DIM = leafNode.DIM;
3168 const auto SIZE = leafNode.SIZE;
3169
3170 std::vector<VoxelState> voxelState(SIZE, NOT_VISITED);
3171
3172 std::vector<std::pair<Index, VoxelState>> offsetStack;
3173 offsetStack.reserve(SIZE);
3174
3175 for (Index offset=0; offset<SIZE; offset++) {
3176 const auto value = leafNode.getValue(offset);
3177
3178 // We do not assign anything for voxel close to the mesh
3179 // This condition is aligned with the condition in traceVoxelLine
3180 if (std::abs(value) <= 0.75) {
3181 voxelState[offset] = NOT_ASSIGNED;
3182 } else if (voxelState[offset] == NOT_VISITED) {
3183
3184 auto coord = leafNode.offsetToGlobalCoord(offset);
3185
3186 if (interiorTest(coord)){
3187 // Yes we assigne positive values to interior points
3188 // this is aligned with how meshToVolume works internally
3189 offsetStack.push_back({offset, POSITIVE});
3190 voxelState[offset] = POSITIVE;
3191 } else {
3192 offsetStack.push_back({offset, NEGATIVE});
3193 voxelState[offset] = NEGATIVE;
3194 }
3195
3196 while(!offsetStack.empty()){
3197
3198 auto [off, state] = offsetStack[offsetStack.size()-1];
3199 offsetStack.pop_back();
3200
3201 if (state == NEGATIVE) {
3202 leafNode.setValueOnly(off, -leafNode.getValue(off));
3203 }
3204
3205 // iterate over all neighbours and assign identical state
3206 // if they have not been visited and if they are far away
3207 // from the mesh (the condition is same as in traceVoxelLine)
3208 for (int dim=2; dim>=0; dim--){
3209 for (int i = -1; i <=1; ++(++i)){
3210 int dimIdx = (off >> dim * Log2Dim) % DIM;
3211 auto neighOff = off + (1 << dim * Log2Dim) * i;
3212 if ((0 < dimIdx) &&
3213 (dimIdx < (int)DIM - 1) &&
3214 (voxelState[neighOff] == NOT_VISITED)) {
3215
3216 if (std::abs(leafNode.getValue(neighOff)) <= 0.75) {
3217 voxelState[neighOff] = NOT_ASSIGNED;
3218 } else {
3219 offsetStack.push_back({neighOff, state});
3220 voxelState[neighOff] = state;
3221 }
3222 }
3223 }
3224 }
3225 }
3226 }
3227 }
3228}
3229
3230////////////////////////////////////////
3231
3232/// @brief Sets the sign of voxel values of `tree` based on the `interiorTest`
3233///
3234/// Inside is set to positive and outside to negative. This is in reverse to the usual
3235/// level set convention, but `meshToVolume` uses the opposite convention at certain point.
3236///
3237/// InteriorTest has to be a function `Coord -> bool` which evaluates true
3238/// inside of the mesh and false outside.
3239///
3240/// Furthermore, InteriorTest does not have to be thread-safe, but it has to be
3241/// copy constructible and evaluating different coppy has to be thread-safe.
3242///
3243/// Example of a interior test
3244///
3245/// auto acc = tree->getAccessor();
3246///
3247/// auto test = [acc = grid.getConstAccessor()](const Cood& coord) -> bool {
3248/// return acc->get(coord) <= 0 ? true : false;
3249/// }
3250template <typename FloatTreeT, typename InteriorTest>
3251void
3252evaluateInteriorTest(FloatTreeT& tree, InteriorTest interiorTest, InteriorTestStrategy interiorTestStrategy)
3253{
3254 static_assert(std::is_invocable_r<bool, InteriorTest, Coord>::value,
3255 "InteriorTest has to be a function `Coord -> bool`!");
3256 static_assert(std::is_copy_constructible_v<InteriorTest>,
3257 "InteriorTest has to be copyable!");
3258
3259 using LeafT = typename FloatTreeT::LeafNodeType;
3260
3261 if (interiorTestStrategy == EVAL_EVERY_VOXEL) {
3262
3263 auto op = [interiorTest](auto& node) {
3264 using Node = std::decay_t<decltype(node)>;
3265
3266 if constexpr (std::is_same_v<Node, LeafT>) {
3267
3268 for (auto iter = node.beginValueAll(); iter; ++iter) {
3269 if (!interiorTest(iter.getCoord())) {
3270 iter.setValue(-*iter);
3271 }
3272 }
3273
3274 } else {
3275 for (auto iter = node.beginChildOff(); iter; ++iter) {
3276 if (!interiorTest(iter.getCoord())) {
3277 iter.setValue(-*iter);
3278 }
3279 }
3280 }
3281 };
3282
3283 openvdb::tree::NodeManager nodes(tree);
3284 nodes.foreachBottomUp(op);
3285 }
3286
3287 if (interiorTestStrategy == EVAL_EVERY_TILE) {
3288
3289 auto op = [interiorTest](auto& node) {
3290 using Node = std::decay_t<decltype(node)>;
3291
3292 if constexpr (std::is_same_v<Node, LeafT>) {
3293 // // leaf node
3294 LeafT& leaf = static_cast<LeafT&>(node);
3295
3296 floodFillLeafNode(leaf, interiorTest);
3297
3298 } else {
3299 for (auto iter = node.beginChildOff(); iter; ++iter) {
3300 if (!interiorTest(iter.getCoord())) {
3301 iter.setValue(-*iter);
3302 }
3303 }
3304 }
3305 };
3306
3307 openvdb::tree::NodeManager nodes(tree);
3308 nodes.foreachBottomUp(op);
3309 }
3310} // void evaluateInteriorTest()
3311
3312////////////////////////////////////////
3313
3314
3315template <typename GridType, typename MeshDataAdapter, typename Interrupter, typename InteriorTest>
3316typename GridType::Ptr
3318 Interrupter& interrupter,
3319 const MeshDataAdapter& mesh,
3320 const math::Transform& transform,
3321 float exteriorBandWidth,
3322 float interiorBandWidth,
3323 int flags,
3324 typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid,
3325 InteriorTest interiorTest,
3326 InteriorTestStrategy interiorTestStrat)
3327{
3328 using GridTypePtr = typename GridType::Ptr;
3329 using TreeType = typename GridType::TreeType;
3330 using LeafNodeType = typename TreeType::LeafNodeType;
3331 using ValueType = typename GridType::ValueType;
3332
3333 using Int32GridType = typename GridType::template ValueConverter<Int32>::Type;
3334 using Int32TreeType = typename Int32GridType::TreeType;
3335
3336 using BoolTreeType = typename TreeType::template ValueConverter<bool>::Type;
3337
3338 //////////
3339
3340 // Setup
3341
3342 GridTypePtr distGrid(new GridType(std::numeric_limits<ValueType>::max()));
3343 distGrid->setTransform(transform.copy());
3344
3345 ValueType exteriorWidth = ValueType(exteriorBandWidth);
3346 ValueType interiorWidth = ValueType(interiorBandWidth);
3347
3348 // Note: inf interior width is all right, this value makes the converter fill
3349 // interior regions with distance values.
3350 if (!std::isfinite(exteriorWidth) || std::isnan(interiorWidth)) {
3351 std::stringstream msg;
3352 msg << "Illegal narrow band width: exterior = " << exteriorWidth
3353 << ", interior = " << interiorWidth;
3354 OPENVDB_LOG_DEBUG(msg.str());
3355 return distGrid;
3356 }
3357
3358 const ValueType voxelSize = ValueType(transform.voxelSize()[0]);
3359
3360 if (!std::isfinite(voxelSize) || math::isZero(voxelSize)) {
3361 std::stringstream msg;
3362 msg << "Illegal transform, voxel size = " << voxelSize;
3363 OPENVDB_LOG_DEBUG(msg.str());
3364 return distGrid;
3365 }
3366
3367 // Convert narrow band width from voxel units to world space units.
3368 exteriorWidth *= voxelSize;
3369 // Avoid the unit conversion if the interior band width is set to
3370 // inf or std::numeric_limits<float>::max().
3371 if (interiorWidth < std::numeric_limits<ValueType>::max()) {
3372 interiorWidth *= voxelSize;
3373 }
3374
3375 const bool computeSignedDistanceField = (flags & UNSIGNED_DISTANCE_FIELD) == 0;
3376 const bool removeIntersectingVoxels = (flags & DISABLE_INTERSECTING_VOXEL_REMOVAL) == 0;
3377 const bool renormalizeValues = (flags & DISABLE_RENORMALIZATION) == 0;
3378 const bool trimNarrowBand = (flags & DISABLE_NARROW_BAND_TRIMMING) == 0;
3379
3380 Int32GridType* indexGrid = nullptr;
3381
3382 typename Int32GridType::Ptr temporaryIndexGrid;
3383
3384 if (polygonIndexGrid) {
3385 indexGrid = polygonIndexGrid;
3386 } else {
3387 temporaryIndexGrid.reset(new Int32GridType(Int32(util::INVALID_IDX)));
3388 indexGrid = temporaryIndexGrid.get();
3389 }
3390
3391 indexGrid->newTree();
3392 indexGrid->setTransform(transform.copy());
3393
3394 if (computeSignedDistanceField) {
3395 distGrid->setGridClass(GRID_LEVEL_SET);
3396 } else {
3397 distGrid->setGridClass(GRID_UNKNOWN);
3398 interiorWidth = ValueType(0.0);
3399 }
3400
3401 TreeType& distTree = distGrid->tree();
3402 Int32TreeType& indexTree = indexGrid->tree();
3403
3404
3405 //////////
3406
3407 // Voxelize mesh
3408
3409 {
3410 using VoxelizationDataType = mesh_to_volume_internal::VoxelizationData<TreeType>;
3411 using DataTable = tbb::enumerable_thread_specific<typename VoxelizationDataType::Ptr>;
3412
3413 DataTable data;
3414 using Voxelizer =
3415 mesh_to_volume_internal::VoxelizePolygons<TreeType, MeshDataAdapter, Interrupter>;
3416
3417 const tbb::blocked_range<size_t> polygonRange(0, mesh.polygonCount());
3418
3419 tbb::parallel_for(polygonRange, Voxelizer(data, mesh, &interrupter));
3420
3421 for (typename DataTable::iterator i = data.begin(); i != data.end(); ++i) {
3422 VoxelizationDataType& dataItem = **i;
3423 mesh_to_volume_internal::combineData(
3424 distTree, indexTree, dataItem.distTree, dataItem.indexTree);
3425 }
3426 }
3427
3428 // The progress estimates are based on the observed average time for a few different
3429 // test cases and is only intended to provide some rough progression feedback to the user.
3430 if (interrupter.wasInterrupted(30)) return distGrid;
3431
3432
3433 //////////
3434
3435 // Classify interior and exterior regions
3436
3437 if (computeSignedDistanceField) {
3438
3439 /// If interior test is not provided
3440 if constexpr (std::is_same_v<InteriorTest, std::nullptr_t>) {
3441 // Determines the inside/outside state for the narrow band of voxels.
3442 (void) interiorTest; // Trigger usage.
3443 traceExteriorBoundaries(distTree);
3444 } else {
3445 evaluateInteriorTest(distTree, interiorTest, interiorTestStrat);
3446 }
3447
3448 /// Do not fix intersecting voxels if we have evaluated interior test for every voxel.
3449 bool signInitializedForEveryVoxel =
3450 /// interior test was provided i.e. not null
3451 !std::is_same_v<InteriorTest, std::nullptr_t> &&
3452 /// interior test was evaluated for every voxel
3453 interiorTestStrat == EVAL_EVERY_VOXEL;
3454
3455 if (!signInitializedForEveryVoxel) {
3456
3457 std::vector<LeafNodeType*> nodes;
3458 nodes.reserve(distTree.leafCount());
3459 distTree.getNodes(nodes);
3460
3461 const tbb::blocked_range<size_t> nodeRange(0, nodes.size());
3462
3463 using SignOp =
3464 mesh_to_volume_internal::ComputeIntersectingVoxelSign<TreeType, MeshDataAdapter>;
3465
3466 tbb::parallel_for(nodeRange, SignOp(nodes, distTree, indexTree, mesh));
3467
3468 if (interrupter.wasInterrupted(45)) return distGrid;
3469
3470 // Remove voxels created by self intersecting portions of the mesh.
3471 if (removeIntersectingVoxels) {
3472
3473 tbb::parallel_for(nodeRange,
3474 mesh_to_volume_internal::ValidateIntersectingVoxels<TreeType>(distTree, nodes));
3475
3476 tbb::parallel_for(nodeRange,
3477 mesh_to_volume_internal::RemoveSelfIntersectingSurface<TreeType>(
3478 nodes, distTree, indexTree));
3479
3480 tools::pruneInactive(distTree, /*threading=*/true);
3481 tools::pruneInactive(indexTree, /*threading=*/true);
3482 }
3483 }
3484 }
3485
3486 if (interrupter.wasInterrupted(50)) return distGrid;
3487
3488 if (distTree.activeVoxelCount() == 0) {
3489 distTree.clear();
3490 distTree.root().setBackground(exteriorWidth, /*updateChildNodes=*/false);
3491 return distGrid;
3492 }
3493
3494 // Transform values (world space scaling etc.).
3495 {
3496 std::vector<LeafNodeType*> nodes;
3497 nodes.reserve(distTree.leafCount());
3498 distTree.getNodes(nodes);
3499
3500 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3501 mesh_to_volume_internal::TransformValues<TreeType>(
3502 nodes, voxelSize, !computeSignedDistanceField));
3503 }
3504
3505 // Propagate sign information into tile regions.
3506 if (computeSignedDistanceField) {
3507 distTree.root().setBackground(exteriorWidth, /*updateChildNodes=*/false);
3508 tools::signedFloodFillWithValues(distTree, exteriorWidth, -interiorWidth);
3509 } else {
3510 tools::changeBackground(distTree, exteriorWidth);
3511 }
3512
3513 if (interrupter.wasInterrupted(54)) return distGrid;
3514
3515
3516 //////////
3517
3518 // Expand the narrow band region
3519
3520 const ValueType minBandWidth = voxelSize * ValueType(2.0);
3521
3522 if (interiorWidth > minBandWidth || exteriorWidth > minBandWidth) {
3523
3524 // Create the initial voxel mask.
3525 BoolTreeType maskTree(false);
3526
3527 {
3528 std::vector<LeafNodeType*> nodes;
3529 nodes.reserve(distTree.leafCount());
3530 distTree.getNodes(nodes);
3531
3532 mesh_to_volume_internal::ConstructVoxelMask<TreeType> op(maskTree, distTree, nodes);
3533 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, nodes.size()), op);
3534 }
3535
3536 // Progress estimation
3537 unsigned maxIterations = std::numeric_limits<unsigned>::max();
3538
3539 float progress = 54.0f, step = 0.0f;
3540 double estimated =
3541 2.0 * std::ceil((std::max(interiorWidth, exteriorWidth) - minBandWidth) / voxelSize);
3542
3543 if (estimated < double(maxIterations)) {
3544 maxIterations = unsigned(estimated);
3545 step = 40.0f / float(maxIterations);
3546 }
3547
3548 std::vector<typename BoolTreeType::LeafNodeType*> maskNodes;
3549
3550 unsigned count = 0;
3551 while (true) {
3552
3553 if (interrupter.wasInterrupted(int(progress))) return distGrid;
3554
3555 const size_t maskNodeCount = maskTree.leafCount();
3556 if (maskNodeCount == 0) break;
3557
3558 maskNodes.clear();
3559 maskNodes.reserve(maskNodeCount);
3560 maskTree.getNodes(maskNodes);
3561
3562 const tbb::blocked_range<size_t> range(0, maskNodes.size());
3563
3564 tbb::parallel_for(range,
3565 mesh_to_volume_internal::DiffLeafNodeMask<TreeType>(distTree, maskNodes));
3566
3567 mesh_to_volume_internal::expandNarrowband(distTree, indexTree, maskTree, maskNodes,
3568 mesh, exteriorWidth, interiorWidth, voxelSize);
3569
3570 if ((++count) >= maxIterations) break;
3571 progress += step;
3572 }
3573 }
3574
3575 if (interrupter.wasInterrupted(94)) return distGrid;
3576
3577 if (!polygonIndexGrid) indexGrid->clear();
3578
3579
3580 /////////
3581
3582 // Renormalize distances to smooth out bumps caused by self intersecting
3583 // and overlapping portions of the mesh and renormalize the level set.
3584
3585 if (computeSignedDistanceField && renormalizeValues) {
3586
3587 std::vector<LeafNodeType*> nodes;
3588 nodes.reserve(distTree.leafCount());
3589 distTree.getNodes(nodes);
3590
3591 std::unique_ptr<ValueType[]> buffer{new ValueType[LeafNodeType::SIZE * nodes.size()]};
3592
3593 const ValueType offset = ValueType(0.8 * voxelSize);
3594
3595 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3596 mesh_to_volume_internal::OffsetValues<TreeType>(nodes, -offset));
3597
3598 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3599 mesh_to_volume_internal::Renormalize<TreeType>(
3600 distTree, nodes, buffer.get(), voxelSize));
3601
3602 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3603 mesh_to_volume_internal::MinCombine<TreeType>(nodes, buffer.get()));
3604
3605 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3606 mesh_to_volume_internal::OffsetValues<TreeType>(
3607 nodes, offset - mesh_to_volume_internal::Tolerance<ValueType>::epsilon()));
3608 }
3609
3610 if (interrupter.wasInterrupted(99)) return distGrid;
3611
3612
3613 /////////
3614
3615 // Remove active voxels that exceed the narrow band limits
3616
3617 if (trimNarrowBand && std::min(interiorWidth, exteriorWidth) < voxelSize * ValueType(4.0)) {
3618
3619 std::vector<LeafNodeType*> nodes;
3620 nodes.reserve(distTree.leafCount());
3621 distTree.getNodes(nodes);
3622
3623 tbb::parallel_for(tbb::blocked_range<size_t>(0, nodes.size()),
3624 mesh_to_volume_internal::InactivateValues<TreeType>(
3625 nodes, exteriorWidth, computeSignedDistanceField ? interiorWidth : exteriorWidth));
3626
3627 tools::pruneLevelSet(
3628 distTree, exteriorWidth, computeSignedDistanceField ? -interiorWidth : -exteriorWidth);
3629 }
3630
3631 return distGrid;
3632}
3633
3634
3635template <typename GridType, typename MeshDataAdapter, typename InteriorTest>
3636typename GridType::Ptr
3638 const MeshDataAdapter& mesh,
3639 const math::Transform& transform,
3640 float exteriorBandWidth,
3641 float interiorBandWidth,
3642 int flags,
3643 typename GridType::template ValueConverter<Int32>::Type * polygonIndexGrid,
3644 InteriorTest /*interiorTest*/,
3645 InteriorTestStrategy /*interiorTestStrat*/)
3646{
3647 util::NullInterrupter nullInterrupter;
3648 return meshToVolume<GridType>(nullInterrupter, mesh, transform,
3649 exteriorBandWidth, interiorBandWidth, flags, polygonIndexGrid);
3650}
3651
3652
3653////////////////////////////////////////
3654
3655
3656//{
3657/// @cond OPENVDB_DOCS_INTERNAL
3658
3659/// @internal This overload is enabled only for grids with a scalar, floating-point ValueType.
3660template<typename GridType, typename Interrupter>
3661inline typename std::enable_if<std::is_floating_point<typename GridType::ValueType>::value,
3662 typename GridType::Ptr>::type
3663doMeshConversion(
3664 Interrupter& interrupter,
3665 const openvdb::math::Transform& xform,
3666 const std::vector<Vec3s>& points,
3667 const std::vector<Vec3I>& triangles,
3668 const std::vector<Vec4I>& quads,
3669 float exBandWidth,
3670 float inBandWidth,
3671 bool unsignedDistanceField = false)
3672{
3673 if (points.empty()) {
3674 return typename GridType::Ptr(new GridType(typename GridType::ValueType(exBandWidth)));
3675 }
3676
3677 const size_t numPoints = points.size();
3678 std::unique_ptr<Vec3s[]> indexSpacePoints{new Vec3s[numPoints]};
3679
3680 // transform points to local grid index space
3681 tbb::parallel_for(tbb::blocked_range<size_t>(0, numPoints),
3682 mesh_to_volume_internal::TransformPoints<Vec3s>(
3683 &points[0], indexSpacePoints.get(), xform));
3684
3685 const int conversionFlags = unsignedDistanceField ? UNSIGNED_DISTANCE_FIELD : 0;
3686
3687 if (quads.empty()) {
3688
3689 QuadAndTriangleDataAdapter<Vec3s, Vec3I>
3690 mesh(indexSpacePoints.get(), numPoints, &triangles[0], triangles.size());
3691
3692 return meshToVolume<GridType>(
3693 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3694
3695 } else if (triangles.empty()) {
3696
3697 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3698 mesh(indexSpacePoints.get(), numPoints, &quads[0], quads.size());
3699
3700 return meshToVolume<GridType>(
3701 interrupter, mesh, xform, exBandWidth, inBandWidth, conversionFlags);
3702 }
3703
3704 // pack primitives
3705
3706 const size_t numPrimitives = triangles.size() + quads.size();
3707 std::unique_ptr<Vec4I[]> prims{new Vec4I[numPrimitives]};
3708
3709 for (size_t n = 0, N = triangles.size(); n < N; ++n) {
3710 const Vec3I& triangle = triangles[n];
3711 Vec4I& prim = prims[n];
3712 prim[0] = triangle[0];
3713 prim[1] = triangle[1];
3714 prim[2] = triangle[2];
3715 prim[3] = util::INVALID_IDX;
3716 }
3717
3718 const size_t offset = triangles.size();
3719 for (size_t n = 0, N = quads.size(); n < N; ++n) {
3720 prims[offset + n] = quads[n];
3721 }
3722
3723 QuadAndTriangleDataAdapter<Vec3s, Vec4I>
3724 mesh(indexSpacePoints.get(), numPoints, prims.get(), numPrimitives);
3725
3726 return meshToVolume<GridType>(interrupter, mesh, xform,
3727 exBandWidth, inBandWidth, conversionFlags);
3728}
3729
3730
3731/// @internal This overload is enabled only for grids that do not have a scalar,
3732/// floating-point ValueType.
3733template<typename GridType, typename Interrupter>
3734inline typename std::enable_if<!std::is_floating_point<typename GridType::ValueType>::value,
3735 typename GridType::Ptr>::type
3736doMeshConversion(
3737 Interrupter&,
3738 const math::Transform& /*xform*/,
3739 const std::vector<Vec3s>& /*points*/,
3740 const std::vector<Vec3I>& /*triangles*/,
3741 const std::vector<Vec4I>& /*quads*/,
3742 float /*exBandWidth*/,
3743 float /*inBandWidth*/,
3744 bool /*unsignedDistanceField*/ = false)
3745{
3746 OPENVDB_THROW(TypeError,
3747 "mesh to volume conversion is supported only for scalar floating-point grids");
3748}
3749
3750/// @endcond
3751//}
3752
3753
3754////////////////////////////////////////
3755
3756
3757template<typename GridType>
3758typename GridType::Ptr
3760 const openvdb::math::Transform& xform,
3761 const std::vector<Vec3s>& points,
3762 const std::vector<Vec3I>& triangles,
3763 float halfWidth)
3764{
3765 util::NullInterrupter nullInterrupter;
3766 return meshToLevelSet<GridType>(nullInterrupter, xform, points, triangles, halfWidth);
3767}
3768
3769
3770template<typename GridType, typename Interrupter>
3771typename GridType::Ptr
3773 Interrupter& interrupter,
3774 const openvdb::math::Transform& xform,
3775 const std::vector<Vec3s>& points,
3776 const std::vector<Vec3I>& triangles,
3777 float halfWidth)
3778{
3779 std::vector<Vec4I> quads(0);
3780 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3781 halfWidth, halfWidth);
3782}
3783
3784
3785template<typename GridType>
3786typename GridType::Ptr
3788 const openvdb::math::Transform& xform,
3789 const std::vector<Vec3s>& points,
3790 const std::vector<Vec4I>& quads,
3791 float halfWidth)
3792{
3793 util::NullInterrupter nullInterrupter;
3794 return meshToLevelSet<GridType>(nullInterrupter, xform, points, quads, halfWidth);
3795}
3796
3797
3798template<typename GridType, typename Interrupter>
3799typename GridType::Ptr
3801 Interrupter& interrupter,
3802 const openvdb::math::Transform& xform,
3803 const std::vector<Vec3s>& points,
3804 const std::vector<Vec4I>& quads,
3805 float halfWidth)
3806{
3807 std::vector<Vec3I> triangles(0);
3808 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3809 halfWidth, halfWidth);
3810}
3811
3812
3813template<typename GridType>
3814typename GridType::Ptr
3816 const openvdb::math::Transform& xform,
3817 const std::vector<Vec3s>& points,
3818 const std::vector<Vec3I>& triangles,
3819 const std::vector<Vec4I>& quads,
3820 float halfWidth)
3821{
3822 util::NullInterrupter nullInterrupter;
3823 return meshToLevelSet<GridType>(
3824 nullInterrupter, xform, points, triangles, quads, halfWidth);
3825}
3826
3827
3828template<typename GridType, typename Interrupter>
3829typename GridType::Ptr
3831 Interrupter& interrupter,
3832 const openvdb::math::Transform& xform,
3833 const std::vector<Vec3s>& points,
3834 const std::vector<Vec3I>& triangles,
3835 const std::vector<Vec4I>& quads,
3836 float halfWidth)
3837{
3838 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3839 halfWidth, halfWidth);
3840}
3841
3842
3843template<typename GridType>
3844typename GridType::Ptr
3846 const openvdb::math::Transform& xform,
3847 const std::vector<Vec3s>& points,
3848 const std::vector<Vec3I>& triangles,
3849 const std::vector<Vec4I>& quads,
3850 float exBandWidth,
3851 float inBandWidth)
3852{
3853 util::NullInterrupter nullInterrupter;
3854 return meshToSignedDistanceField<GridType>(
3855 nullInterrupter, xform, points, triangles, quads, exBandWidth, inBandWidth);
3856}
3857
3858
3859template<typename GridType, typename Interrupter>
3860typename GridType::Ptr
3862 Interrupter& interrupter,
3863 const openvdb::math::Transform& xform,
3864 const std::vector<Vec3s>& points,
3865 const std::vector<Vec3I>& triangles,
3866 const std::vector<Vec4I>& quads,
3867 float exBandWidth,
3868 float inBandWidth)
3869{
3870 return doMeshConversion<GridType>(interrupter, xform, points, triangles,
3871 quads, exBandWidth, inBandWidth);
3872}
3873
3874
3875template<typename GridType>
3876typename GridType::Ptr
3878 const openvdb::math::Transform& xform,
3879 const std::vector<Vec3s>& points,
3880 const std::vector<Vec3I>& triangles,
3881 const std::vector<Vec4I>& quads,
3882 float bandWidth)
3883{
3884 util::NullInterrupter nullInterrupter;
3885 return meshToUnsignedDistanceField<GridType>(
3886 nullInterrupter, xform, points, triangles, quads, bandWidth);
3887}
3888
3889
3890template<typename GridType, typename Interrupter>
3891typename GridType::Ptr
3893 Interrupter& interrupter,
3894 const openvdb::math::Transform& xform,
3895 const std::vector<Vec3s>& points,
3896 const std::vector<Vec3I>& triangles,
3897 const std::vector<Vec4I>& quads,
3898 float bandWidth)
3899{
3900 return doMeshConversion<GridType>(interrupter, xform, points, triangles, quads,
3901 bandWidth, bandWidth, true);
3902}
3903
3904
3905////////////////////////////////////////////////////////////////////////////////
3906
3907
3908// Required by several of the tree nodes
3909inline std::ostream&
3910operator<<(std::ostream& ostr, const MeshToVoxelEdgeData::EdgeData& rhs)
3911{
3912 ostr << "{[ " << rhs.mXPrim << ", " << rhs.mXDist << "]";
3913 ostr << " [ " << rhs.mYPrim << ", " << rhs.mYDist << "]";
3914 ostr << " [ " << rhs.mZPrim << ", " << rhs.mZDist << "]}";
3915 return ostr;
3916}
3917
3918// Required by math::Abs
3919inline MeshToVoxelEdgeData::EdgeData
3921{
3922 return x;
3923}
3924
3925
3926////////////////////////////////////////
3927
3928
3930{
3931public:
3932
3934 const std::vector<Vec3s>& pointList,
3935 const std::vector<Vec4I>& polygonList);
3936
3937 void run(bool threaded = true);
3938
3939 GenEdgeData(GenEdgeData& rhs, tbb::split);
3940 inline void operator() (const tbb::blocked_range<size_t> &range);
3941 inline void join(GenEdgeData& rhs);
3942
3943 inline TreeType& tree() { return mTree; }
3944
3945private:
3946 void operator=(const GenEdgeData&) {}
3947
3948 struct Primitive { Vec3d a, b, c, d; Int32 index; };
3949
3950 template<bool IsQuad>
3951 inline void voxelize(const Primitive&);
3952
3953 template<bool IsQuad>
3954 inline bool evalPrimitive(const Coord&, const Primitive&);
3955
3956 inline bool rayTriangleIntersection( const Vec3d& origin, const Vec3d& dir,
3957 const Vec3d& a, const Vec3d& b, const Vec3d& c, double& t);
3958
3959
3960 TreeType mTree;
3961 Accessor mAccessor;
3962
3963 const std::vector<Vec3s>& mPointList;
3964 const std::vector<Vec4I>& mPolygonList;
3965
3966 // Used internally for acceleration
3967 using IntTreeT = TreeType::ValueConverter<Int32>::Type;
3968 IntTreeT mLastPrimTree;
3969 tree::ValueAccessor<IntTreeT> mLastPrimAccessor;
3970}; // class MeshToVoxelEdgeData::GenEdgeData
3971
3972
3973inline
3974MeshToVoxelEdgeData::GenEdgeData::GenEdgeData(
3975 const std::vector<Vec3s>& pointList,
3976 const std::vector<Vec4I>& polygonList)
3977 : mTree(EdgeData())
3978 , mAccessor(mTree)
3979 , mPointList(pointList)
3980 , mPolygonList(polygonList)
3981 , mLastPrimTree(Int32(util::INVALID_IDX))
3982 , mLastPrimAccessor(mLastPrimTree)
3983{
3984}
3985
3986
3987inline
3989 : mTree(EdgeData())
3990 , mAccessor(mTree)
3991 , mPointList(rhs.mPointList)
3992 , mPolygonList(rhs.mPolygonList)
3993 , mLastPrimTree(Int32(util::INVALID_IDX))
3994 , mLastPrimAccessor(mLastPrimTree)
3995{
3996}
3997
3998
3999inline void
4001{
4002 if (threaded) {
4003 tbb::parallel_reduce(tbb::blocked_range<size_t>(0, mPolygonList.size()), *this);
4004 } else {
4005 (*this)(tbb::blocked_range<size_t>(0, mPolygonList.size()));
4006 }
4007}
4008
4009
4010inline void
4012{
4013 using RootNodeType = TreeType::RootNodeType;
4014 using NodeChainType = RootNodeType::NodeChainType;
4015 static_assert(NodeChainType::Size > 1, "expected tree height > 1");
4016 using InternalNodeType = typename NodeChainType::template Get<1>;
4017
4018 Coord ijk;
4019 Index offset;
4020
4021 rhs.mTree.clearAllAccessors();
4022
4023 TreeType::LeafIter leafIt = rhs.mTree.beginLeaf();
4024 for ( ; leafIt; ++leafIt) {
4025 ijk = leafIt->origin();
4026
4027 TreeType::LeafNodeType* lhsLeafPt = mTree.probeLeaf(ijk);
4028
4029 if (!lhsLeafPt) {
4030
4031 mAccessor.addLeaf(rhs.mAccessor.probeLeaf(ijk));
4032 InternalNodeType* node = rhs.mAccessor.getNode<InternalNodeType>();
4033 node->stealNode<TreeType::LeafNodeType>(ijk, EdgeData(), false);
4034 rhs.mAccessor.clear();
4035
4036 } else {
4037
4038 TreeType::LeafNodeType::ValueOnCIter it = leafIt->cbeginValueOn();
4039 for ( ; it; ++it) {
4040
4041 offset = it.pos();
4042 const EdgeData& rhsValue = it.getValue();
4043
4044 if (!lhsLeafPt->isValueOn(offset)) {
4045 lhsLeafPt->setValueOn(offset, rhsValue);
4046 } else {
4047
4048 EdgeData& lhsValue = const_cast<EdgeData&>(lhsLeafPt->getValue(offset));
4049
4050 if (rhsValue.mXDist < lhsValue.mXDist) {
4051 lhsValue.mXDist = rhsValue.mXDist;
4052 lhsValue.mXPrim = rhsValue.mXPrim;
4053 }
4054
4055 if (rhsValue.mYDist < lhsValue.mYDist) {
4056 lhsValue.mYDist = rhsValue.mYDist;
4057 lhsValue.mYPrim = rhsValue.mYPrim;
4058 }
4059
4060 if (rhsValue.mZDist < lhsValue.mZDist) {
4061 lhsValue.mZDist = rhsValue.mZDist;
4062 lhsValue.mZPrim = rhsValue.mZPrim;
4063 }
4064
4065 }
4066 } // end value iteration
4067 }
4068 } // end leaf iteration
4069}
4070
4071
4072inline void
4073MeshToVoxelEdgeData::GenEdgeData::operator()(const tbb::blocked_range<size_t> &range)
4074{
4075 Primitive prim;
4076
4077 for (size_t n = range.begin(); n < range.end(); ++n) {
4078
4079 const Vec4I& verts = mPolygonList[n];
4080
4081 prim.index = Int32(n);
4082 prim.a = Vec3d(mPointList[verts[0]]);
4083 prim.b = Vec3d(mPointList[verts[1]]);
4084 prim.c = Vec3d(mPointList[verts[2]]);
4085
4086 if (util::INVALID_IDX != verts[3]) {
4087 prim.d = Vec3d(mPointList[verts[3]]);
4088 voxelize<true>(prim);
4089 } else {
4090 voxelize<false>(prim);
4091 }
4092 }
4093}
4094
4095
4096template<bool IsQuad>
4097inline void
4098MeshToVoxelEdgeData::GenEdgeData::voxelize(const Primitive& prim)
4099{
4100 std::deque<Coord> coordList;
4101 Coord ijk, nijk;
4102
4103 ijk = Coord::floor(prim.a);
4104 coordList.push_back(ijk);
4105
4106 evalPrimitive<IsQuad>(ijk, prim);
4107
4108 while (!coordList.empty()) {
4109
4110 ijk = coordList.back();
4111 coordList.pop_back();
4112
4113 for (Int32 i = 0; i < 26; ++i) {
4114 nijk = ijk + util::COORD_OFFSETS[i];
4115
4116 if (prim.index != mLastPrimAccessor.getValue(nijk)) {
4117 mLastPrimAccessor.setValue(nijk, prim.index);
4118 if(evalPrimitive<IsQuad>(nijk, prim)) coordList.push_back(nijk);
4119 }
4120 }
4121 }
4122}
4123
4124
4125template<bool IsQuad>
4126inline bool
4127MeshToVoxelEdgeData::GenEdgeData::evalPrimitive(const Coord& ijk, const Primitive& prim)
4128{
4129 Vec3d uvw, org(ijk[0], ijk[1], ijk[2]);
4130 bool intersecting = false;
4131 double t;
4132
4133 EdgeData edgeData;
4134 mAccessor.probeValue(ijk, edgeData);
4135
4136 // Evaluate first triangle
4137 double dist = (org -
4138 closestPointOnTriangleToPoint(prim.a, prim.c, prim.b, org, uvw)).lengthSqr();
4139
4140 if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.c, prim.b, t)) {
4141 if (t < edgeData.mXDist) {
4142 edgeData.mXDist = float(t);
4143 edgeData.mXPrim = prim.index;
4144 intersecting = true;
4145 }
4146 }
4147
4148 if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.c, prim.b, t)) {
4149 if (t < edgeData.mYDist) {
4150 edgeData.mYDist = float(t);
4151 edgeData.mYPrim = prim.index;
4152 intersecting = true;
4153 }
4154 }
4155
4156 if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.c, prim.b, t)) {
4157 if (t < edgeData.mZDist) {
4158 edgeData.mZDist = float(t);
4159 edgeData.mZPrim = prim.index;
4160 intersecting = true;
4161 }
4162 }
4163
4164 if (IsQuad) {
4165 // Split quad into a second triangle and calculate distance.
4166 double secondDist = (org -
4167 closestPointOnTriangleToPoint(prim.a, prim.d, prim.c, org, uvw)).lengthSqr();
4168
4169 if (secondDist < dist) dist = secondDist;
4170
4171 if (rayTriangleIntersection(org, Vec3d(1.0, 0.0, 0.0), prim.a, prim.d, prim.c, t)) {
4172 if (t < edgeData.mXDist) {
4173 edgeData.mXDist = float(t);
4174 edgeData.mXPrim = prim.index;
4175 intersecting = true;
4176 }
4177 }
4178
4179 if (rayTriangleIntersection(org, Vec3d(0.0, 1.0, 0.0), prim.a, prim.d, prim.c, t)) {
4180 if (t < edgeData.mYDist) {
4181 edgeData.mYDist = float(t);
4182 edgeData.mYPrim = prim.index;
4183 intersecting = true;
4184 }
4185 }
4186
4187 if (rayTriangleIntersection(org, Vec3d(0.0, 0.0, 1.0), prim.a, prim.d, prim.c, t)) {
4188 if (t < edgeData.mZDist) {
4189 edgeData.mZDist = float(t);
4190 edgeData.mZPrim = prim.index;
4191 intersecting = true;
4192 }
4193 }
4194 }
4195
4196 if (intersecting) mAccessor.setValue(ijk, edgeData);
4197
4198 return (dist < 0.86602540378443861);
4199}
4200
4201
4202inline bool
4203MeshToVoxelEdgeData::GenEdgeData::rayTriangleIntersection(
4204 const Vec3d& origin, const Vec3d& dir,
4205 const Vec3d& a, const Vec3d& b, const Vec3d& c,
4206 double& t)
4207{
4208 // Check if ray is parallel with triangle
4209
4210 Vec3d e1 = b - a;
4211 Vec3d e2 = c - a;
4212 Vec3d s1 = dir.cross(e2);
4213
4214 double divisor = s1.dot(e1);
4215 if (!(std::abs(divisor) > 0.0)) return false;
4216
4217 // Compute barycentric coordinates
4218
4219 double inv_divisor = 1.0 / divisor;
4220 Vec3d d = origin - a;
4221 double b1 = d.dot(s1) * inv_divisor;
4222
4223 if (b1 < 0.0 || b1 > 1.0) return false;
4224
4225 Vec3d s2 = d.cross(e1);
4226 double b2 = dir.dot(s2) * inv_divisor;
4227
4228 if (b2 < 0.0 || (b1 + b2) > 1.0) return false;
4229
4230 // Compute distance to intersection point
4231
4232 t = e2.dot(s2) * inv_divisor;
4233 return (t < 0.0) ? false : true;
4234}
4235
4236
4237////////////////////////////////////////
4238
4239
4240inline
4245
4246
4247inline void
4249 const std::vector<Vec3s>& pointList,
4250 const std::vector<Vec4I>& polygonList)
4251{
4252 GenEdgeData converter(pointList, polygonList);
4253 converter.run();
4254
4255 mTree.clear();
4256 mTree.merge(converter.tree());
4257}
4258
4259
4260inline void
4262 Accessor& acc,
4263 const Coord& ijk,
4264 std::vector<Vec3d>& points,
4265 std::vector<Index32>& primitives)
4266{
4267 EdgeData data;
4268 Vec3d point;
4269
4270 Coord coord = ijk;
4271
4272 if (acc.probeValue(coord, data)) {
4273
4274 if (data.mXPrim != util::INVALID_IDX) {
4275 point[0] = double(coord[0]) + data.mXDist;
4276 point[1] = double(coord[1]);
4277 point[2] = double(coord[2]);
4278
4279 points.push_back(point);
4280 primitives.push_back(data.mXPrim);
4281 }
4282
4283 if (data.mYPrim != util::INVALID_IDX) {
4284 point[0] = double(coord[0]);
4285 point[1] = double(coord[1]) + data.mYDist;
4286 point[2] = double(coord[2]);
4287
4288 points.push_back(point);
4289 primitives.push_back(data.mYPrim);
4290 }
4291
4292 if (data.mZPrim != util::INVALID_IDX) {
4293 point[0] = double(coord[0]);
4294 point[1] = double(coord[1]);
4295 point[2] = double(coord[2]) + data.mZDist;
4296
4297 points.push_back(point);
4298 primitives.push_back(data.mZPrim);
4299 }
4300
4301 }
4302
4303 coord[0] += 1;
4304
4305 if (acc.probeValue(coord, data)) {
4306
4307 if (data.mYPrim != util::INVALID_IDX) {
4308 point[0] = double(coord[0]);
4309 point[1] = double(coord[1]) + data.mYDist;
4310 point[2] = double(coord[2]);
4311
4312 points.push_back(point);
4313 primitives.push_back(data.mYPrim);
4314 }
4315
4316 if (data.mZPrim != util::INVALID_IDX) {
4317 point[0] = double(coord[0]);
4318 point[1] = double(coord[1]);
4319 point[2] = double(coord[2]) + data.mZDist;
4320
4321 points.push_back(point);
4322 primitives.push_back(data.mZPrim);
4323 }
4324 }
4325
4326 coord[2] += 1;
4327
4328 if (acc.probeValue(coord, data)) {
4329 if (data.mYPrim != util::INVALID_IDX) {
4330 point[0] = double(coord[0]);
4331 point[1] = double(coord[1]) + data.mYDist;
4332 point[2] = double(coord[2]);
4333
4334 points.push_back(point);
4335 primitives.push_back(data.mYPrim);
4336 }
4337 }
4338
4339 coord[0] -= 1;
4340
4341 if (acc.probeValue(coord, data)) {
4342
4343 if (data.mXPrim != util::INVALID_IDX) {
4344 point[0] = double(coord[0]) + data.mXDist;
4345 point[1] = double(coord[1]);
4346 point[2] = double(coord[2]);
4347
4348 points.push_back(point);
4349 primitives.push_back(data.mXPrim);
4350 }
4351
4352 if (data.mYPrim != util::INVALID_IDX) {
4353 point[0] = double(coord[0]);
4354 point[1] = double(coord[1]) + data.mYDist;
4355 point[2] = double(coord[2]);
4356
4357 points.push_back(point);
4358 primitives.push_back(data.mYPrim);
4359 }
4360 }
4361
4362
4363 coord[1] += 1;
4364
4365 if (acc.probeValue(coord, data)) {
4366
4367 if (data.mXPrim != util::INVALID_IDX) {
4368 point[0] = double(coord[0]) + data.mXDist;
4369 point[1] = double(coord[1]);
4370 point[2] = double(coord[2]);
4371
4372 points.push_back(point);
4373 primitives.push_back(data.mXPrim);
4374 }
4375 }
4376
4377 coord[2] -= 1;
4378
4379 if (acc.probeValue(coord, data)) {
4380
4381 if (data.mXPrim != util::INVALID_IDX) {
4382 point[0] = double(coord[0]) + data.mXDist;
4383 point[1] = double(coord[1]);
4384 point[2] = double(coord[2]);
4385
4386 points.push_back(point);
4387 primitives.push_back(data.mXPrim);
4388 }
4389
4390 if (data.mZPrim != util::INVALID_IDX) {
4391 point[0] = double(coord[0]);
4392 point[1] = double(coord[1]);
4393 point[2] = double(coord[2]) + data.mZDist;
4394
4395 points.push_back(point);
4396 primitives.push_back(data.mZPrim);
4397 }
4398 }
4399
4400 coord[0] += 1;
4401
4402 if (acc.probeValue(coord, data)) {
4403
4404 if (data.mZPrim != util::INVALID_IDX) {
4405 point[0] = double(coord[0]);
4406 point[1] = double(coord[1]);
4407 point[2] = double(coord[2]) + data.mZDist;
4408
4409 points.push_back(point);
4410 primitives.push_back(data.mZPrim);
4411 }
4412 }
4413}
4414
4415
4416template<typename GridType, typename VecType>
4417typename GridType::Ptr
4419 const openvdb::math::Transform& xform,
4420 typename VecType::ValueType halfWidth)
4421{
4422 const Vec3s pmin = Vec3s(xform.worldToIndex(bbox.min()));
4423 const Vec3s pmax = Vec3s(xform.worldToIndex(bbox.max()));
4424
4425 Vec3s points[8];
4426 points[0] = Vec3s(pmin[0], pmin[1], pmin[2]);
4427 points[1] = Vec3s(pmin[0], pmin[1], pmax[2]);
4428 points[2] = Vec3s(pmax[0], pmin[1], pmax[2]);
4429 points[3] = Vec3s(pmax[0], pmin[1], pmin[2]);
4430 points[4] = Vec3s(pmin[0], pmax[1], pmin[2]);
4431 points[5] = Vec3s(pmin[0], pmax[1], pmax[2]);
4432 points[6] = Vec3s(pmax[0], pmax[1], pmax[2]);
4433 points[7] = Vec3s(pmax[0], pmax[1], pmin[2]);
4434
4435 Vec4I faces[6];
4436 faces[0] = Vec4I(0, 1, 2, 3); // bottom
4437 faces[1] = Vec4I(7, 6, 5, 4); // top
4438 faces[2] = Vec4I(4, 5, 1, 0); // front
4439 faces[3] = Vec4I(6, 7, 3, 2); // back
4440 faces[4] = Vec4I(0, 3, 7, 4); // left
4441 faces[5] = Vec4I(1, 5, 6, 2); // right
4442
4443 QuadAndTriangleDataAdapter<Vec3s, Vec4I> mesh(points, 8, faces, 6);
4444
4445 return meshToVolume<GridType>(mesh, xform, static_cast<float>(halfWidth), static_cast<float>(halfWidth));
4446}
4447
4448
4449////////////////////////////////////////
4450
4451
4452// Explicit Template Instantiation
4453
4454#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
4455
4456#ifdef OPENVDB_INSTANTIATE_MESHTOVOLUME
4458#endif
4459
4460#define _FUNCTION(TreeT) \
4461 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \
4462 const QuadAndTriangleDataAdapter<Vec3s, Vec3I>&, const openvdb::math::Transform&, \
4463 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy)
4465#undef _FUNCTION
4466
4467#define _FUNCTION(TreeT) \
4468 Grid<TreeT>::Ptr meshToVolume<Grid<TreeT>>(util::NullInterrupter&, \
4469 const QuadAndTriangleDataAdapter<Vec3s, Vec4I>&, const openvdb::math::Transform&, \
4470 float, float, int, Grid<TreeT>::ValueConverter<Int32>::Type*, std::nullptr_t, InteriorTestStrategy)
4472#undef _FUNCTION
4473
4474#define _FUNCTION(TreeT) \
4475 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4476 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec3I>&, \
4477 float)
4479#undef _FUNCTION
4480
4481#define _FUNCTION(TreeT) \
4482 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4483 const openvdb::math::Transform&, const std::vector<Vec3s>&, const std::vector<Vec4I>&, \
4484 float)
4486#undef _FUNCTION
4487
4488#define _FUNCTION(TreeT) \
4489 Grid<TreeT>::Ptr meshToLevelSet<Grid<TreeT>>(util::NullInterrupter&, \
4490 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4491 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float)
4493#undef _FUNCTION
4494
4495#define _FUNCTION(TreeT) \
4496 Grid<TreeT>::Ptr meshToSignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \
4497 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4498 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float, float)
4500#undef _FUNCTION
4501
4502#define _FUNCTION(TreeT) \
4503 Grid<TreeT>::Ptr meshToUnsignedDistanceField<Grid<TreeT>>(util::NullInterrupter&, \
4504 const openvdb::math::Transform&, const std::vector<Vec3s>&, \
4505 const std::vector<Vec3I>&, const std::vector<Vec4I>&, float)
4507#undef _FUNCTION
4508
4509#define _FUNCTION(TreeT) \
4510 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3s>&, \
4511 const openvdb::math::Transform&, float)
4513#undef _FUNCTION
4514
4515#define _FUNCTION(TreeT) \
4516 Grid<TreeT>::Ptr createLevelSetBox<Grid<TreeT>>(const math::BBox<Vec3d>&, \
4517 const openvdb::math::Transform&, double)
4519#undef _FUNCTION
4520
4521#define _FUNCTION(TreeT) \
4522 void traceExteriorBoundaries(TreeT&)
4524#undef _FUNCTION
4525
4526#endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
4527
4528
4529} // namespace tools
4530} // namespace OPENVDB_VERSION_NAME
4531} // namespace openvdb
4532
4533#endif // OPENVDB_TOOLS_MESH_TO_VOLUME_HAS_BEEN_INCLUDED
Efficient multi-threaded replacement of the background values in tree.
OPENVDB_API std::ostream & operator<<(std::ostream &os, half h)
Output h to os, formatted as a float.
Defined various multi-threaded utility functions for trees.
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Axis-aligned bounding box.
Definition BBox.h:24
const Vec3T & max() const
Return a const reference to the maximum point of this bounding box.
Definition BBox.h:64
const Vec3T & min() const
Return a const reference to the minimum point of this bounding box.
Definition BBox.h:62
Signed (x, y, z) 32-bit integer coordinates.
Definition Coord.h:25
static Coord floor(const Vec3< T > &xyz)
Return the largest integer coordinates that are not greater than xyz (node centered conversion).
Definition Coord.h:56
Definition Transform.h:40
Vec3d voxelSize() const
Return the size of a voxel using the linear component of the map.
Definition Transform.h:93
Ptr copy() const
Definition Transform.h:50
T dot(const Vec3< T > &v) const
Dot product.
Definition Vec3.h:191
Vec3< T > cross(const Vec3< T > &v) const
Return the cross product of "this" vector and v;.
Definition Vec3.h:220
Definition Vec4.h:25
TreeType & tree()
Definition MeshToVolume.h:3943
void operator()(const tbb::blocked_range< size_t > &range)
Definition MeshToVolume.h:4073
void run(bool threaded=true)
Definition MeshToVolume.h:4000
void join(GenEdgeData &rhs)
Definition MeshToVolume.h:4011
GenEdgeData(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Definition MeshToVolume.h:3974
Extracts and stores voxel edge intersection data from a mesh.
Definition MeshToVolume.h:456
void convert(const std::vector< Vec3s > &pointList, const std::vector< Vec4I > &polygonList)
Threaded method to extract voxel edge data, the closest intersection point and corresponding primitiv...
Definition MeshToVolume.h:4248
MeshToVoxelEdgeData()
Definition MeshToVolume.h:4241
Accessor getAccessor()
Definition MeshToVolume.h:517
void getEdgeData(Accessor &acc, const Coord &ijk, std::vector< Vec3d > &points, std::vector< Index32 > &primitives)
Returns intersection points with corresponding primitive indices for the given ijk voxel.
Definition MeshToVolume.h:4261
Base class for tree-traversal iterators over all leaf nodes (but not leaf voxels)
Definition TreeIterator.h:1187
Templated block class to hold specific data types and a fixed number of values determined by Log2Dim....
Definition LeafNode.h:38
const ValueType & getValue(const Coord &xyz) const
Return the value of the voxel at the given coordinates.
Definition LeafNode.h:1045
void setValueOnly(const Coord &xyz, const ValueType &val)
Set the value of the voxel at the given coordinates but don't change its active state.
Definition LeafNode.h:1103
static const Index DIM
Definition LeafNode.h:50
Coord offsetToGlobalCoord(Index n) const
Return the global coordinates for a linear table offset.
Definition LeafNode.h:1034
static const Index SIZE
Definition LeafNode.h:53
Base class for tree-traversal iterators over tile and voxel values.
Definition TreeIterator.h:617
const ValueT & getValue() const
Return the tile or voxel value to which this iterator is currently pointing.
Definition TreeIterator.h:692
Definition Tree.h:178
void clearAllAccessors()
Clear all registered accessors.
Definition Tree.h:1356
void merge(Tree &other, MergePolicy=MERGE_ACTIVE_STATES)
Efficiently merge another tree into this tree using one of several schemes.
Definition Tree.h:1669
_RootNodeType RootNodeType
Definition Tree.h:183
LeafNodeType * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains voxel (x, y, z). If no such node exists,...
Definition Tree.h:1543
void clear()
Remove all tiles from this tree and all nodes other than the root node.
Definition Tree.h:1297
typename RootNodeType::LeafNodeType LeafNodeType
Definition Tree.h:186
LeafIter beginLeaf()
Return an iterator over all leaf nodes in this tree.
Definition Tree.h:975
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition ValueAccessor.h:367
void clear() override final
Remove all the cached nodes and invalidate the corresponding hash-keys.
Definition ValueAccessor.h:880
LeafNodeT * probeLeaf(const Coord &xyz)
Return a pointer to the leaf node that contains the voxel coordinate xyz. If no LeafNode exists,...
Definition ValueAccessor.h:836
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
Definition ValueAccessor.h:492
NodeT * getNode()
Return the node of type NodeT that has been cached on this accessor. If this accessor does not cache ...
Definition ValueAccessor.h:848
Convert polygonal meshes that consist of quads and/or triangles into signed or unsigned distance fiel...
Partitions points into BucketLog2Dim aligned buckets using a parallel radix-based sorting algorithm.
#define OPENVDB_LOG_DEBUG(message)
In debug builds only, log a debugging message of the form 'someVar << "text" << .....
Definition logging.h:266
BBox< Coord > CoordBBox
Definition NanoVDB.h:2535
bool operator<(const Tuple< SIZE, T0 > &t0, const Tuple< SIZE, T1 > &t1)
Definition Tuple.h:174
Vec3< double > Vec3d
Definition Vec3.h:664
Axis
Definition Math.h:901
@ Z_AXIS
Definition Math.h:904
@ X_AXIS
Definition Math.h:902
@ Y_AXIS
Definition Math.h:903
Vec3< float > Vec3s
Definition Vec3.h:663
OPENVDB_API Vec3d closestPointOnTriangleToPoint(const Vec3d &a, const Vec3d &b, const Vec3d &c, const Vec3d &p, Vec3d &uvw)
Closest Point on Triangle to Point. Given a triangle abc and a point p, return the point on abc close...
const std::enable_if<!VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition Composite.h:110
MeshToVolumeFlags
Mesh to volume conversion flags.
Definition MeshToVolume.h:59
@ DISABLE_INTERSECTING_VOXEL_REMOVAL
Definition MeshToVolume.h:69
@ DISABLE_RENORMALIZATION
Definition MeshToVolume.h:73
@ DISABLE_NARROW_BAND_TRIMMING
Definition MeshToVolume.h:77
@ UNSIGNED_DISTANCE_FIELD
Definition MeshToVolume.h:65
void floodFillLeafNode(tree::LeafNode< T, Log2Dim > &leafNode, const InteriorTest &interiorTest)
Definition MeshToVolume.h:3146
GridType::Ptr meshToLevelSet(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, float halfWidth=float(LEVEL_SET_HALF_WIDTH))
Convert a triangle mesh to a level set volume.
Definition MeshToVolume.h:3759
void traceExteriorBoundaries(FloatTreeT &tree)
Traces the exterior voxel boundary of closed objects in the input volume tree. Exterior voxels are ma...
Definition MeshToVolume.h:3048
GridType::Ptr meshToUnsignedDistanceField(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float bandWidth)
Convert a triangle and quad mesh to an unsigned distance field.
Definition MeshToVolume.h:3877
GridType::Ptr meshToVolume(const MeshDataAdapter &mesh, const math::Transform &transform, float exteriorBandWidth=3.0f, float interiorBandWidth=3.0f, int flags=0, typename GridType::template ValueConverter< Int32 >::Type *polygonIndexGrid=nullptr, InteriorTest interiorTest=nullptr, InteriorTestStrategy interiorTestStrat=EVAL_EVERY_VOXEL)
Definition MeshToVolume.h:3637
GridType::Ptr createLevelSetBox(const math::BBox< VecType > &bbox, const openvdb::math::Transform &xform, typename VecType::ValueType halfWidth=LEVEL_SET_HALF_WIDTH)
Return a grid of type GridType containing a narrow-band level set representation of a box.
Definition MeshToVolume.h:4418
GridType::Ptr meshToSignedDistanceField(const openvdb::math::Transform &xform, const std::vector< Vec3s > &points, const std::vector< Vec3I > &triangles, const std::vector< Vec4I > &quads, float exBandWidth, float inBandWidth)
Convert a triangle and quad mesh to a signed distance field with an asymmetrical narrow band.
Definition MeshToVolume.h:3845
void evaluateInteriorTest(FloatTreeT &tree, InteriorTest interiorTest, InteriorTestStrategy interiorTestStrategy)
Sets the sign of voxel values of tree based on the interiorTest
Definition MeshToVolume.h:3252
InteriorTestStrategy
Different staregies how to determine sign of an SDF when using interior test.
Definition MeshToVolume.h:83
@ EVAL_EVERY_VOXEL
Definition MeshToVolume.h:87
@ EVAL_EVERY_TILE
Evaluates interior test at least once per tile and flood fills within the tile.
Definition MeshToVolume.h:90
constexpr Index32 INVALID_IDX
Definition Util.h:19
constexpr Coord COORD_OFFSETS[26]
coordinate offset table for neighboring voxels
Definition Util.h:22
bool wasInterrupted(T *i, int percent=-1)
Definition NullInterrupter.h:49
Index32 Index
Definition Types.h:54
math::Vec4< Index32 > Vec4I
Definition Types.h:88
@ GRID_LEVEL_SET
Definition Types.h:455
@ GRID_UNKNOWN
Definition Types.h:454
uint32_t Index32
Definition Types.h:52
math::Vec3< Index32 > Vec3I
Definition Types.h:73
int32_t Int32
Definition Types.h:56
Definition Exceptions.h:13
Definition Coord.h:589
#define OPENVDB_THROW(exception, message)
Definition Exceptions.h:74
Internal edge data type.
Definition MeshToVolume.h:462
EdgeData(float dist=1.0)
Definition MeshToVolume.h:463
bool operator==(const EdgeData &rhs) const
Definition MeshToVolume.h:481
Index32 mXPrim
Definition MeshToVolume.h:487
float mZDist
Definition MeshToVolume.h:486
EdgeData operator+(const T &) const
Definition MeshToVolume.h:476
float mYDist
Definition MeshToVolume.h:486
float mXDist
Definition MeshToVolume.h:486
Index32 mZPrim
Definition MeshToVolume.h:487
EdgeData operator-(const T &) const
Definition MeshToVolume.h:477
Index32 mYPrim
Definition MeshToVolume.h:487
EdgeData operator-() const
Definition MeshToVolume.h:478
Contiguous quad and triangle data adapter class.
Definition MeshToVolume.h:186
size_t polygonCount() const
Definition MeshToVolume.h:206
void getIndexSpacePoint(size_t n, size_t v, Vec3d &pos) const
Returns position pos in local grid index space for polygon n and vertex v.
Definition MeshToVolume.h:216
QuadAndTriangleDataAdapter(const std::vector< PointType > &points, const std::vector< PolygonType > &polygons)
Definition MeshToVolume.h:188
size_t pointCount() const
Definition MeshToVolume.h:207
QuadAndTriangleDataAdapter(const PointType *pointArray, size_t pointArraySize, const PolygonType *polygonArray, size_t polygonArraySize)
Definition MeshToVolume.h:197
size_t vertexCount(size_t n) const
Vertex count for polygon n.
Definition MeshToVolume.h:210
Base class for interrupters.
Definition NullInterrupter.h:26
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition version.h.in:121
#define OPENVDB_USE_VERSION_NAMESPACE
Definition version.h.in:212
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition version.h.in:157