OpenVDB  9.0.0
LevelSetSphere.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 LevelSetSphere.h
5 ///
6 /// @brief Generate a narrow-band level set of sphere.
7 ///
8 /// @note By definition a level set has a fixed narrow band width
9 /// (the half width is defined by LEVEL_SET_HALF_WIDTH in Types.h),
10 /// whereas an SDF can have a variable narrow band width.
11 
12 #ifndef OPENVDB_TOOLS_LEVELSETSPHERE_HAS_BEEN_INCLUDED
13 #define OPENVDB_TOOLS_LEVELSETSPHERE_HAS_BEEN_INCLUDED
14 
15 #include <openvdb/Grid.h>
16 #include <openvdb/Types.h>
17 #include <openvdb/math/Math.h>
19 #include "SignedFloodFill.h"
20 #include <openvdb/openvdb.h>
21 #include <type_traits>
22 
23 #include <tbb/enumerable_thread_specific.h>
24 #include <tbb/parallel_for.h>
25 #include <tbb/parallel_reduce.h>
26 #include <tbb/blocked_range.h>
27 #include <thread>
28 
29 namespace openvdb {
31 namespace OPENVDB_VERSION_NAME {
32 namespace tools {
33 
34 /// @brief Return a grid of type @c GridType containing a narrow-band level set
35 /// representation of a sphere.
36 ///
37 /// @param radius radius of the sphere in world units
38 /// @param center center of the sphere in world units
39 /// @param voxelSize voxel size in world units
40 /// @param halfWidth half the width of the narrow band, in voxel units
41 /// @param interrupt a pointer adhering to the util::NullInterrupter interface
42 /// @param threaded if true multi-threading is enabled (true by default)
43 ///
44 /// @note @c GridType::ValueType must be a floating-point scalar.
45 /// @note The leapfrog algorithm employed in this method is best suited
46 /// for a single large sphere. For multiple small spheres consider
47 /// using the faster algorithm in ParticlesToLevelSet.h
48 template<typename GridType, typename InterruptT>
49 typename GridType::Ptr
50 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
51  float halfWidth = float(LEVEL_SET_HALF_WIDTH),
52  InterruptT* interrupt = nullptr, bool threaded = true);
53 
54 /// @brief Return a grid of type @c GridType containing a narrow-band level set
55 /// representation of a sphere.
56 ///
57 /// @param radius radius of the sphere in world units
58 /// @param center center of the sphere in world units
59 /// @param voxelSize voxel size in world units
60 /// @param halfWidth half the width of the narrow band, in voxel units
61 /// @param threaded if true multi-threading is enabled (true by default)
62 ///
63 /// @note @c GridType::ValueType must be a floating-point scalar.
64 /// @note The leapfrog algorithm employed in this method is best suited
65 /// for a single large sphere. For multiple small spheres consider
66 /// using the faster algorithm in ParticlesToLevelSet.h
67 template<typename GridType>
68 typename GridType::Ptr
69 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
70  float halfWidth = float(LEVEL_SET_HALF_WIDTH), bool threaded = true)
71 {
72  return createLevelSetSphere<GridType, util::NullInterrupter>(radius,center,voxelSize,halfWidth,nullptr,threaded);
73 }
74 
75 
76 ////////////////////////////////////////
77 
78 
79 /// @brief Generates a signed distance field (or narrow band level
80 /// set) to a single sphere.
81 ///
82 /// @note The leapfrog algorithm employed in this class is best
83 /// suited for a single large sphere. For multiple small spheres consider
84 /// using the faster algorithm in tools/ParticlesToLevelSet.h
85 template<typename GridT, typename InterruptT = util::NullInterrupter>
87 {
88 public:
89  using TreeT = typename GridT::TreeType;
90  using ValueT = typename GridT::ValueType;
91  using Vec3T = typename math::Vec3<ValueT>;
93  "level set grids must have scalar, floating-point value types");
94 
95  /// @brief Constructor
96  ///
97  /// @param radius radius of the sphere in world units
98  /// @param center center of the sphere in world units
99  /// @param interrupt pointer to optional interrupter. Use template
100  /// argument util::NullInterrupter if no interruption is desired.
101  ///
102  /// @note If the radius of the sphere is smaller than
103  /// 1.5*voxelSize, i.e. the sphere is smaller than the Nyquist
104  /// frequency of the grid, it is ignored!
105  LevelSetSphere(ValueT radius, const Vec3T &center, InterruptT* interrupt = nullptr)
106  : mRadius(radius), mCenter(center), mInterrupt(interrupt)
107  {
108  if (mRadius<=0) OPENVDB_THROW(ValueError, "radius must be positive");
109  }
110 
111  /// @return a narrow-band level set of the sphere
112  ///
113  /// @param voxelSize Size of voxels in world units
114  /// @param halfWidth Half-width of narrow-band in voxel units
115  /// @param threaded If true multi-threading is enabled (true by default)
116  typename GridT::Ptr getLevelSet(ValueT voxelSize, ValueT halfWidth, bool threaded = true)
117  {
118  mGrid = createLevelSet<GridT>(voxelSize, halfWidth);
119  this->rasterSphere(voxelSize, halfWidth, threaded);
120  mGrid->setGridClass(GRID_LEVEL_SET);
121  return mGrid;
122  }
123 
124 private:
125  void rasterSphere(ValueT dx, ValueT w, bool threaded)
126  {
127  if (!(dx>0.0f)) OPENVDB_THROW(ValueError, "voxel size must be positive");
128  if (!(w>1)) OPENVDB_THROW(ValueError, "half-width must be larger than one");
129 
130  // Define radius of sphere and narrow-band in voxel units
131  const ValueT r0 = mRadius/dx, rmax = r0 + w;
132 
133  // Radius below the Nyquist frequency
134  if (r0 < 1.5f) return;
135 
136  // Define center of sphere in voxel units
137  const Vec3T c(mCenter[0]/dx, mCenter[1]/dx, mCenter[2]/dx);
138 
139  // Define bounds of the voxel coordinates
140  const int imin=math::Floor(c[0]-rmax), imax=math::Ceil(c[0]+rmax);
141  const int jmin=math::Floor(c[1]-rmax), jmax=math::Ceil(c[1]+rmax);
142  const int kmin=math::Floor(c[2]-rmax), kmax=math::Ceil(c[2]+rmax);
143 
144  // Allocate a ValueAccessor for accelerated random access
145  typename GridT::Accessor accessor = mGrid->getAccessor();
146 
147  if (mInterrupt) mInterrupt->start("Generating level set of sphere");
148 
149  tbb::enumerable_thread_specific<TreeT> pool(mGrid->tree());
150 
151  auto kernel = [&](const tbb::blocked_range<int>& r) {
152  openvdb::Coord ijk;
153  int &i = ijk[0], &j = ijk[1], &k = ijk[2], m=1;
154  TreeT &tree = pool.local();
155  typename GridT::Accessor acc(tree);
156  // Compute signed distances to sphere using leapfrogging in k
157  for (i = r.begin(); i <= r.end(); ++i) {
158  if (util::wasInterrupted(mInterrupt)) return;
159  const auto x2 = math::Pow2(ValueT(i) - c[0]);
160  for (j = jmin; j <= jmax; ++j) {
161  const auto x2y2 = math::Pow2(ValueT(j) - c[1]) + x2;
162  for (k = kmin; k <= kmax; k += m) {
163  m = 1;
164  // Distance in voxel units to sphere
165  const auto v = math::Sqrt(x2y2 + math::Pow2(ValueT(k)-c[2]))-r0;
166  const auto d = math::Abs(v);
167  if (d < w) { // inside narrow band
168  acc.setValue(ijk, dx*v);// distance in world units
169  } else { // outside narrow band
170  m += math::Floor(d-w);// leapfrog
171  }
172  }//end leapfrog over k
173  }//end loop over j
174  }//end loop over i
175  };// kernel
176 
177  if (threaded) {
178  // The code blow is making use of a TLS container to minimize the number of concurrent trees
179  // initially populated by tbb::parallel_for and subsequently merged by tbb::parallel_reduce.
180  // Experiments have demonstrated this approach to outperform others, including serial reduction
181  // and a custom concurrent reduction implementation.
182  tbb::parallel_for(tbb::blocked_range<int>(imin, imax, 128), kernel);
183  using RangeT = tbb::blocked_range<typename tbb::enumerable_thread_specific<TreeT>::iterator>;
184  struct Op {
185  const bool mDelete;
186  TreeT *mTree;
187  Op(TreeT &tree) : mDelete(false), mTree(&tree) {}
188  Op(const Op& other, tbb::split) : mDelete(true), mTree(new TreeT(other.mTree->background())) {}
189  ~Op() { if (mDelete) delete mTree; }
190  void operator()(RangeT &r) { for (auto i=r.begin(); i!=r.end(); ++i) this->merge(*i);}
191  void join(Op &other) { this->merge(*(other.mTree)); }
192  void merge(TreeT &tree) { mTree->merge(tree, openvdb::MERGE_ACTIVE_STATES); }
193  } op( mGrid->tree() );
194  tbb::parallel_reduce(RangeT(pool.begin(), pool.end(), 4), op);
195  } else {
196  kernel(tbb::blocked_range<int>(imin, imax));//serial
197  mGrid->tree().merge(*pool.begin(), openvdb::MERGE_ACTIVE_STATES);
198  }
199 
200  // Define consistent signed distances outside the narrow-band
201  tools::signedFloodFill(mGrid->tree(), threaded);
202 
203  if (mInterrupt) mInterrupt->end();
204  }
205 
206  const ValueT mRadius;
207  const Vec3T mCenter;
208  InterruptT* mInterrupt;
209  typename GridT::Ptr mGrid;
210 };// LevelSetSphere
211 
212 
213 ////////////////////////////////////////
214 
215 
216 template<typename GridType, typename InterruptT>
217 typename GridType::Ptr
218 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
219  float halfWidth, InterruptT* interrupt, bool threaded)
220 {
221  // GridType::ValueType is required to be a floating-point scalar.
223  "level set grids must have scalar, floating-point value types");
224 
225  using ValueT = typename GridType::ValueType;
226  LevelSetSphere<GridType, InterruptT> factory(ValueT(radius), center, interrupt);
227  return factory.getLevelSet(ValueT(voxelSize), ValueT(halfWidth), threaded);
228 }
229 
230 
231 ////////////////////////////////////////
232 
233 
234 // Explicit Template Instantiation
235 
236 #ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
237 
238 #ifdef OPENVDB_INSTANTIATE_LEVELSETSPHERE
240 #endif
241 
242 #define _FUNCTION(TreeT) \
243  Grid<TreeT>::Ptr createLevelSetSphere<Grid<TreeT>>(float, const openvdb::Vec3f&, float, float, \
244  util::NullInterrupter*, bool)
246 #undef _FUNCTION
247 
248 #endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
249 
250 
251 } // namespace tools
252 } // namespace OPENVDB_VERSION_NAME
253 } // namespace openvdb
254 
255 #endif // OPENVDB_TOOLS_LEVELSETSPHERE_HAS_BEEN_INCLUDED
ValueT value
Definition: GridBuilder.h:1287
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Propagate the signs of distance values from the active voxels in the narrow band to the inactive valu...
Definition: Exceptions.h:65
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:25
Definition: Vec3.h:24
Generates a signed distance field (or narrow band level set) to a single sphere.
Definition: LevelSetSphere.h:87
LevelSetSphere(ValueT radius, const Vec3T &center, InterruptT *interrupt=nullptr)
Constructor.
Definition: LevelSetSphere.h:105
typename GridT::TreeType TreeT
Definition: LevelSetSphere.h:89
GridT::Ptr getLevelSet(ValueT voxelSize, ValueT halfWidth, bool threaded=true)
Definition: LevelSetSphere.h:116
typename math::Vec3< ValueT > Vec3T
Definition: LevelSetSphere.h:91
typename GridT::ValueType ValueT
Definition: LevelSetSphere.h:90
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:859
Coord Abs(const Coord &xyz)
Definition: Coord.h:514
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
int Floor(float x)
Return the floor of x.
Definition: Math.h:851
GridType::Ptr createLevelSetSphere(float radius, const openvdb::Vec3f &center, float voxelSize, float halfWidth=float(LEVEL_SET_HALF_WIDTH), bool threaded=true)
Return a grid of type GridType containing a narrow-band level set representation of a sphere.
Definition: LevelSetSphere.h:69
void signedFloodFill(TreeOrLeafManagerT &tree, bool threaded=true, size_t grainSize=1, Index minLevel=0)
Set the values of all inactive voxels and tiles of a narrow-band level set from the signs of the acti...
Definition: SignedFloodFill.h:267
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
static const Real LEVEL_SET_HALF_WIDTH
Definition: Types.h:343
@ GRID_LEVEL_SET
Definition: Types.h:337
@ MERGE_ACTIVE_STATES
Definition: Types.h:389
Definition: Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202
#define OPENVDB_REAL_TREE_INSTANTIATE(Function)
Definition: version.h.in:147