42 #ifndef OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
43 #define OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED
57 template<
typename GridT,
58 typename InterruptT = util::NullInterrupter>
72 :
BaseType(grid, interrupt), mTask(0), mMask(NULL)
79 :
BaseType(other), mTask(other.mTask), mMask(other.mMask)
90 if (mTask) mTask(const_cast<LevelSetFilter*>(
this), range);
100 void laplacian(
const GridT* mask = NULL);
108 void gaussian(
int width = 1,
const GridT* mask = NULL);
113 void offset(ValueType offset,
const GridT* mask = NULL);
121 void median(
int width = 1,
const GridT* mask = NULL);
128 void mean(
int width = 1,
const GridT* mask = NULL);
131 typedef typename GridT::ConstAccessor AccT;
132 typedef typename TreeType::LeafNodeType LeafT;
133 typedef typename LeafT::ValueOnIter VoxelIterT;
134 typedef typename LeafT::ValueOnCIter VoxelCIterT;
136 typedef typename RangeType::Iterator LeafIterT;
139 typename boost::function<void (LevelSetFilter*, const RangeType&)> mTask;
145 const int n = BaseType::getGrainSize();
147 tbb::parallel_for(BaseType::leafs().leafRange(n), *
this);
149 (*this)(BaseType::leafs().leafRange());
151 if (swap) BaseType::leafs().swapLeafBuffer(1, n==0);
157 AffineCombine(
const GridT& grid,
const GridT* mask)
158 : mGrid(&grid), mMask(mask),
mAcc(mask ? new AccT(mask->tree()) : NULL) {}
159 ~AffineCombine() {
delete mAcc; }
160 inline ValueType alpha(
const Coord& xyz)
const
163 const Vec3R world = mGrid->indexToWorld(xyz);
164 const Vec3R voxel = mMask->worldToIndex(world);
165 return tools::BoxSampler::sample<AccT>(*
mAcc, voxel);
167 inline ValueType operator()(
const Coord& xyz, ValueType unfiltered, ValueType filtered)
const
169 if (mMask==NULL)
return filtered;
170 const ValueType alpha = this->alpha(xyz);
171 return alpha*filtered + (1-alpha)*unfiltered;
181 template <
size_t Axis>
183 Avg(
const GridT& grid,
Int32 w) : acc(grid.tree()), width(w), frac(1/ValueType(2*w+1)) {}
184 ValueType operator()(
Coord xyz) {
185 ValueType sum = zeroVal<ValueType>();
187 for (i -= width; i <= j; ++i) sum += acc.getValue(xyz);
190 typename GridT::ConstAccessor acc;
192 const ValueType frac;
196 template <
typename AvgT>
197 void doBox(
const RangeType& r,
Int32 w);
198 void doBoxX(
const RangeType& r,
Int32 w) { this->doBox<Avg<0> >(r,w); }
199 void doBoxZ(
const RangeType& r,
Int32 w) { this->doBox<Avg<1> >(r,w); }
200 void doBoxY(
const RangeType& r,
Int32 w) { this->doBox<Avg<2> >(r,w); }
201 void doMedian(
const RangeType&,
int);
202 void doMeanCurvature(
const RangeType&);
203 void doLaplacian(
const RangeType&);
204 void doOffset(
const RangeType&, ValueType);
211 template<
typename Gr
idT,
typename InterruptT>
217 BaseType::startInterrupter(
"Median-value flow of level set");
219 BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
221 mTask = boost::bind(&LevelSetFilter::doMedian, _1, _2,
std::max(1, width));
226 BaseType::endInterrupter();
230 template<
typename Gr
idT,
typename InterruptT>
236 BaseType::startInterrupter(
"Mean-value flow of level set");
240 BaseType::endInterrupter();
243 template<
typename Gr
idT,
typename InterruptT>
249 BaseType::startInterrupter(
"Gaussian flow of level set");
251 for (
int n=0; n<4; ++n) this->box(width);
253 BaseType::endInterrupter();
256 template<
typename Gr
idT,
typename InterruptT>
260 BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
264 mTask = boost::bind(&LevelSetFilter::doBoxX, _1, _2, width);
267 mTask = boost::bind(&LevelSetFilter::doBoxY, _1, _2, width);
270 mTask = boost::bind(&LevelSetFilter::doBoxZ, _1, _2, width);
277 template<
typename Gr
idT,
typename InterruptT>
283 BaseType::startInterrupter(
"Mean-curvature flow of level set");
285 BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
287 mTask = boost::bind(&LevelSetFilter::doMeanCurvature, _1, _2);
292 BaseType::endInterrupter();
295 template<
typename Gr
idT,
typename InterruptT>
301 BaseType::startInterrupter(
"Laplacian flow of level set");
303 BaseType::leafs().rebuildAuxBuffers(1, BaseType::getGrainSize()==0);
305 mTask = boost::bind(&LevelSetFilter::doLaplacian, _1, _2);
310 BaseType::endInterrupter();
314 template<
typename Gr
idT,
typename InterruptT>
320 BaseType::startInterrupter(
"Offsetting level set");
322 BaseType::leafs().removeAuxBuffers();
326 while (offset-dist >
ValueType(0.001)*CFL && BaseType::checkInterrupter()) {
330 mTask = boost::bind(&LevelSetFilter::doOffset, _1, _2, copysign(delta, value));
336 BaseType::endInterrupter();
343 template<
typename Gr
idT,
typename InterruptT>
347 BaseType::checkInterrupter();
349 const ValueType dx = BaseType::voxelSize(),dt =
math::Pow2(dx) / ValueType(3.0);
351 AffineCombine comb(BaseType::grid(), mMask);
352 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
353 BufferT& buffer = leafIter.buffer(1);
354 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
355 stencil.moveTo(iter);
356 ValueType A = stencil.getCenterValue(), B = A + dt * stencil.meanCurvatureNormGrad();
357 buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
369 template<
typename Gr
idT,
typename InterruptT>
371 LevelSetFilter<GridT, InterruptT>::doLaplacian(
const RangeType& range)
373 BaseType::checkInterrupter();
375 const ValueType dx = BaseType::voxelSize(), dt =
math::Pow2(dx) / ValueType(6.0);
376 math::GradStencil<GridType> stencil(BaseType::grid(), dx);
377 AffineCombine comb(BaseType::grid(), mMask);
378 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
379 BufferT& buffer = leafIter.buffer(1);
380 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
381 stencil.moveTo(iter);
382 const ValueType A = stencil.getCenterValue(), B = A + dt * stencil.laplacian();
383 buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
389 template<
typename Gr
idT,
typename InterruptT>
391 LevelSetFilter<GridT, InterruptT>::doOffset(
const RangeType& range, ValueType offset)
393 BaseType::checkInterrupter();
395 AffineCombine comb(BaseType::grid(), mMask);
396 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
397 for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
398 iter.setValue(*iter + comb.alpha(iter.getCoord())*offset);
402 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
403 for (VoxelIterT iter = leafIter->beginValueOn(); iter; ++iter) {
404 iter.setValue(*iter + offset);
411 template<
typename Gr
idT,
typename InterruptT>
413 LevelSetFilter<GridT, InterruptT>::doMedian(
const RangeType& range,
int width)
415 BaseType::checkInterrupter();
416 typename math::DenseStencil<GridType> stencil(BaseType::grid(), width);
417 AffineCombine comb(BaseType::grid(), mMask);
418 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
419 BufferT& buffer = leafIter.buffer(1);
420 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
421 stencil.moveTo(iter);
422 const ValueType A = stencil.getCenterValue(), B = stencil.median();
423 buffer.setValue(iter.pos(), comb(stencil.getCenterCoord(), A, B));
429 template<
typename Gr
idT,
typename InterruptT>
430 template <
typename AvgT>
432 LevelSetFilter<GridT, InterruptT>::doBox(
const RangeType& range,
Int32 w)
434 BaseType::checkInterrupter();
435 AvgT avg(BaseType::grid(), w);
436 AffineCombine comb(BaseType::grid(), mMask);
437 for (LeafIterT leafIter=range.begin(); leafIter; ++leafIter) {
438 BufferT& buffer = leafIter.buffer(1);
439 for (VoxelCIterT iter = leafIter->cbeginValueOn(); iter; ++iter) {
440 const Coord xyz = iter.getCoord();
441 buffer.setValue(iter.pos(), comb(xyz, *iter, avg(xyz)));
450 #endif // OPENVDB_TOOLS_LEVELSETFILTER_HAS_BEEN_INCLUDED