OpenVDB 10.0.1
Loading...
Searching...
No Matches
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
31namespace openvdb {
33namespace OPENVDB_VERSION_NAME {
34namespace 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
50template<typename GridType, typename InterruptT>
51typename GridType::Ptr
52createLevelSetSphere(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
69template<typename GridType>
70typename GridType::Ptr
71createLevelSetSphere(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
87template<typename GridT, typename InterruptT = util::NullInterrupter>
89{
90public:
91 using TreeT = typename GridT::TreeType;
92 using ValueT = typename GridT::ValueType;
93 using Vec3T = typename math::Vec3<ValueT>;
94 static_assert(std::is_floating_point<ValueT>::value,
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
126private:
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
218template<typename GridType, typename InterruptT>
219typename GridType::Ptr
220createLevelSetSphere(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.
224 static_assert(std::is_floating_point<typename GridType::ValueType>::value,
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
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
GridType::Ptr createLevelSetSphere(float radius, const openvdb::Vec3f &center, float voxelSize, float halfWidth=float(LEVEL_SET_HALF_WIDTH), InterruptT *interrupt=nullptr, bool threaded=true)
Return a grid of type GridType containing a narrow-band level set representation of a sphere.
Definition LevelSetSphere.h:220
@ 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