OpenVDB  8.0.1
tools/PointAdvect.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 //
9 
10 #ifndef OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
11 #define OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
12 
13 #include <openvdb/openvdb.h>
14 #include <openvdb/math/Math.h> // min
15 #include <openvdb/Types.h> // Vec3 types and version number
16 #include <openvdb/Grid.h> // grid
18 #include "Interpolation.h" // sampling
19 #include "VelocityFields.h" // VelocityIntegrator
20 #include <tbb/blocked_range.h> // threading
21 #include <tbb/parallel_for.h> // threading
22 #include <tbb/task.h> // for cancel
23 #include <vector>
24 
25 
26 namespace openvdb {
28 namespace OPENVDB_VERSION_NAME {
29 namespace tools {
30 
34 template<typename CptGridT = Vec3fGrid>
36 {
37 public:
38  using CptGridType = CptGridT;
39  using CptAccessor = typename CptGridType::ConstAccessor;
40  using CptValueType = typename CptGridType::ValueType;
41 
43  mCptIterations(0)
44  {
45  }
46  ClosestPointProjector(const CptGridType& cptGrid, int n):
47  mCptGrid(&cptGrid),
48  mCptAccessor(cptGrid.getAccessor()),
49  mCptIterations(n)
50  {
51  }
53  mCptGrid(other.mCptGrid),
54  mCptAccessor(mCptGrid->getAccessor()),
55  mCptIterations(other.mCptIterations)
56  {
57  }
58  void setConstraintIterations(unsigned int cptIterations) { mCptIterations = cptIterations; }
59  unsigned int numIterations() { return mCptIterations; }
60 
61  // point constraint
62  template <typename LocationType>
63  inline void projectToConstraintSurface(LocationType& W) const
64  {
69  CptValueType result(W[0], W[1],W[2]);
70  for (unsigned int i = 0; i < mCptIterations; ++i) {
71  const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
72  BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
73  }
74  W[0] = result[0];
75  W[1] = result[1];
76  W[2] = result[2];
77  }
78 
79 private:
80  const CptGridType* mCptGrid; // Closest-Point-Transform vector field
81  CptAccessor mCptAccessor;
82  unsigned int mCptIterations;
83 };// end of ClosestPointProjector class
84 
86 
87 
109 template<typename GridT = Vec3fGrid,
110  typename PointListT = std::vector<typename GridT::ValueType>,
111  bool StaggeredVelocity = false,
112  typename InterrupterType = util::NullInterrupter>
114 {
115 public:
116  using GridType = GridT;
117  using PointListType = PointListT;
118  using LocationType = typename PointListT::value_type;
120 
121  PointAdvect(const GridT& velGrid, InterrupterType* interrupter = nullptr):
122  mVelGrid(&velGrid),
123  mPoints(nullptr),
124  mIntegrationOrder(1),
125  mThreaded(true),
126  mInterrupter(interrupter)
127  {
128  }
129  PointAdvect(const PointAdvect &other) :
130  mVelGrid(other.mVelGrid),
131  mPoints(other.mPoints),
132  mDt(other.mDt),
133  mAdvIterations(other.mAdvIterations),
134  mIntegrationOrder(other.mIntegrationOrder),
135  mThreaded(other.mThreaded),
136  mInterrupter(other.mInterrupter)
137  {
138  }
139  virtual ~PointAdvect()
140  {
141  }
143  bool earlyOut() const { return (mIntegrationOrder==0);}
145  void setThreaded(bool threaded) { mThreaded = threaded; }
146  bool getThreaded() { return mThreaded; }
147  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
148 
150  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
151  {
152  if (this->earlyOut()) return; // nothing to do!
153  mPoints = &points;
154  mDt = dt;
155  mAdvIterations = advIterations;
156 
157  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
158  if (mThreaded) {
159  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *this);
160  } else {
161  (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
162  }
163  if (mInterrupter) mInterrupter->end();
164  }
165 
167  void operator() (const tbb::blocked_range<size_t> &range) const
168  {
169  if (mInterrupter && mInterrupter->wasInterrupted()) {
170  tbb::task::self().cancel_group_execution();
171  }
172 
173  VelocityFieldIntegrator velField(*mVelGrid);
174  switch (mIntegrationOrder) {
175  case 1:
176  {
177  for (size_t n = range.begin(); n != range.end(); ++n) {
178  LocationType& X0 = (*mPoints)[n];
179  // loop over number of time steps
180  for (unsigned int i = 0; i < mAdvIterations; ++i) {
181  velField.template rungeKutta<1>(mDt, X0);
182  }
183  }
184  }
185  break;
186  case 2:
187  {
188  for (size_t n = range.begin(); n != range.end(); ++n) {
189  LocationType& X0 = (*mPoints)[n];
190  // loop over number of time steps
191  for (unsigned int i = 0; i < mAdvIterations; ++i) {
192  velField.template rungeKutta<2>(mDt, X0);
193  }
194  }
195  }
196  break;
197  case 3:
198  {
199  for (size_t n = range.begin(); n != range.end(); ++n) {
200  LocationType& X0 = (*mPoints)[n];
201  // loop over number of time steps
202  for (unsigned int i = 0; i < mAdvIterations; ++i) {
203  velField.template rungeKutta<3>(mDt, X0);
204  }
205  }
206  }
207  break;
208  case 4:
209  {
210  for (size_t n = range.begin(); n != range.end(); ++n) {
211  LocationType& X0 = (*mPoints)[n];
212  // loop over number of time steps
213  for (unsigned int i = 0; i < mAdvIterations; ++i) {
214  velField.template rungeKutta<4>(mDt, X0);
215  }
216  }
217  }
218  break;
219  }
220  }
221 
222 private:
223  // the velocity field
224  const GridType* mVelGrid;
225 
226  // vertex list of all the points
227  PointListT* mPoints;
228 
229  // time integration parameters
230  float mDt; // time step
231  unsigned int mAdvIterations; // number of time steps
232  unsigned int mIntegrationOrder;
233 
234  // operational parameters
235  bool mThreaded;
236  InterrupterType* mInterrupter;
237 
238 };//end of PointAdvect class
239 
240 
241 template<typename GridT = Vec3fGrid,
242  typename PointListT = std::vector<typename GridT::ValueType>,
243  bool StaggeredVelocity = false,
244  typename CptGridType = GridT,
245  typename InterrupterType = util::NullInterrupter>
247 {
248 public:
249  using GridType = GridT;
250  using LocationType = typename PointListT::value_type;
253  using PointListType = PointListT;
254 
256  const GridType& cptGrid, int cptn, InterrupterType* interrupter = nullptr):
257  mVelGrid(&velGrid),
258  mCptGrid(&cptGrid),
259  mCptIter(cptn),
260  mInterrupter(interrupter)
261  {
262  }
264  mVelGrid(other.mVelGrid),
265  mCptGrid(other.mCptGrid),
266  mCptIter(other.mCptIter),
267  mPoints(other.mPoints),
268  mDt(other.mDt),
269  mAdvIterations(other.mAdvIterations),
270  mIntegrationOrder(other.mIntegrationOrder),
271  mThreaded(other.mThreaded),
272  mInterrupter(other.mInterrupter)
273  {
274  }
276 
277  void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
278  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
279 
280  void setThreaded(bool threaded) { mThreaded = threaded; }
281  bool getThreaded() { return mThreaded; }
282 
284  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
285  {
286  mPoints = &points;
287  mDt = dt;
288 
289  if (mIntegrationOrder==0 && mCptIter == 0) {
290  return; // nothing to do!
291  }
292  (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
293 
294  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
295  const size_t N = mPoints->size();
296 
297  if (mThreaded) {
298  tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
299  } else {
300  (*this)(tbb::blocked_range<size_t>(0, N));
301  }
302  if (mInterrupter) mInterrupter->end();
303  }
304 
305 
307  void operator() (const tbb::blocked_range<size_t> &range) const
308  {
309  if (mInterrupter && mInterrupter->wasInterrupted()) {
310  tbb::task::self().cancel_group_execution();
311  }
312 
313  VelocityIntegratorType velField(*mVelGrid);
314  ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
315  switch (mIntegrationOrder) {
316  case 0://pure CPT projection
317  {
318  for (size_t n = range.begin(); n != range.end(); ++n) {
319  LocationType& X0 = (*mPoints)[n];
320  for (unsigned int i = 0; i < mAdvIterations; ++i) {
321  cptField.projectToConstraintSurface(X0);
322  }
323  }
324  }
325  break;
326  case 1://1'th order advection and CPT projection
327  {
328  for (size_t n = range.begin(); n != range.end(); ++n) {
329  LocationType& X0 = (*mPoints)[n];
330  for (unsigned int i = 0; i < mAdvIterations; ++i) {
331  velField.template rungeKutta<1>(mDt, X0);
332  cptField.projectToConstraintSurface(X0);
333  }
334  }
335  }
336  break;
337  case 2://2'nd order advection and CPT projection
338  {
339  for (size_t n = range.begin(); n != range.end(); ++n) {
340  LocationType& X0 = (*mPoints)[n];
341  for (unsigned int i = 0; i < mAdvIterations; ++i) {
342  velField.template rungeKutta<2>(mDt, X0);
343  cptField.projectToConstraintSurface(X0);
344  }
345  }
346  }
347  break;
348 
349  case 3://3'rd order advection and CPT projection
350  {
351  for (size_t n = range.begin(); n != range.end(); ++n) {
352  LocationType& X0 = (*mPoints)[n];
353  for (unsigned int i = 0; i < mAdvIterations; ++i) {
354  velField.template rungeKutta<3>(mDt, X0);
355  cptField.projectToConstraintSurface(X0);
356  }
357  }
358  }
359  break;
360  case 4://4'th order advection and CPT projection
361  {
362  for (size_t n = range.begin(); n != range.end(); ++n) {
363  LocationType& X0 = (*mPoints)[n];
364  for (unsigned int i = 0; i < mAdvIterations; ++i) {
365  velField.template rungeKutta<4>(mDt, X0);
366  cptField.projectToConstraintSurface(X0);
367  }
368  }
369  }
370  break;
371  }
372  }
373 
374 private:
375  const GridType* mVelGrid; // the velocity field
376  const GridType* mCptGrid;
377  int mCptIter;
378  PointListT* mPoints; // vertex list of all the points
379 
380  // time integration parameters
381  float mDt; // time step
382  unsigned int mAdvIterations; // number of time steps
383  unsigned int mIntegrationOrder; // order of Runge-Kutta integration
384  // operational parameters
385  bool mThreaded;
386  InterrupterType* mInterrupter;
387 };// end of ConstrainedPointAdvect class
388 
389 } // namespace tools
390 } // namespace OPENVDB_VERSION_NAME
391 } // namespace openvdb
392 
393 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
GridT GridType
Definition: tools/PointAdvect.h:116
PointListT PointListType
Definition: tools/PointAdvect.h:117
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:215
void setConstraintIterations(unsigned int cptIterations)
Definition: tools/PointAdvect.h:58
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:121
PointAdvect(const PointAdvect &other)
Definition: tools/PointAdvect.h:129
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...
Definition: tools/PointAdvect.h:246
void setThreaded(bool threaded)
get & set
Definition: tools/PointAdvect.h:145
Definition: tools/PointAdvect.h:113
void setThreaded(bool threaded)
Definition: tools/PointAdvect.h:280
unsigned int numIterations()
Definition: tools/PointAdvect.h:59
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:255
typename CptGridType::ConstAccessor CptAccessor
Definition: tools/PointAdvect.h:39
virtual ~ConstrainedPointAdvect()
Definition: tools/PointAdvect.h:275
Definition: tools/PointAdvect.h:35
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:118
ClosestPointProjector(const ClosestPointProjector &other)
Definition: tools/PointAdvect.h:52
bool getThreaded()
Definition: tools/PointAdvect.h:146
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:101
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
virtual ~PointAdvect()
Definition: tools/PointAdvect.h:139
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:150
Definition: openvdb/Exceptions.h:13
void setConstraintIterations(unsigned int cptIter)
Definition: tools/PointAdvect.h:277
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:147
Vec3SGrid Vec3fGrid
Definition: openvdb.h:56
typename CptGridType::ValueType CptValueType
Definition: tools/PointAdvect.h:40
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: tools/PointAdvect.h:263
GridT GridType
Definition: tools/PointAdvect.h:249
void projectToConstraintSurface(LocationType &W) const
Definition: tools/PointAdvect.h:63
bool getThreaded()
Definition: tools/PointAdvect.h:281
ClosestPointProjector()
Definition: tools/PointAdvect.h:42
CptGridT CptGridType
Definition: tools/PointAdvect.h:38
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:278
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:153
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:250
bool earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: tools/PointAdvect.h:143
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: tools/PointAdvect.h:46
math::Vec3< Real > Vec3R
Definition: openvdb/Types.h:50
PointListT PointListType
Definition: tools/PointAdvect.h:253
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:284