OpenVDB  8.0.1
LevelSetSphere.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
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 <type_traits>
21 
22 #include <tbb/enumerable_thread_specific.h>
23 #include <tbb/parallel_for.h>
24 #include <tbb/parallel_reduce.h>
25 #include <tbb/blocked_range.h>
26 #include <thread>
27 
28 namespace openvdb {
30 namespace OPENVDB_VERSION_NAME {
31 namespace tools {
32 
47 template<typename GridType, typename InterruptT>
48 typename GridType::Ptr
49 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
50  float halfWidth = float(LEVEL_SET_HALF_WIDTH),
51  InterruptT* interrupt = nullptr, bool threaded = true);
52 
66 template<typename GridType>
67 typename GridType::Ptr
68 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
69  float halfWidth = float(LEVEL_SET_HALF_WIDTH), bool threaded = true)
70 {
71  return createLevelSetSphere<GridType, util::NullInterrupter>(radius,center,voxelSize,halfWidth,nullptr,threaded);
72 }
73 
74 
76 
77 
84 template<typename GridT, typename InterruptT = util::NullInterrupter>
86 {
87 public:
88  using TreeT = typename GridT::TreeType;
89  using ValueT = typename GridT::ValueType;
90  using Vec3T = typename math::Vec3<ValueT>;
91  static_assert(std::is_floating_point<ValueT>::value,
92  "level set grids must have scalar, floating-point value types");
93 
104  LevelSetSphere(ValueT radius, const Vec3T &center, InterruptT* interrupt = nullptr)
105  : mRadius(radius), mCenter(center), mInterrupt(interrupt)
106  {
107  if (mRadius<=0) OPENVDB_THROW(ValueError, "radius must be positive");
108  }
109 
115  typename GridT::Ptr getLevelSet(ValueT voxelSize, ValueT halfWidth, bool threaded = true)
116  {
117  mGrid = createLevelSet<GridT>(voxelSize, halfWidth);
118  this->rasterSphere(voxelSize, halfWidth, threaded);
119  mGrid->setGridClass(GRID_LEVEL_SET);
120  return mGrid;
121  }
122 
123 private:
124  void rasterSphere(ValueT dx, ValueT w, bool threaded)
125  {
126  if (!(dx>0.0f)) OPENVDB_THROW(ValueError, "voxel size must be positive");
127  if (!(w>1)) OPENVDB_THROW(ValueError, "half-width must be larger than one");
128 
129  // Define radius of sphere and narrow-band in voxel units
130  const ValueT r0 = mRadius/dx, rmax = r0 + w;
131 
132  // Radius below the Nyquist frequency
133  if (r0 < 1.5f) return;
134 
135  // Define center of sphere in voxel units
136  const Vec3T c(mCenter[0]/dx, mCenter[1]/dx, mCenter[2]/dx);
137 
138  // Define bounds of the voxel coordinates
139  const int imin=math::Floor(c[0]-rmax), imax=math::Ceil(c[0]+rmax);
140  const int jmin=math::Floor(c[1]-rmax), jmax=math::Ceil(c[1]+rmax);
141  const int kmin=math::Floor(c[2]-rmax), kmax=math::Ceil(c[2]+rmax);
142 
143  // Allocate a ValueAccessor for accelerated random access
144  typename GridT::Accessor accessor = mGrid->getAccessor();
145 
146  if (mInterrupt) mInterrupt->start("Generating level set of sphere");
147 
148  tbb::enumerable_thread_specific<TreeT> pool(mGrid->tree());
149 
150  auto kernel = [&](const tbb::blocked_range<int>& r) {
151  openvdb::Coord ijk;
152  int &i = ijk[0], &j = ijk[1], &k = ijk[2], m=1;
153  TreeT &tree = pool.local();
154  typename GridT::Accessor acc(tree);
155  // Compute signed distances to sphere using leapfrogging in k
156  for (i = r.begin(); i <= r.end(); ++i) {
157  if (util::wasInterrupted(mInterrupt)) return;
158  const auto x2 = math::Pow2(ValueT(i) - c[0]);
159  for (j = jmin; j <= jmax; ++j) {
160  const auto x2y2 = math::Pow2(ValueT(j) - c[1]) + x2;
161  for (k = kmin; k <= kmax; k += m) {
162  m = 1;
163  // Distance in voxel units to sphere
164  const auto v = math::Sqrt(x2y2 + math::Pow2(ValueT(k)-c[2]))-r0;
165  const auto d = math::Abs(v);
166  if (d < w) { // inside narrow band
167  acc.setValue(ijk, dx*v);// distance in world units
168  } else { // outside narrow band
169  m += math::Floor(d-w);// leapfrog
170  }
171  }//end leapfrog over k
172  }//end loop over j
173  }//end loop over i
174  };// kernel
175 
176  if (threaded) {
177  // The code blow is making use of a TLS container to minimize the number of concurrent trees
178  // initially populated by tbb::parallel_for and subsequently merged by tbb::parallel_reduce.
179  // Experiments have demonstrated this approach to outperform others, including serial reduction
180  // and a custom concurrent reduction implementation.
181  tbb::parallel_for(tbb::blocked_range<int>(imin, imax, 128), kernel);
182  using RangeT = tbb::blocked_range<typename tbb::enumerable_thread_specific<TreeT>::iterator>;
183  struct Op {
184  const bool mDelete;
185  TreeT *mTree;
186  Op(TreeT &tree) : mDelete(false), mTree(&tree) {}
187  Op(const Op& other, tbb::split) : mDelete(true), mTree(new TreeT(other.mTree->background())) {}
188  ~Op() { if (mDelete) delete mTree; }
189  void operator()(RangeT &r) { for (auto i=r.begin(); i!=r.end(); ++i) this->merge(*i);}
190  void join(Op &other) { this->merge(*(other.mTree)); }
191  void merge(TreeT &tree) { mTree->merge(tree, openvdb::MERGE_ACTIVE_STATES); }
192  } op( mGrid->tree() );
193  tbb::parallel_reduce(RangeT(pool.begin(), pool.end(), 4), op);
194  } else {
195  kernel(tbb::blocked_range<int>(imin, imax));//serial
196  mGrid->tree().merge(*pool.begin(), openvdb::MERGE_ACTIVE_STATES);
197  }
198 
199  // Define consistent signed distances outside the narrow-band
200  tools::signedFloodFill(mGrid->tree(), threaded);
201 
202  if (mInterrupt) mInterrupt->end();
203  }
204 
205  const ValueT mRadius;
206  const Vec3T mCenter;
207  InterruptT* mInterrupt;
208  typename GridT::Ptr mGrid;
209 };// LevelSetSphere
210 
211 
213 
214 
215 template<typename GridType, typename InterruptT>
216 typename GridType::Ptr
217 createLevelSetSphere(float radius, const openvdb::Vec3f& center, float voxelSize,
218  float halfWidth, InterruptT* interrupt, bool threaded)
219 {
220  // GridType::ValueType is required to be a floating-point scalar.
221  static_assert(std::is_floating_point<typename GridType::ValueType>::value,
222  "level set grids must have scalar, floating-point value types");
223 
224  using ValueT = typename GridType::ValueType;
225  LevelSetSphere<GridType, InterruptT> factory(ValueT(radius), center, interrupt);
226  return factory.getLevelSet(ValueT(voxelSize), ValueT(halfWidth), threaded);
227 }
228 
229 } // namespace tools
230 } // namespace OPENVDB_VERSION_NAME
231 } // namespace openvdb
232 
233 #endif // OPENVDB_TOOLS_LEVELSETSPHERE_HAS_BEEN_INCLUDED
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: openvdb/Exceptions.h:65
Signed (x, y, z) 32-bit integer coordinates.
Definition: Coord.h:26
Definition: Vec3.h:24
Generates a signed distance field (or narrow band level set) to a single sphere.
Definition: LevelSetSphere.h:86
LevelSetSphere(ValueT radius, const Vec3T &center, InterruptT *interrupt=nullptr)
Constructor.
Definition: LevelSetSphere.h:104
typename GridT::TreeType TreeT
Definition: LevelSetSphere.h:88
GridT::Ptr getLevelSet(ValueT voxelSize, ValueT halfWidth, bool threaded=true)
Definition: LevelSetSphere.h:115
typename math::Vec3< ValueT > Vec3T
Definition: LevelSetSphere.h:90
typename GridT::ValueType ValueT
Definition: LevelSetSphere.h:89
int Floor(float x)
Return the floor of x.
Definition: Math.h:851
int Ceil(float x)
Return the ceiling of x.
Definition: Math.h:859
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
float Sqrt(float x)
Return the square root of a floating-point value.
Definition: Math.h:764
Type Pow2(Type x)
Return x2.
Definition: Math.h:551
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:68
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:266
bool wasInterrupted(T *i, int percent=-1)
Definition: NullInterrupter.h:49
static const Real LEVEL_SET_HALF_WIDTH
Definition: openvdb/Types.h:321
@ GRID_LEVEL_SET
Definition: openvdb/Types.h:315
@ MERGE_ACTIVE_STATES
Definition: openvdb/Types.h:367
Definition: openvdb/Exceptions.h:13
#define OPENVDB_THROW(exception, message)
Definition: openvdb/Exceptions.h:74
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:101
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:153