OpenVDB  2.0.0
PointAdvect.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 //
35 
36 #ifndef OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
37 #define OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
38 
39 #include <openvdb/openvdb.h>
40 #include <openvdb/math/Math.h> // min
41 #include <openvdb/Types.h> // Vec3 types and version number
42 #include <openvdb/Grid.h> // grid
43 #include <openvdb/util/NullInterrupter.h>
44 #include "Interpolation.h" // sampling
45 
46 #include <boost/static_assert.hpp>
47 #include <tbb/blocked_range.h> // threading
48 #include <tbb/parallel_for.h> // threading
49 #include <tbb/task.h> // for cancel
50 
51 
52 namespace openvdb {
54 namespace OPENVDB_VERSION_NAME {
55 namespace tools {
56 
60 template<typename CptGridT = Vec3fGrid>
62 {
63 public:
64  typedef CptGridT CptGridType;
65  typedef typename CptGridType::ConstAccessor CptAccessor;
66  typedef typename CptGridType::ValueType CptValueType;
67 
69  mCptIterations(0)
70  {
71  }
72  ClosestPointProjector(const CptGridType& cptGrid, int n):
73  mCptGrid(&cptGrid),
74  mCptAccessor(cptGrid.getAccessor()),
75  mCptIterations(n)
76  {
77  }
79  mCptGrid(other.mCptGrid),
80  mCptAccessor(mCptGrid->getAccessor()),
81  mCptIterations(other.mCptIterations)
82  {
83  }
84  void setConstraintIterations(unsigned int cptIterations) { mCptIterations = cptIterations; }
85  unsigned int numIterations() { return mCptIterations;};
86 
87  // point constraint
88  template <typename LocationType>
89  inline void projectToConstraintSurface(LocationType& W) const
90  {
95  CptValueType result(W[0], W[1],W[2]);
96  for (unsigned int i = 0; i < mCptIterations; ++i) {
97  const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
98  BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
99  }
100  W[0] = result[0];
101  W[1] = result[1];
102  W[2] = result[2];
103  }
104 
105 private:
106  const CptGridType* mCptGrid; // Closest-Point-Transform vector field
107  CptAccessor mCptAccessor;
108  unsigned int mCptIterations;
109 };// end of ClosestPointProjector class
110 
111 
115 template<typename GridT = Vec3fGrid, bool StaggeredVelocity = false>
117 {
118 public:
119  typedef typename GridT::ConstAccessor VelAccessor;
120  typedef typename GridT::ValueType VelValueType;
121 
122  VelocitySampler(const GridT& velGrid):
123  mVelGrid(&velGrid),
124  mVelAccessor(mVelGrid->getAccessor())
125  {
126  }
128  mVelGrid(other.mVelGrid),
129  mVelAccessor(mVelGrid->getAccessor())
130  {
131  }
133  {
134  }
137  template <typename LocationType>
138  inline void sample(const LocationType& W, VelValueType& result) const
139  {
140  const Vec3R location = mVelGrid->worldToIndex(Vec3R(W[0], W[1], W[2]));
142  if (StaggeredVelocity) {
143  // the velocity Grid stores data in MAC-style staggered layout
144  StaggeredBoxSampler::sample<VelAccessor>(mVelAccessor, location, result);
145  } else {
146  // the velocity Grid uses collocated data
147  BoxSampler::sample<VelAccessor>(mVelAccessor, location, result);
148  }
149  }
150 
151 private:
152  // holding the Grids for the transforms
153  const GridT* mVelGrid; // Velocity vector field
154  VelAccessor mVelAccessor;
155 };// end of VelocitySampler class
156 
157 
160 template<typename GridT = Vec3fGrid, bool StaggeredVelocity = false>
162 {
163 public:
164  typedef typename GridT::ValueType VecType;
165  typedef typename VecType::ValueType ElementType;
166 
167  VelocityIntegrator(const GridT& velGrid):
168  mVelField(velGrid)
169  {
170  }
171  // variable order Runge-Kutta time integration for a single time step
172  template<int Order, typename LocationType>
173  void rungeKutta(const float dt, LocationType& loc) {
174  VecType P(loc[0],loc[1],loc[2]), V0, V1, V2, V3;
175 
176  BOOST_STATIC_ASSERT((Order < 5) && (Order > -1));
178  if (Order == 0) {
179  // do nothing
180  return ;
181  } else if (Order == 1) {
182  mVelField.sample(P, V0);
183  P = dt*V0;
184 
185  } else if (Order == 2) {
186  mVelField.sample(P, V0);
187  mVelField.sample(P + ElementType(0.5) * ElementType(dt) * V0, V1);
188  P = dt*V1;
189 
190  } else if (Order == 3) {
191  mVelField.sample(P, V0);
192  mVelField.sample(P+ElementType(0.5)*ElementType(dt)*V0, V1);
193  mVelField.sample(P+dt*(ElementType(2.0)*V1-V0), V2);
194  P = dt*(V0 + ElementType(4.0)*V1 + V2)*ElementType(1.0/6.0);
195 
196  } else if (Order == 4) {
197  mVelField.sample(P, V0);
198  mVelField.sample(P+ElementType(0.5)*ElementType(dt)*V0, V1);
199  mVelField.sample(P+ElementType(0.5)*ElementType(dt)*V1, V2);
200  mVelField.sample(P+ dt*V2, V3);
201  P = dt*(V0 + ElementType(2.0)*(V1 + V2) + V3)*ElementType(1.0/6.0);
202 
203  }
204  loc += LocationType(P[0], P[1], P[2]);
205 
206  }
207 private:
209 };// end of VelocityIntegrator class
210 
211 
213 
214 
236 template<typename GridT = Vec3fGrid,
237  typename PointListT = std::vector<typename GridT::ValueType>,
238  bool StaggeredVelocity = false,
239  typename InterrupterType = util::NullInterrupter>
241 {
242 public:
243  typedef GridT GridType;
244  typedef PointListT PointListType;
245  typedef typename PointListT::value_type LocationType;
247 
248  PointAdvect(const GridT& velGrid, InterrupterType* interrupter=NULL) :
249  mVelGrid(&velGrid),
250  mPoints(NULL),
251  mIntegrationOrder(1),
252  mThreaded(true),
253  mInterrupter(interrupter)
254  {
255  }
256  PointAdvect(const PointAdvect &other) :
257  mVelGrid(other.mVelGrid),
258  mPoints(other.mPoints),
259  mDt(other.mDt),
260  mAdvIterations(other.mAdvIterations),
261  mIntegrationOrder(other.mIntegrationOrder),
262  mThreaded(other.mThreaded),
263  mInterrupter(other.mInterrupter)
264  {
265  }
266  virtual ~PointAdvect()
267  {
268  }
270  bool earlyOut() const { return (mIntegrationOrder==0);}
272  void setThreaded(bool threaded) { mThreaded = threaded; }
273  bool getThreaded() { return mThreaded; }
274  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
275 
277  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
278  {
279  if (this->earlyOut()) return; // nothing to do!
280  mPoints = &points;
281  mDt = dt;
282  mAdvIterations = advIterations;
283 
284  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
285  if (mThreaded) {
286  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *this);
287  } else {
288  (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
289  }
290  if (mInterrupter) mInterrupter->end();
291  }
292 
294  void operator() (const tbb::blocked_range<size_t> &range) const
295  {
296  if (mInterrupter && mInterrupter->wasInterrupted()) {
297  tbb::task::self().cancel_group_execution();
298  }
299 
300  VelocityFieldIntegrator velField(*mVelGrid);
301  switch (mIntegrationOrder) {
302  case 1:
303  {
304  for (size_t n = range.begin(); n != range.end(); ++n) {
305  LocationType& X0 = (*mPoints)[n];
306  // loop over number of time steps
307  for (unsigned int i = 0; i < mAdvIterations; ++i) {
308  velField.template rungeKutta<1>(mDt, X0);
309  }
310  }
311  }
312  break;
313  case 2:
314  {
315  for (size_t n = range.begin(); n != range.end(); ++n) {
316  LocationType& X0 = (*mPoints)[n];
317  // loop over number of time steps
318  for (unsigned int i = 0; i < mAdvIterations; ++i) {
319  velField.template rungeKutta<2>(mDt, X0);
320  }
321  }
322  }
323  break;
324  case 3:
325  {
326  for (size_t n = range.begin(); n != range.end(); ++n) {
327  LocationType& X0 = (*mPoints)[n];
328  // loop over number of time steps
329  for (unsigned int i = 0; i < mAdvIterations; ++i) {
330  velField.template rungeKutta<3>(mDt, X0);
331  }
332  }
333  }
334  break;
335  case 4:
336  {
337  for (size_t n = range.begin(); n != range.end(); ++n) {
338  LocationType& X0 = (*mPoints)[n];
339  // loop over number of time steps
340  for (unsigned int i = 0; i < mAdvIterations; ++i) {
341  velField.template rungeKutta<4>(mDt, X0);
342  }
343  }
344  }
345  break;
346  }
347  }
348 
349 private:
350  // the velocity field
351  const GridType* mVelGrid;
352 
353  // vertex list of all the points
354  PointListT* mPoints;
355 
356  // time integration parameters
357  float mDt; // time step
358  unsigned int mAdvIterations; // number of time steps
359  unsigned int mIntegrationOrder;
360 
361  // operational parameters
362  bool mThreaded;
363  InterrupterType* mInterrupter;
364 
365 };//end of PointAdvect class
366 
367 
368 template<typename GridT = Vec3fGrid,
369  typename PointListT = std::vector<typename GridT::ValueType>,
370  bool StaggeredVelocity = false,
371  typename CptGridType = GridT,
372  typename InterrupterType = util::NullInterrupter>
374 {
375 public:
376  typedef GridT GridType;
377  typedef typename PointListT::value_type LocationType;
380  typedef PointListT PointListType;
381 
383  const GridType& cptGrid, int cptn, InterrupterType* interrupter = NULL):
384  mVelGrid(&velGrid),
385  mCptGrid(&cptGrid),
386  mCptIter(cptn),
387  mInterrupter(interrupter)
388  {
389  }
391  mVelGrid(other.mVelGrid),
392  mCptGrid(other.mCptGrid),
393  mCptIter(other.mCptIter),
394  mPoints(other.mPoints),
395  mDt(other.mDt),
396  mAdvIterations(other.mAdvIterations),
397  mIntegrationOrder(other.mIntegrationOrder),
398  mThreaded(other.mThreaded),
399  mInterrupter(other.mInterrupter)
400  {
401  }
403 
404  void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
405  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
406 
407  void setThreaded(bool threaded) { mThreaded = threaded; }
408  bool getThreaded() { return mThreaded; }
409 
411  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
412  {
413  mPoints = &points;
414  mDt = dt;
415 
416  if (mIntegrationOrder==0 && mCptIter == 0) {
417  return; // nothing to do!
418  }
419  (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
420 
421  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
422  const int N = mPoints->size();
423 
424  if (mThreaded) {
425  tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
426  } else {
427  (*this)(tbb::blocked_range<size_t>(0, N));
428  }
429  if (mInterrupter) mInterrupter->end();
430  }
431 
432 
434  void operator() (const tbb::blocked_range<size_t> &range) const
435  {
436  if (mInterrupter && mInterrupter->wasInterrupted()) {
437  tbb::task::self().cancel_group_execution();
438  }
439 
440  VelocityIntegratorType velField(*mVelGrid);
441  ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
442  switch (mIntegrationOrder) {
443  case 0://pure CPT projection
444  {
445  for (size_t n = range.begin(); n != range.end(); ++n) {
446  LocationType& X0 = (*mPoints)[n];
447  for (unsigned int i = 0; i < mAdvIterations; ++i) {
448  cptField.projectToConstraintSurface(X0);
449  }
450  }
451  }
452  break;
453  case 1://1'th order advection and CPT projection
454  {
455  for (size_t n = range.begin(); n != range.end(); ++n) {
456  LocationType& X0 = (*mPoints)[n];
457  for (unsigned int i = 0; i < mAdvIterations; ++i) {
458  velField.template rungeKutta<1>(mDt, X0);
459  cptField.projectToConstraintSurface(X0);
460  }
461  }
462  }
463  break;
464  case 2://2'nd order advection and CPT projection
465  {
466  for (size_t n = range.begin(); n != range.end(); ++n) {
467  LocationType& X0 = (*mPoints)[n];
468  for (unsigned int i = 0; i < mAdvIterations; ++i) {
469  velField.template rungeKutta<2>(mDt, X0);
470  cptField.projectToConstraintSurface(X0);
471  }
472  }
473  }
474  break;
475 
476  case 3://3'rd order advection and CPT projection
477  {
478  for (size_t n = range.begin(); n != range.end(); ++n) {
479  LocationType& X0 = (*mPoints)[n];
480  for (unsigned int i = 0; i < mAdvIterations; ++i) {
481  velField.template rungeKutta<3>(mDt, X0);
482  cptField.projectToConstraintSurface(X0);
483  }
484  }
485  }
486  break;
487  case 4://4'th order advection and CPT projection
488  {
489  for (size_t n = range.begin(); n != range.end(); ++n) {
490  LocationType& X0 = (*mPoints)[n];
491  for (unsigned int i = 0; i < mAdvIterations; ++i) {
492  velField.template rungeKutta<4>(mDt, X0);
493  cptField.projectToConstraintSurface(X0);
494  }
495  }
496  }
497  break;
498  }
499  }
500 
501 private:
502  const GridType* mVelGrid; // the velocity field
503  const GridType* mCptGrid;
504  int mCptIter;
505  PointListT* mPoints; // vertex list of all the points
506 
507  // time integration parameters
508  float mDt; // time step
509  unsigned int mAdvIterations; // number of time steps
510  unsigned int mIntegrationOrder; // order of Runge-Kutta integration
511  // operational parameters
512  bool mThreaded;
513  InterrupterType* mInterrupter;
514 };// end of ConstrainedPointAdvect class
515 
516 } // namespace tools
517 } // namespace OPENVDB_VERSION_NAME
518 } // namespace openvdb
519 
520 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
521 
522 // Copyright (c) 2012-2013 DreamWorks Animation LLC
523 // All rights reserved. This software is distributed under the
524 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
Performs runge-kutta time integration of variable order in a static velocity field.
Definition: PointAdvect.h:161
GridT::ValueType VelValueType
Definition: PointAdvect.h:120
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained advection of a list of points over a time = dt * advIterations.
Definition: PointAdvect.h:277
bool getThreaded()
Definition: PointAdvect.h:408
CptGridType::ConstAccessor CptAccessor
Definition: PointAdvect.h:65
GridT GridType
Definition: PointAdvect.h:376
unsigned int numIterations()
Definition: PointAdvect.h:85
virtual ~PointAdvect()
Definition: PointAdvect.h:266
void sample(const LocationType &W, VelValueType &result) const
Definition: PointAdvect.h:138
void projectToConstraintSurface(LocationType &W) const
Definition: PointAdvect.h:89
void setConstraintIterations(unsigned int cptIterations)
Definition: PointAdvect.h:84
VecType::ValueType ElementType
Definition: PointAdvect.h:165
void rungeKutta(const float dt, LocationType &loc)
Definition: PointAdvect.h:173
Definition: PointAdvect.h:116
void setIntegrationOrder(unsigned int order)
Definition: PointAdvect.h:274
PointAdvect(const PointAdvect &other)
Definition: PointAdvect.h:256
GridT::ValueType VecType
Definition: PointAdvect.h:164
VelocitySampler(const GridT &velGrid)
Definition: PointAdvect.h:122
void setConstraintIterations(unsigned int cptIter)
Definition: PointAdvect.h:404
~VelocitySampler()
Definition: PointAdvect.h:132
virtual ~ConstrainedPointAdvect()
Definition: PointAdvect.h:402
VelocitySampler(const VelocitySampler &other)
Definition: PointAdvect.h:127
GridT GridType
Definition: PointAdvect.h:243
ClosestPointProjector< CptGridType > ClosestPointProjectorType
Definition: PointAdvect.h:379
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: PointAdvect.h:390
#define OPENVDB_VERSION_NAME
Definition: version.h:45
math::Vec3< Real > Vec3R
Definition: Types.h:74
CptGridT CptGridType
Definition: PointAdvect.h:64
void setThreaded(bool threaded)
get &amp; set
Definition: PointAdvect.h:272
ClosestPointProjector(const ClosestPointProjector &other)
Definition: PointAdvect.h:78
VelocityIntegrator(const GridT &velGrid)
Definition: PointAdvect.h:167
void setIntegrationOrder(unsigned int order)
Definition: PointAdvect.h:405
Definition: PointAdvect.h:240
PointListT PointListType
Definition: PointAdvect.h:244
bool getThreaded()
Definition: PointAdvect.h:273
CptGridType::ValueType CptValueType
Definition: PointAdvect.h:66
PointListT::value_type LocationType
Definition: PointAdvect.h:245
VelocityIntegrator< GridT, StaggeredVelocity > VelocityIntegratorType
Definition: PointAdvect.h:378
PointListT::value_type LocationType
Definition: PointAdvect.h:377
bool earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: PointAdvect.h:270
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=NULL)
Definition: PointAdvect.h:382
Vec3SGrid Vec3fGrid
Definition: openvdb.h:79
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained Advection a list of points over a time = dt * advIterations.
Definition: PointAdvect.h:411
PointListT PointListType
Definition: PointAdvect.h:380
ClosestPointProjector()
Definition: PointAdvect.h:68
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=NULL)
Definition: PointAdvect.h:248
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
GridT::ConstAccessor VelAccessor
Definition: PointAdvect.h:119
void setThreaded(bool threaded)
Definition: PointAdvect.h:407
VelocityIntegrator< GridT, StaggeredVelocity > VelocityFieldIntegrator
Definition: PointAdvect.h:246
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: PointAdvect.h:72
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52