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