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