OpenVDB  7.0.0
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
Definition: tools/PointAdvect.h:35
void setThreaded(bool threaded)
get & set
Definition: tools/PointAdvect.h:145
GridT GridType
Definition: tools/PointAdvect.h:116
Definition: tools/PointAdvect.h:246
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...
ClosestPointProjector(const ClosestPointProjector &other)
Definition: tools/PointAdvect.h:52
PointAdvect(const PointAdvect &other)
Definition: tools/PointAdvect.h:129
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
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:278
void setThreaded(bool threaded)
Definition: tools/PointAdvect.h:280
void setConstraintIterations(unsigned int cptIterations)
Definition: tools/PointAdvect.h:58
void setConstraintIterations(unsigned int cptIter)
Definition: tools/PointAdvect.h:277
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:121
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:214
virtual ~ConstrainedPointAdvect()
Definition: tools/PointAdvect.h:275
Vec3SGrid Vec3fGrid
Definition: openvdb.h:56
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=nullptr)
Definition: tools/PointAdvect.h:255
bool getThreaded()
Definition: tools/PointAdvect.h:146
typename CptGridType::ValueType CptValueType
Definition: tools/PointAdvect.h:40
ClosestPointProjector()
Definition: tools/PointAdvect.h:42
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
CptGridT CptGridType
Definition: tools/PointAdvect.h:38
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:118
PointListT PointListType
Definition: tools/PointAdvect.h:253
typename PointListT::value_type LocationType
Definition: tools/PointAdvect.h:250
Definition: Exceptions.h:13
math::Vec3< Real > Vec3R
Definition: Types.h:49
Dummy NOOP interrupter class defining interface.
Definition: NullInterrupter.h:25
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
void projectToConstraintSurface(LocationType &W) const
Definition: tools/PointAdvect.h:63
GridT GridType
Definition: tools/PointAdvect.h:249
bool getThreaded()
Definition: tools/PointAdvect.h:281
void setIntegrationOrder(unsigned int order)
Definition: tools/PointAdvect.h:147
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
unsigned int numIterations()
Definition: tools/PointAdvect.h:59
virtual ~PointAdvect()
Definition: tools/PointAdvect.h:139
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: tools/PointAdvect.h:263
typename CptGridType::ConstAccessor CptAccessor
Definition: tools/PointAdvect.h:39
PointListT PointListType
Definition: tools/PointAdvect.h:117
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
Definition: tools/PointAdvect.h:113