91 using TreeT =
typename GridT::TreeType;
92 using ValueT =
typename GridT::ValueType;
94 static_assert(std::is_floating_point<ValueT>::value,
95 "level set grids must have scalar, floating-point value types");
108 : mRadius(radius), mCenter(center), mInterrupt(interrupt)
120 mGrid = createLevelSet<GridT>(voxelSize, halfWidth);
121 this->rasterSphere(voxelSize, halfWidth, threaded);
127 void rasterSphere(ValueT dx, ValueT w,
bool threaded)
133 const ValueT r0 = mRadius/dx, rmax = r0 + w;
136 if (r0 < 1.5f)
return;
139 const Vec3T c(mCenter[0]/dx, mCenter[1]/dx, mCenter[2]/dx);
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);
147 typename GridT::Accessor accessor = mGrid->getAccessor();
149 if (mInterrupt) mInterrupt->start(
"Generating level set of sphere");
151 tbb::enumerable_thread_specific<TreeT> pool(mGrid->tree());
153 auto kernel = [&](
const tbb::blocked_range<int>& r) {
155 int &i = ijk[0], &j = ijk[1], &k = ijk[2], m=1;
156 TreeT &tree = pool.local();
157 typename GridT::Accessor acc(tree);
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) {
167 const auto v = math::Sqrt(x2y2 + math::Pow2(ValueT(k)-c[2]))-r0;
168 const auto d = math::Abs(v);
170 acc.setValue(ijk, dx*v);
172 m += math::Floor(d-w);
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>;
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()(
const RangeT &r) {
for (
auto i=r.begin(); i!=r.end(); ++i) this->merge(*i);}
193 void join(Op &other) { this->merge(*(other.mTree)); }
195 } op( mGrid->tree() );
196 tbb::parallel_reduce(RangeT(pool.begin(), pool.end(), 4), op);
198 kernel(tbb::blocked_range<int>(imin, imax));
203 tools::signedFloodFill(mGrid->tree(), threaded);
205 if (mInterrupt) mInterrupt->end();
208 const ValueT mRadius;
210 InterruptT* mInterrupt;
211 typename GridT::Ptr mGrid;
221 float halfWidth, InterruptT* interrupt,
bool threaded)
224 static_assert(std::is_floating_point<typename GridType::ValueType>::value,
225 "level set grids must have scalar, floating-point value types");
227 using ValueT =
typename GridType::ValueType;
229 return factory.
getLevelSet(ValueT(voxelSize), ValueT(halfWidth), threaded);