OpenVDB  6.2.0
tools/PointAdvect.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 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
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) 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 earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: tools/PointAdvect.h:170
Definition: tools/PointAdvect.h:62
ClosestPointProjector(const ClosestPointProjector &other)
Definition: tools/PointAdvect.h:79
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Defines two simple wrapper classes for advection velocity fields as well as VelocitySampler and Veloc...
PointListT PointListType
Definition: tools/PointAdvect.h:144
GridT GridType
Definition: tools/PointAdvect.h:276
bool getThreaded()
Definition: tools/PointAdvect.h:308
void setConstraintIterations(unsigned int cptIterations)
Definition: tools/PointAdvect.h:85
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:305
math::Vec3< Real > Vec3R
Definition: Types.h:79
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
Vec3SGrid Vec3fGrid
Definition: openvdb.h:83
Definition: tools/PointAdvect.h:273
CptGridT CptGridType
Definition: tools/PointAdvect.h:65
PointAdvect(const PointAdvect &other)
Definition: tools/PointAdvect.h:156
void projectToConstraintSurface(LocationType &W) const
Definition: tools/PointAdvect.h:90
PointListT PointListType
Definition: tools/PointAdvect.h:280
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
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:128
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:145
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: tools/PointAdvect.h:290
bool getThreaded()
Definition: tools/PointAdvect.h:173
Definition: Exceptions.h:40
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:277
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:52
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:174
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:282
typename CptGridType::ValueType CptValueType
Definition: tools/PointAdvect.h:67
ClosestPointProjector()
Definition: tools/PointAdvect.h:69
GridT GridType
Definition: tools/PointAdvect.h:143
Definition: tools/PointAdvect.h:140
void setThreaded(bool threaded)
Definition: tools/PointAdvect.h:307
virtual ~PointAdvect()
Definition: tools/PointAdvect.h:166
unsigned int numIterations()
Definition: tools/PointAdvect.h:86
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:180
void setConstraintIterations(unsigned int cptIter)
Definition: tools/PointAdvect.h:304
void setThreaded(bool threaded)
get & set
Definition: tools/PointAdvect.h:172
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: tools/PointAdvect.h:73
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:241
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:148
typename CptGridType::ConstAccessor CptAccessor
Definition: tools/PointAdvect.h:66
virtual ~ConstrainedPointAdvect()
Definition: tools/PointAdvect.h:302