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