OpenVDB  12.0.0
PointAdvect.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: Apache-2.0
3 //
4 /// @author Ken Museth, D.J. Hill (openvdb port, added staggered grid support)
5 ///
6 /// @file tools/PointAdvect.h
7 ///
8 /// @brief Class PointAdvect advects points (with position) in a static velocity field
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/Types.h> // Vec3 types and version number
15 #include <openvdb/Grid.h> // grid
16 #include <openvdb/math/Math.h> // min
18 #include <openvdb/thread/Threading.h>
19 #include "Interpolation.h" // sampling
20 #include "VelocityFields.h" // VelocityIntegrator
21 
22 #include <tbb/blocked_range.h> // threading
23 #include <tbb/parallel_for.h> // threading
24 
25 #include <vector>
26 
27 
28 namespace openvdb {
30 namespace OPENVDB_VERSION_NAME {
31 namespace tools {
32 
33 /// Class that holds a Vec3 grid, to be interpreted as the closest point to a constraint
34 /// surface. Supports a method to allow a point to be projected onto the closest point
35 /// on the constraint surface. Uses Caching.
36 template<typename CptGridT = Vec3fGrid>
38 {
39 public:
40  using CptGridType = CptGridT;
41  using CptAccessor = typename CptGridType::ConstAccessor;
42  using CptValueType = typename CptGridType::ValueType;
43 
45  mCptIterations(0)
46  {
47  }
48  ClosestPointProjector(const CptGridType& cptGrid, int n):
49  mCptGrid(&cptGrid),
50  mCptAccessor(cptGrid.getAccessor()),
51  mCptIterations(n)
52  {
53  }
55  mCptGrid(other.mCptGrid),
56  mCptAccessor(mCptGrid->getAccessor()),
57  mCptIterations(other.mCptIterations)
58  {
59  }
60  void setConstraintIterations(unsigned int cptIterations) { mCptIterations = cptIterations; }
61  unsigned int numIterations() { return mCptIterations; }
62 
63  // point constraint
64  template <typename LocationType>
65  inline void projectToConstraintSurface(LocationType& W) const
66  {
67  /// Entries in the CPT tree are the closest point to the constraint surface.
68  /// The interpolation step in sample introduces error so that the result
69  /// of a single sample may not lie exactly on the surface. The iterations
70  /// in the loop exist to minimize this error.
71  CptValueType result(W[0], W[1],W[2]);
72  for (unsigned int i = 0; i < mCptIterations; ++i) {
73  const Vec3R location = mCptGrid->worldToIndex(Vec3R(result[0], result[1], result[2]));
74  BoxSampler::sample<CptAccessor>(mCptAccessor, location, result);
75  }
76  W[0] = result[0];
77  W[1] = result[1];
78  W[2] = result[2];
79  }
80 
81 private:
82  const CptGridType* mCptGrid; // Closest-Point-Transform vector field
83  CptAccessor mCptAccessor;
84  unsigned int mCptIterations;
85 };// end of ClosestPointProjector class
86 
87 ////////////////////////////////////////
88 
89 
90 /// Performs passive or constrained advection of points in a velocity field
91 /// represented by an OpenVDB grid and an optional closest-point-transform (CPT)
92 /// represented in another OpenVDB grid. Note the CPT is assumed to be
93 /// in world coordinates and NOT index coordinates!
94 /// Supports both collocated velocity grids and staggered velocity grids
95 ///
96 /// The @c PointListT template argument refers to any class with the following
97 /// interface (e.g., std::vector<openvdb::Vec3f>):
98 /// @code
99 /// class PointList {
100 /// ...
101 /// public:
102 /// using value_type = internal_vector3_type; // must support [] component access
103 /// openvdb::Index size() const; // number of points in list
104 /// value_type& operator[](int n); // world space position of nth point
105 /// };
106 /// @endcode
107 ///
108 /// @note All methods (except size) are assumed to be thread-safe and
109 /// the positions are returned as non-const references since the
110 /// advection method needs to modify them!
111 template<typename GridT = Vec3fGrid,
112  typename PointListT = std::vector<typename GridT::ValueType>,
113  bool StaggeredVelocity = false,
114  typename InterrupterType = util::NullInterrupter>
116 {
117 public:
118  using GridType = GridT;
119  using PointListType = PointListT;
120  using LocationType = typename PointListT::value_type;
122 
123  PointAdvect(const GridT& velGrid, InterrupterType* interrupter = nullptr):
124  mVelGrid(&velGrid),
125  mPoints(nullptr),
126  mIntegrationOrder(1),
127  mThreaded(true),
128  mInterrupter(interrupter)
129  {
130  }
131  PointAdvect(const PointAdvect &other) :
132  mVelGrid(other.mVelGrid),
133  mPoints(other.mPoints),
134  mDt(other.mDt),
135  mAdvIterations(other.mAdvIterations),
136  mIntegrationOrder(other.mIntegrationOrder),
137  mThreaded(other.mThreaded),
138  mInterrupter(other.mInterrupter)
139  {
140  }
141  virtual ~PointAdvect()
142  {
143  }
144  /// If the order of the integration is set to zero no advection is performed
145  bool earlyOut() const { return (mIntegrationOrder==0);}
146  /// get & set
147  void setThreaded(bool threaded) { mThreaded = threaded; }
148  bool getThreaded() { return mThreaded; }
149  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
150 
151  /// Constrained advection of a list of points over a time = dt * advIterations
152  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
153  {
154  if (this->earlyOut()) return; // nothing to do!
155  mPoints = &points;
156  mDt = dt;
157  mAdvIterations = advIterations;
158 
159  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
160  if (mThreaded) {
161  tbb::parallel_for(tbb::blocked_range<size_t>(0, mPoints->size()), *this);
162  } else {
163  (*this)(tbb::blocked_range<size_t>(0, mPoints->size()));
164  }
165  if (mInterrupter) mInterrupter->end();
166  }
167 
168  /// Never call this method directly - it is use by TBB and has to be public!
169  void operator() (const tbb::blocked_range<size_t> &range) const
170  {
171  if (mInterrupter && mInterrupter->wasInterrupted()) {
172  thread::cancelGroupExecution();
173  }
174 
175  VelocityFieldIntegrator velField(*mVelGrid);
176  switch (mIntegrationOrder) {
177  case 1:
178  {
179  for (size_t n = range.begin(); n != range.end(); ++n) {
180  LocationType& X0 = (*mPoints)[n];
181  // loop over number of time steps
182  for (unsigned int i = 0; i < mAdvIterations; ++i) {
183  velField.template rungeKutta<1>(mDt, X0);
184  }
185  }
186  }
187  break;
188  case 2:
189  {
190  for (size_t n = range.begin(); n != range.end(); ++n) {
191  LocationType& X0 = (*mPoints)[n];
192  // loop over number of time steps
193  for (unsigned int i = 0; i < mAdvIterations; ++i) {
194  velField.template rungeKutta<2>(mDt, X0);
195  }
196  }
197  }
198  break;
199  case 3:
200  {
201  for (size_t n = range.begin(); n != range.end(); ++n) {
202  LocationType& X0 = (*mPoints)[n];
203  // loop over number of time steps
204  for (unsigned int i = 0; i < mAdvIterations; ++i) {
205  velField.template rungeKutta<3>(mDt, X0);
206  }
207  }
208  }
209  break;
210  case 4:
211  {
212  for (size_t n = range.begin(); n != range.end(); ++n) {
213  LocationType& X0 = (*mPoints)[n];
214  // loop over number of time steps
215  for (unsigned int i = 0; i < mAdvIterations; ++i) {
216  velField.template rungeKutta<4>(mDt, X0);
217  }
218  }
219  }
220  break;
221  }
222  }
223 
224 private:
225  // the velocity field
226  const GridType* mVelGrid;
227 
228  // vertex list of all the points
229  PointListT* mPoints;
230 
231  // time integration parameters
232  float mDt; // time step
233  unsigned int mAdvIterations; // number of time steps
234  unsigned int mIntegrationOrder;
235 
236  // operational parameters
237  bool mThreaded;
238  InterrupterType* mInterrupter;
239 
240 };//end of PointAdvect class
241 
242 
243 template<typename GridT = Vec3fGrid,
244  typename PointListT = std::vector<typename GridT::ValueType>,
245  bool StaggeredVelocity = false,
246  typename CptGridType = GridT,
247  typename InterrupterType = util::NullInterrupter>
249 {
250 public:
251  using GridType = GridT;
252  using LocationType = typename PointListT::value_type;
255  using PointListType = PointListT;
256 
258  const GridType& cptGrid, int cptn, InterrupterType* interrupter = nullptr):
259  mVelGrid(&velGrid),
260  mCptGrid(&cptGrid),
261  mCptIter(cptn),
262  mInterrupter(interrupter)
263  {
264  }
266  mVelGrid(other.mVelGrid),
267  mCptGrid(other.mCptGrid),
268  mCptIter(other.mCptIter),
269  mPoints(other.mPoints),
270  mDt(other.mDt),
271  mAdvIterations(other.mAdvIterations),
272  mIntegrationOrder(other.mIntegrationOrder),
273  mThreaded(other.mThreaded),
274  mInterrupter(other.mInterrupter)
275  {
276  }
278 
279  void setConstraintIterations(unsigned int cptIter) {mCptIter = cptIter;}
280  void setIntegrationOrder(unsigned int order) {mIntegrationOrder = order;}
281 
282  void setThreaded(bool threaded) { mThreaded = threaded; }
283  bool getThreaded() { return mThreaded; }
284 
285  /// Constrained Advection a list of points over a time = dt * advIterations
286  void advect(PointListT& points, float dt, unsigned int advIterations = 1)
287  {
288  mPoints = &points;
289  mDt = dt;
290 
291  if (mIntegrationOrder==0 && mCptIter == 0) {
292  return; // nothing to do!
293  }
294  (mIntegrationOrder>0) ? mAdvIterations = advIterations : mAdvIterations = 1;
295 
296  if (mInterrupter) mInterrupter->start("Advecting points by OpenVDB velocity field: ");
297  const size_t N = mPoints->size();
298 
299  if (mThreaded) {
300  tbb::parallel_for(tbb::blocked_range<size_t>(0, N), *this);
301  } else {
302  (*this)(tbb::blocked_range<size_t>(0, N));
303  }
304  if (mInterrupter) mInterrupter->end();
305  }
306 
307 
308  /// Never call this method directly - it is use by TBB and has to be public!
309  void operator() (const tbb::blocked_range<size_t> &range) const
310  {
311  if (mInterrupter && mInterrupter->wasInterrupted()) {
312  thread::cancelGroupExecution();
313  }
314 
315  VelocityIntegratorType velField(*mVelGrid);
316  ClosestPointProjectorType cptField(*mCptGrid, mCptIter);
317  switch (mIntegrationOrder) {
318  case 0://pure CPT projection
319  {
320  for (size_t n = range.begin(); n != range.end(); ++n) {
321  LocationType& X0 = (*mPoints)[n];
322  for (unsigned int i = 0; i < mAdvIterations; ++i) {
323  cptField.projectToConstraintSurface(X0);
324  }
325  }
326  }
327  break;
328  case 1://1'th order advection and CPT projection
329  {
330  for (size_t n = range.begin(); n != range.end(); ++n) {
331  LocationType& X0 = (*mPoints)[n];
332  for (unsigned int i = 0; i < mAdvIterations; ++i) {
333  velField.template rungeKutta<1>(mDt, X0);
334  cptField.projectToConstraintSurface(X0);
335  }
336  }
337  }
338  break;
339  case 2://2'nd order advection and CPT projection
340  {
341  for (size_t n = range.begin(); n != range.end(); ++n) {
342  LocationType& X0 = (*mPoints)[n];
343  for (unsigned int i = 0; i < mAdvIterations; ++i) {
344  velField.template rungeKutta<2>(mDt, X0);
345  cptField.projectToConstraintSurface(X0);
346  }
347  }
348  }
349  break;
350 
351  case 3://3'rd order advection and CPT projection
352  {
353  for (size_t n = range.begin(); n != range.end(); ++n) {
354  LocationType& X0 = (*mPoints)[n];
355  for (unsigned int i = 0; i < mAdvIterations; ++i) {
356  velField.template rungeKutta<3>(mDt, X0);
357  cptField.projectToConstraintSurface(X0);
358  }
359  }
360  }
361  break;
362  case 4://4'th order advection and CPT projection
363  {
364  for (size_t n = range.begin(); n != range.end(); ++n) {
365  LocationType& X0 = (*mPoints)[n];
366  for (unsigned int i = 0; i < mAdvIterations; ++i) {
367  velField.template rungeKutta<4>(mDt, X0);
368  cptField.projectToConstraintSurface(X0);
369  }
370  }
371  }
372  break;
373  }
374  }
375 
376 private:
377  const GridType* mVelGrid; // the velocity field
378  const GridType* mCptGrid;
379  int mCptIter;
380  PointListT* mPoints; // vertex list of all the points
381 
382  // time integration parameters
383  float mDt; // time step
384  unsigned int mAdvIterations; // number of time steps
385  unsigned int mIntegrationOrder; // order of Runge-Kutta integration
386  // operational parameters
387  bool mThreaded;
388  InterrupterType* mInterrupter;
389 };// end of ConstrainedPointAdvect class
390 
391 } // namespace tools
392 } // namespace OPENVDB_VERSION_NAME
393 } // namespace openvdb
394 
395 #endif // OPENVDB_TOOLS_POINT_ADVECT_HAS_BEEN_INCLUDED
typename CptGridType::ValueType CptValueType
Definition: PointAdvect.h:42
void projectToConstraintSurface(LocationType &W) const
Definition: PointAdvect.h:65
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: PointAdvect.h:119
void setConstraintIterations(unsigned int cptIter)
Definition: PointAdvect.h:279
void advect(PointListT &points, float dt, unsigned int advIterations=1)
Constrained Advection a list of points over a time = dt * advIterations.
Definition: PointAdvect.h:286
Base class for interrupters.
Definition: NullInterrupter.h:25
GridT GridType
Definition: PointAdvect.h:251
Performs Runge-Kutta time integration of variable order in a static velocity field.
Definition: VelocityFields.h:215
PointAdvect(const GridT &velGrid, InterrupterType *interrupter=nullptr)
Definition: PointAdvect.h:123
virtual ~PointAdvect()
Definition: PointAdvect.h:141
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:152
typename PointListT::value_type LocationType
Definition: PointAdvect.h:252
ConstrainedPointAdvect(const GridType &velGrid, const GridType &cptGrid, int cptn, InterrupterType *interrupter=nullptr)
Definition: PointAdvect.h:257
Definition: PointAdvect.h:115
ClosestPointProjector(const ClosestPointProjector &other)
Definition: PointAdvect.h:54
OutGridT XformOp bool threaded
Definition: ValueTransformer.h:140
ConstrainedPointAdvect(const ConstrainedPointAdvect &other)
Definition: PointAdvect.h:265
bool getThreaded()
Definition: PointAdvect.h:148
Definition: Exceptions.h:13
PointAdvect(const PointAdvect &other)
Definition: PointAdvect.h:131
CptGridT CptGridType
Definition: PointAdvect.h:40
bool getThreaded()
Definition: PointAdvect.h:283
virtual ~ConstrainedPointAdvect()
Definition: PointAdvect.h:277
ClosestPointProjector()
Definition: PointAdvect.h:44
typename PointListT::value_type LocationType
Definition: PointAdvect.h:120
void setThreaded(bool threaded)
get & set
Definition: PointAdvect.h:147
bool earlyOut() const
If the order of the integration is set to zero no advection is performed.
Definition: PointAdvect.h:145
void setConstraintIterations(unsigned int cptIterations)
Definition: PointAdvect.h:60
void setIntegrationOrder(unsigned int order)
Definition: PointAdvect.h:149
Vec3SGrid Vec3fGrid
Definition: openvdb.h:85
PointListT PointListType
Definition: PointAdvect.h:255
typename CptGridType::ConstAccessor CptAccessor
Definition: PointAdvect.h:41
void setIntegrationOrder(unsigned int order)
Definition: PointAdvect.h:280
GridT GridType
Definition: PointAdvect.h:118
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:121
ClosestPointProjector(const CptGridType &cptGrid, int n)
Definition: PointAdvect.h:48
void setThreaded(bool threaded)
Definition: PointAdvect.h:282
math::Vec3< Real > Vec3R
Definition: Types.h:72
unsigned int numIterations()
Definition: PointAdvect.h:61
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:218