OpenVDB  5.2.0
tools/PointAdvect.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2018 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 //
36 
37 #ifndef OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
38 #define OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
39 
40 #include <openvdb/openvdb.h>
41 #include <openvdb/math/Math.h> // min
42 #include <openvdb/Types.h> // Vec3 types and version number
43 #include <openvdb/Grid.h> // grid
44 #include <openvdb/util/NullInterrupter.h>
45 #include "Interpolation.h" // sampling
46 #include "VelocityFields.h" // VelocityIntegrator
47 #include <tbb/blocked_range.h> // threading
48 #include <tbb/parallel_for.h> // threading
49 #include <tbb/task.h> // for cancel
50 #include <vector>
51 
52 
53 namespace openvdb {
55 namespace OPENVDB_VERSION_NAME {
56 namespace tools {
57 
61 template<typename CptGridT = Vec3fGrid>
63 {
64 public:
65  using CptGridType = CptGridT;
66  using CptAccessor = typename CptGridType::ConstAccessor;
67  using CptValueType = typename CptGridType::ValueType;
68 
70  mCptIterations(0)
71  {
72  }
73  ClosestPointProjector(const CptGridType& cptGrid, int n):
74  mCptGrid(&cptGrid),
75  mCptAccessor(cptGrid.getAccessor()),
76  mCptIterations(n)
77  {
78  }
80  mCptGrid(other.mCptGrid),
81  mCptAccessor(mCptGrid->getAccessor()),
82  mCptIterations(other.mCptIterations)
83  {
84  }
85  void setConstraintIterations(unsigned int cptIterations) { mCptIterations = cptIterations; }
86  unsigned int numIterations() { return mCptIterations; }
87 
88  // point constraint
89  template <typename LocationType>
90  inline void projectToConstraintSurface(LocationType& W) const
91  {
96  CptValueType result(W[0], W[1],W[2]);
97  for (unsigned int i = 0; i < mCptIterations; ++i) {
98  const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
99  BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
100  }
101  W[0] = result[0];
102  W[1] = result[1];
103  W[2] = result[2];
104  }
105 
106 private:
107  const CptGridType* mCptGrid; // Closest-Point-Transform vector field
108  CptAccessor mCptAccessor;
109  unsigned int mCptIterations;
110 };// end of ClosestPointProjector class
111 
113 
114 
136 template<typename GridT = Vec3fGrid,
137  typename PointListT = std::vector<typename GridT::ValueType>,
138  bool StaggeredVelocity = false,
139  typename InterrupterType = util::NullInterrupter>
141 {
142 public:
143  using GridType = GridT;
144  using PointListType = PointListT;
145  using LocationType = typename PointListT::value_type;
147 
148  PointAdvect(const GridT& velGrid, InterrupterType* interrupter = nullptr):
149  mVelGrid(&velGrid),
150  mPoints(nullptr),
151  mIntegrationOrder(1),
152  mThreaded(true),
153  mInterrupter(interrupter)
154  {
155  }
156  PointAdvect(const PointAdvect &other) :
157  mVelGrid(other.mVelGrid),
158  mPoints(other.mPoints),
159  mDt(other.mDt),
160  mAdvIterations(other.mAdvIterations),
161  mIntegrationOrder(other.mIntegrationOrder),
162  mThreaded(other.mThreaded),
163  mInterrupter(other.mInterrupter)
164  {
165  }
166  virtual ~PointAdvect()
167  {
168  }
170  bool earlyOut() const { return (mIntegrationOrder==0);}
172  void setThreaded(bool threaded) { mThreaded = threaded; }
173  bool getThreaded() { return mThreaded; }
174  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
175 
177  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
178  {
179  if (this->earlyOut()) return; // nothing to do!
180  mPoints = &points;
181  mDt = dt;
182  mAdvIterations = advIterations;
183 
184  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
185  if (mThreaded) {
186  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *this);
187  } else {
188  (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
189  }
190  if (mInterrupter) mInterrupter->end();
191  }
192 
194  void operator() (const tbb::blocked_range<size_t> &range) const
195  {
196  if (mInterrupter && mInterrupter->wasInterrupted()) {
197  tbb::task::self().cancel_group_execution();
198  }
199 
200  VelocityFieldIntegrator velField(*mVelGrid);
201  switch (mIntegrationOrder) {
202  case 1:
203  {
204  for (size_t n = range.begin(); n != range.end(); ++n) {
205  LocationType& X0 = (*mPoints)[n];
206  // loop over number of time steps
207  for (unsigned int i = 0; i < mAdvIterations; ++i) {
208  velField.template rungeKutta<1>(mDt, X0);
209  }
210  }
211  }
212  break;
213  case 2:
214  {
215  for (size_t n = range.begin(); n != range.end(); ++n) {
216  LocationType& X0 = (*mPoints)[n];
217  // loop over number of time steps
218  for (unsigned int i = 0; i < mAdvIterations; ++i) {
219  velField.template rungeKutta<2>(mDt, X0);
220  }
221  }
222  }
223  break;
224  case 3:
225  {
226  for (size_t n = range.begin(); n != range.end(); ++n) {
227  LocationType& X0 = (*mPoints)[n];
228  // loop over number of time steps
229  for (unsigned int i = 0; i < mAdvIterations; ++i) {
230  velField.template rungeKutta<3>(mDt, X0);
231  }
232  }
233  }
234  break;
235  case 4:
236  {
237  for (size_t n = range.begin(); n != range.end(); ++n) {
238  LocationType& X0 = (*mPoints)[n];
239  // loop over number of time steps
240  for (unsigned int i = 0; i < mAdvIterations; ++i) {
241  velField.template rungeKutta<4>(mDt, X0);
242  }
243  }
244  }
245  break;
246  }
247  }
248 
249 private:
250  // the velocity field
251  const GridType* mVelGrid;
252 
253  // vertex list of all the points
254  PointListT* mPoints;
255 
256  // time integration parameters
257  float mDt; // time step
258  unsigned int mAdvIterations; // number of time steps
259  unsigned int mIntegrationOrder;
260 
261  // operational parameters
262  bool mThreaded;
263  InterrupterType* mInterrupter;
264 
265 };//end of PointAdvect class
266 
267 
268 template<typename GridT = Vec3fGrid,
269  typename PointListT = std::vector<typename GridT::ValueType>,
270  bool StaggeredVelocity = false,
271  typename CptGridType = GridT,
272  typename InterrupterType = util::NullInterrupter>
274 {
275 public:
276  using GridType = GridT;
277  using LocationType = typename PointListT::value_type;
280  using PointListType = PointListT;
281 
283  const GridType& cptGrid, int cptn, InterrupterType* interrupter = nullptr):
284  mVelGrid(&velGrid),
285  mCptGrid(&cptGrid),
286  mCptIter(cptn),
287  mInterrupter(interrupter)
288  {
289  }
291  mVelGrid(other.mVelGrid),
292  mCptGrid(other.mCptGrid),
293  mCptIter(other.mCptIter),
294  mPoints(other.mPoints),
295  mDt(other.mDt),
296  mAdvIterations(other.mAdvIterations),
297  mIntegrationOrder(other.mIntegrationOrder),
298  mThreaded(other.mThreaded),
299  mInterrupter(other.mInterrupter)
300  {
301  }
303 
304  void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
305  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
306 
307  void setThreaded(bool threaded) { mThreaded = threaded; }
308  bool getThreaded() { return mThreaded; }
309 
311  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
312  {
313  mPoints = &points;
314  mDt = dt;
315 
316  if (mIntegrationOrder==0 && mCptIter == 0) {
317  return; // nothing to do!
318  }
319  (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
320 
321  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
322  const size_t N = mPoints->size();
323 
324  if (mThreaded) {
325  tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
326  } else {
327  (*this)(tbb::blocked_range<size_t>(0, N));
328  }
329  if (mInterrupter) mInterrupter->end();
330  }
331 
332 
334  void operator() (const tbb::blocked_range<size_t> &range) const
335  {
336  if (mInterrupter && mInterrupter->wasInterrupted()) {
337  tbb::task::self().cancel_group_execution();
338  }
339 
340  VelocityIntegratorType velField(*mVelGrid);
341  ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
342  switch (mIntegrationOrder) {
343  case 0://pure CPT projection
344  {
345  for (size_t n = range.begin(); n != range.end(); ++n) {
346  LocationType& X0 = (*mPoints)[n];
347  for (unsigned int i = 0; i < mAdvIterations; ++i) {
348  cptField.projectToConstraintSurface(X0);
349  }
350  }
351  }
352  break;
353  case 1://1'th order advection and CPT projection
354  {
355  for (size_t n = range.begin(); n != range.end(); ++n) {
356  LocationType& X0 = (*mPoints)[n];
357  for (unsigned int i = 0; i < mAdvIterations; ++i) {
358  velField.template rungeKutta<1>(mDt, X0);
359  cptField.projectToConstraintSurface(X0);
360  }
361  }
362  }
363  break;
364  case 2://2'nd order advection and CPT projection
365  {
366  for (size_t n = range.begin(); n != range.end(); ++n) {
367  LocationType& X0 = (*mPoints)[n];
368  for (unsigned int i = 0; i < mAdvIterations; ++i) {
369  velField.template rungeKutta<2>(mDt, X0);
370  cptField.projectToConstraintSurface(X0);
371  }
372  }
373  }
374  break;
375 
376  case 3://3'rd order advection and CPT projection
377  {
378  for (size_t n = range.begin(); n != range.end(); ++n) {
379  LocationType& X0 = (*mPoints)[n];
380  for (unsigned int i = 0; i < mAdvIterations; ++i) {
381  velField.template rungeKutta<3>(mDt, X0);
382  cptField.projectToConstraintSurface(X0);
383  }
384  }
385  }
386  break;
387  case 4://4'th order advection and CPT projection
388  {
389  for (size_t n = range.begin(); n != range.end(); ++n) {
390  LocationType& X0 = (*mPoints)[n];
391  for (unsigned int i = 0; i < mAdvIterations; ++i) {
392  velField.template rungeKutta<4>(mDt, X0);
393  cptField.projectToConstraintSurface(X0);
394  }
395  }
396  }
397  break;
398  }
399  }
400 
401 private:
402  const GridType* mVelGrid; // the velocity field
403  const GridType* mCptGrid;
404  int mCptIter;
405  PointListT* mPoints; // vertex list of all the points
406 
407  // time integration parameters
408  float mDt; // time step
409  unsigned int mAdvIterations; // number of time steps
410  unsigned int mIntegrationOrder; // order of Runge-Kutta integration
411  // operational parameters
412  bool mThreaded;
413  InterrupterType* mInterrupter;
414 };// end of ConstrainedPointAdvect class
415 
416 } // namespace tools
417 } // namespace OPENVDB_VERSION_NAME
418 } // namespace openvdb
419 
420 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
421 
422 // Copyright (c) 2012-2018 DreamWorks Animation LLC
423 // All rights reserved. This software is distributed under the
424 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
bool getThreaded()
Definition: tools/PointAdvect.h:308
virtual ~PointAdvect()
Definition: tools/PointAdvect.h:166
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:241
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained advection of a list of points over a time = dt * advIterations.
Definition: tools/PointAdvect.h:177
unsigned int numIterations()
Definition: tools/PointAdvect.h:86
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
void setConstraintIterations(unsigned int cptIter)
Definition: tools/PointAdvect.h:304
ClosestPointProjector()
Definition: tools/PointAdvect.h:69
void setThreaded(bool threaded)
Definition: tools/PointAdvect.h:307
CptGridT CptGridType
Definition: tools/PointAdvect.h:65
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:305
Definition: tools/PointAdvect.h:273
bool getThreaded()
Definition: tools/PointAdvect.h:173
Definition: tools/PointAdvect.h:62
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: tools/PointAdvect.h:73
typename CptGridType::ValueType CptValueType
Definition: tools/PointAdvect.h:67
Vec3SGrid Vec3fGrid
Definition: openvdb.h:83
void setThreaded(bool threaded)
get & set
Definition: tools/PointAdvect.h:172
virtual ~ConstrainedPointAdvect()
Definition: tools/PointAdvect.h:302
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:136
math::Vec3< Real > Vec3R
Definition: Types.h:79
bool earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: tools/PointAdvect.h:170
void projectToConstraintSurface(LocationType &W) const
Definition: tools/PointAdvect.h:90
Definition: Exceptions.h:40
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: tools/PointAdvect.h:290
GridT GridType
Definition: tools/PointAdvect.h:143
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:148
typename CptGridType::ConstAccessor CptAccessor
Definition: tools/PointAdvect.h:66
Definition: tools/PointAdvect.h:140
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained Advection a list of points over a time = dt * advIterations.
Definition: tools/PointAdvect.h:311
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:174
ClosestPointProjector(const ClosestPointProjector &other)
Definition: tools/PointAdvect.h:79
void setConstraintIterations(unsigned int cptIterations)
Definition: tools/PointAdvect.h:85
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:188
PointAdvect(const PointAdvect &other)
Definition: tools/PointAdvect.h:156
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:282
PointListT PointListType
Definition: tools/PointAdvect.h:144
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:145
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
GridT GridType
Definition: tools/PointAdvect.h:276
PointListT PointListType
Definition: tools/PointAdvect.h:280
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:277