OpenVDB  9.1.1
RayIntersector.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 ///
4 /// @file RayIntersector.h
5 ///
6 /// @author Ken Museth
7 ///
8 /// @brief Accelerated intersection of a ray with a narrow-band level
9 /// set or a generic (e.g. density) volume. This will of course be
10 /// useful for respectively surface and volume rendering.
11 ///
12 /// @details This file defines two main classes,
13 /// LevelSetRayIntersector and VolumeRayIntersector, as well as the
14 /// three support classes LevelSetHDDA, VolumeHDDA and LinearSearchImpl.
15 /// The LevelSetRayIntersector is templated on the LinearSearchImpl class
16 /// and calls instances of the LevelSetHDDA class. The reason to split
17 /// level set ray intersection into three classes is twofold. First
18 /// LevelSetRayIntersector defines the public API for client code and
19 /// LinearSearchImpl defines the actual algorithm used for the
20 /// ray level-set intersection. In other words this design will allow
21 /// for the public API to be fixed while the intersection algorithm
22 /// can change without resolving to (slow) virtual methods. Second,
23 /// LevelSetHDDA, which implements a hierarchical Differential Digital
24 /// Analyzer, relies on partial template specialization, so it has to
25 /// be a standalone class (as opposed to a member class of
26 /// LevelSetRayIntersector). The VolumeRayIntersector is conceptually
27 /// much simpler than the LevelSetRayIntersector, and hence it only
28 /// depends on VolumeHDDA that implements the hierarchical
29 /// Differential Digital Analyzer.
30 
31 
32 #ifndef OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
33 #define OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
34 
35 #include <openvdb/math/DDA.h>
36 #include <openvdb/math/Math.h>
37 #include <openvdb/math/Ray.h>
38 #include <openvdb/math/Stencils.h>
39 #include <openvdb/Grid.h>
40 #include <openvdb/Types.h>
41 #include "Morphology.h"
42 #include <iostream>
43 #include <type_traits>
44 
45 
46 namespace openvdb {
48 namespace OPENVDB_VERSION_NAME {
49 namespace tools {
50 
51 // Helper class that implements the actual search of the zero-crossing
52 // of the level set along the direction of a ray. This particular
53 // implementation uses iterative linear search.
54 template<typename GridT, int Iterations = 0, typename RealT = double>
56 
57 
58 ///////////////////////////////////// LevelSetRayIntersector /////////////////////////////////////
59 
60 
61 /// @brief This class provides the public API for intersecting a ray
62 /// with a narrow-band level set.
63 ///
64 /// @details It wraps a SearchImplT with a simple public API and
65 /// performs the actual hierarchical tree node and voxel traversal.
66 ///
67 /// @warning Use the (default) copy-constructor to make sure each
68 /// computational thread has their own instance of this class. This is
69 /// important since the SearchImplT contains a ValueAccessor that is
70 /// not thread-safe. However copying is very efficient.
71 ///
72 /// @see tools/RayTracer.h for examples of intended usage.
73 ///
74 /// @todo Add TrilinearSearchImpl, as an alternative to LinearSearchImpl,
75 /// that performs analytical 3D trilinear intersection tests, i.e., solves
76 /// cubic equations. This is slower but also more accurate than the 1D
77 /// linear interpolation in LinearSearchImpl.
78 template<typename GridT,
79  typename SearchImplT = LinearSearchImpl<GridT>,
80  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
81  typename RayT = math::Ray<Real> >
83 {
84 public:
85  using GridType = GridT;
86  using RayType = RayT;
87  using RealType = typename RayT::RealType;
88  using Vec3Type = typename RayT::Vec3T;
89  using ValueT = typename GridT::ValueType;
90  using TreeT = typename GridT::TreeType;
91 
92  static_assert(NodeLevel >= -1 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range");
94  "level set grids must have scalar, floating-point value types");
95 
96  /// @brief Constructor
97  /// @param grid level set grid to intersect rays against.
98  /// @param isoValue optional iso-value for the ray-intersection.
99  LevelSetRayIntersector(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>())
100  : mTester(grid, isoValue)
101  {
102  if (!grid.hasUniformVoxels() ) {
104  "LevelSetRayIntersector only supports uniform voxels!");
105  }
106  if (grid.getGridClass() != GRID_LEVEL_SET) {
108  "LevelSetRayIntersector only supports level sets!"
109  "\nUse Grid::setGridClass(openvdb::GRID_LEVEL_SET)");
110  }
111  }
112 
113  /// @brief Return the iso-value used for ray-intersections
114  const ValueT& getIsoValue() const { return mTester.getIsoValue(); }
115 
116  /// @brief Return @c true if the index-space ray intersects the level set.
117  /// @param iRay ray represented in index space.
118  bool intersectsIS(const RayType& iRay) const
119  {
120  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
122  }
123 
124  /// @brief Return @c true if the index-space ray intersects the level set
125  /// @param iRay ray represented in index space.
126  /// @param iTime if an intersection was found it is assigned the time of the
127  /// intersection along the index ray.
128  bool intersectsIS(const RayType& iRay, RealType &iTime) const
129  {
130  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
131  iTime = mTester.getIndexTime();
133  }
134 
135  /// @brief Return @c true if the index-space ray intersects the level set.
136  /// @param iRay ray represented in index space.
137  /// @param xyz if an intersection was found it is assigned the
138  /// intersection point in index space, otherwise it is unchanged.
139  bool intersectsIS(const RayType& iRay, Vec3Type& xyz) const
140  {
141  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
142  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
143  mTester.getIndexPos(xyz);
144  return true;
145  }
146 
147  /// @brief Return @c true if the index-space ray intersects the level set.
148  /// @param iRay ray represented in index space.
149  /// @param xyz if an intersection was found it is assigned the
150  /// intersection point in index space, otherwise it is unchanged.
151  /// @param iTime if an intersection was found it is assigned the time of the
152  /// intersection along the index ray.
153  bool intersectsIS(const RayType& iRay, Vec3Type& xyz, RealType &iTime) const
154  {
155  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
156  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
157  mTester.getIndexPos(xyz);
158  iTime = mTester.getIndexTime();
159  return true;
160  }
161 
162  /// @brief Return @c true if the world-space ray intersects the level set.
163  /// @param wRay ray represented in world space.
164  bool intersectsWS(const RayType& wRay) const
165  {
166  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
168  }
169 
170  /// @brief Return @c true if the world-space ray intersects the level set.
171  /// @param wRay ray represented in world space.
172  /// @param wTime if an intersection was found it is assigned the time of the
173  /// intersection along the world ray.
174  bool intersectsWS(const RayType& wRay, RealType &wTime) const
175  {
176  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
177  wTime = mTester.getWorldTime();
179  }
180 
181  /// @brief Return @c true if the world-space ray intersects the level set.
182  /// @param wRay ray represented in world space.
183  /// @param world if an intersection was found it is assigned the
184  /// intersection point in world space, otherwise it is unchanged
185  bool intersectsWS(const RayType& wRay, Vec3Type& world) const
186  {
187  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
188  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
189  mTester.getWorldPos(world);
190  return true;
191  }
192 
193  /// @brief Return @c true if the world-space ray intersects the level set.
194  /// @param wRay ray represented in world space.
195  /// @param world if an intersection was found it is assigned the
196  /// intersection point in world space, otherwise it is unchanged.
197  /// @param wTime if an intersection was found it is assigned the time of the
198  /// intersection along the world ray.
199  bool intersectsWS(const RayType& wRay, Vec3Type& world, RealType &wTime) const
200  {
201  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
202  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
203  mTester.getWorldPos(world);
204  wTime = mTester.getWorldTime();
205  return true;
206  }
207 
208  /// @brief Return @c true if the world-space ray intersects the level set.
209  /// @param wRay ray represented in world space.
210  /// @param world if an intersection was found it is assigned the
211  /// intersection point in world space, otherwise it is unchanged.
212  /// @param normal if an intersection was found it is assigned the normal
213  /// of the level set surface in world space, otherwise it is unchanged.
214  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal) const
215  {
216  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
217  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
218  mTester.getWorldPosAndNml(world, normal);
219  return true;
220  }
221 
222  /// @brief Return @c true if the world-space ray intersects the level set.
223  /// @param wRay ray represented in world space.
224  /// @param world if an intersection was found it is assigned the
225  /// intersection point in world space, otherwise it is unchanged.
226  /// @param normal if an intersection was found it is assigned the normal
227  /// of the level set surface in world space, otherwise it is unchanged.
228  /// @param wTime if an intersection was found it is assigned the time of the
229  /// intersection along the world ray.
230  bool intersectsWS(const RayType& wRay, Vec3Type& world, Vec3Type& normal, RealType &wTime) const
231  {
232  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
233  if (!math::LevelSetHDDA<TreeT, NodeLevel>::test(mTester)) return false;//missed level set
234  mTester.getWorldPosAndNml(world, normal);
235  wTime = mTester.getWorldTime();
236  return true;
237  }
238 
239 private:
240 
241  mutable SearchImplT mTester;
242 
243 };// LevelSetRayIntersector
244 
245 
246 ////////////////////////////////////// VolumeRayIntersector //////////////////////////////////////
247 
248 
249 /// @brief This class provides the public API for intersecting a ray
250 /// with a generic (e.g. density) volume.
251 /// @details Internally it performs the actual hierarchical tree node traversal.
252 /// @warning Use the (default) copy-constructor to make sure each
253 /// computational thread has their own instance of this class. This is
254 /// important since it contains a ValueAccessor that is
255 /// not thread-safe and a CoordBBox of the active voxels that should
256 /// not be re-computed for each thread. However copying is very efficient.
257 /// @par Example:
258 /// @code
259 /// // Create an instance for the master thread
260 /// VolumeRayIntersector inter(grid);
261 /// // For each additional thread use the copy constructor. This
262 /// // amortizes the overhead of computing the bbox of the active voxels!
263 /// VolumeRayIntersector inter2(inter);
264 /// // Before each ray-traversal set the index ray.
265 /// iter.setIndexRay(ray);
266 /// // or world ray
267 /// iter.setWorldRay(ray);
268 /// // Now you can begin the ray-marching using consecutive calls to VolumeRayIntersector::march
269 /// double t0=0, t1=0;// note the entry and exit times are with respect to the INDEX ray
270 /// while ( inter.march(t0, t1) ) {
271 /// // perform line-integration between t0 and t1
272 /// }}
273 /// @endcode
274 template<typename GridT,
275  int NodeLevel = GridT::TreeType::RootNodeType::ChildNodeType::LEVEL,
276  typename RayT = math::Ray<Real> >
278 {
279 public:
280  using GridType = GridT;
281  using RayType = RayT;
282  using RealType = typename RayT::RealType;
283  using RootType = typename GridT::TreeType::RootNodeType;
285 
286  static_assert(NodeLevel >= 0 && NodeLevel < int(TreeT::DEPTH)-1, "NodeLevel out of range");
287 
288  /// @brief Grid constructor
289  /// @param grid Generic grid to intersect rays against.
290  /// @param dilationCount The number of voxel dilations performed
291  /// on (a boolean copy of) the input grid. This allows the
292  /// intersector to account for the size of interpolation kernels
293  /// in client code.
294  /// @throw RuntimeError if the voxels of the grid are not uniform
295  /// or the grid is empty.
296  VolumeRayIntersector(const GridT& grid, int dilationCount = 0)
297  : mIsMaster(true)
298  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
299  , mGrid(&grid)
300  , mAccessor(*mTree)
301  {
302  if (!grid.hasUniformVoxels() ) {
304  "VolumeRayIntersector only supports uniform voxels!");
305  }
306  if ( grid.empty() ) {
307  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
308  }
309 
310  // Dilate active voxels to better account for the size of interpolation kernels
312 
313  mTree->root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
314 
315  mBBox.max().offset(1);//padding so the bbox of a node becomes (origin,origin + node_dim)
316  }
317 
318  /// @brief Grid and BBox constructor
319  /// @param grid Generic grid to intersect rays against.
320  /// @param bbox The axis-aligned bounding-box in the index space of the grid.
321  /// @warning It is assumed that bbox = (min, min + dim) where min denotes
322  /// to the smallest grid coordinates and dim are the integer length of the bbox.
323  /// @throw RuntimeError if the voxels of the grid are not uniform
324  /// or the grid is empty.
325  VolumeRayIntersector(const GridT& grid, const math::CoordBBox& bbox)
326  : mIsMaster(true)
327  , mTree(new TreeT(grid.tree(), false, TopologyCopy()))
328  , mGrid(&grid)
329  , mAccessor(*mTree)
330  , mBBox(bbox)
331  {
332  if (!grid.hasUniformVoxels() ) {
334  "VolumeRayIntersector only supports uniform voxels!");
335  }
336  if ( grid.empty() ) {
337  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
338  }
339  }
340 
341  /// @brief Shallow copy constructor
342  /// @warning This copy constructor creates shallow copies of data
343  /// members of the instance passed as the argument. For
344  /// performance reasons we are not using shared pointers (their
345  /// mutex-lock impairs multi-threading).
347  : mIsMaster(false)
348  , mTree(other.mTree)//shallow copy
349  , mGrid(other.mGrid)//shallow copy
350  , mAccessor(*mTree)//initialize new (vs deep copy)
351  , mRay(other.mRay)//deep copy
352  , mTmax(other.mTmax)//deep copy
353  , mBBox(other.mBBox)//deep copy
354  {
355  }
356 
357  /// @brief Destructor
358  ~VolumeRayIntersector() { if (mIsMaster) delete mTree; }
359 
360  /// @brief Return @c false if the index ray misses the bbox of the grid.
361  /// @param iRay Ray represented in index space.
362  /// @warning Call this method (or setWorldRay) before the ray
363  /// traversal starts and use the return value to decide if further
364  /// marching is required.
365  inline bool setIndexRay(const RayT& iRay)
366  {
367  mRay = iRay;
368  const bool hit = mRay.clip(mBBox);
369  if (hit) mTmax = mRay.t1();
370  return hit;
371  }
372 
373  /// @brief Return @c false if the world ray misses the bbox of the grid.
374  /// @param wRay Ray represented in world space.
375  /// @warning Call this method (or setIndexRay) before the ray
376  /// traversal starts and use the return value to decide if further
377  /// marching is required.
378  /// @details Since hit times are computed with respect to the ray
379  /// represented in index space of the current grid, it is
380  /// recommended that either the client code uses getIndexPos to
381  /// compute index position from hit times or alternatively keeps
382  /// an instance of the index ray and instead uses setIndexRay to
383  /// initialize the ray.
384  inline bool setWorldRay(const RayT& wRay)
385  {
386  return this->setIndexRay(wRay.worldToIndex(*mGrid));
387  }
388 
389  inline typename RayT::TimeSpan march()
390  {
391  const typename RayT::TimeSpan t = mHDDA.march(mRay, mAccessor);
392  if (t.t1>0) mRay.setTimes(t.t1 + math::Delta<RealType>::value(), mTmax);
393  return t;
394  }
395 
396  /// @brief Return @c true if the ray intersects active values,
397  /// i.e. either active voxels or tiles. Only when a hit is
398  /// detected are t0 and t1 updated with the corresponding entry
399  /// and exit times along the INDEX ray!
400  /// @note Note that t0 and t1 are only resolved at the node level
401  /// (e.g. a LeafNode with active voxels) as opposed to the individual
402  /// active voxels.
403  /// @param t0 If the return value > 0 this is the time of the
404  /// first hit of an active tile or leaf.
405  /// @param t1 If the return value > t0 this is the time of the
406  /// first hit (> t0) of an inactive tile or exit point of the
407  /// BBOX for the leaf nodes.
408  /// @warning t0 and t1 are computed with respect to the ray represented in
409  /// index space of the current grid, not world space!
410  inline bool march(RealType& t0, RealType& t1)
411  {
412  const typename RayT::TimeSpan t = this->march();
413  t.get(t0, t1);
414  return t.valid();
415  }
416 
417  /// @brief Generates a list of hits along the ray.
418  ///
419  /// @param list List of hits represented as time spans.
420  ///
421  /// @note ListType is a list of RayType::TimeSpan and is required to
422  /// have the two methods: clear() and push_back(). Thus, it could
423  /// be std::vector<typename RayType::TimeSpan> or
424  /// std::deque<typename RayType::TimeSpan>.
425  template <typename ListType>
426  inline void hits(ListType& list)
427  {
428  mHDDA.hits(mRay, mAccessor, list);
429  }
430 
431  /// @brief Return the floating-point index position along the
432  /// current index ray at the specified time.
433  inline Vec3R getIndexPos(RealType time) const { return mRay(time); }
434 
435  /// @brief Return the floating-point world position along the
436  /// current index ray at the specified time.
437  inline Vec3R getWorldPos(RealType time) const { return mGrid->indexToWorld(mRay(time)); }
438 
439  inline RealType getWorldTime(RealType time) const
440  {
441  return time*mGrid->transform().baseMap()->applyJacobian(mRay.dir()).length();
442  }
443 
444  /// @brief Return a const reference to the input grid.
445  const GridT& grid() const { return *mGrid; }
446 
447  /// @brief Return a const reference to the (potentially dilated)
448  /// bool tree used to accelerate the ray marching.
449  const TreeT& tree() const { return *mTree; }
450 
451  /// @brief Return a const reference to the BBOX of the grid
452  const math::CoordBBox& bbox() const { return mBBox; }
453 
454  /// @brief Print bbox, statistics, memory usage and other information.
455  /// @param os a stream to which to write textual information
456  /// @param verboseLevel 1: print bbox only; 2: include boolean tree
457  /// statistics; 3: include memory usage
458  void print(std::ostream& os = std::cout, int verboseLevel = 1)
459  {
460  if (verboseLevel>0) {
461  os << "BBox: " << mBBox << std::endl;
462  if (verboseLevel==2) {
463  mTree->print(os, 1);
464  } else if (verboseLevel>2) {
465  mTree->print(os, 2);
466  }
467  }
468  }
469 
470 private:
471  using AccessorT = typename tree::ValueAccessor<const TreeT,/*IsSafe=*/false>;
472 
473  const bool mIsMaster;
474  TreeT* mTree;
475  const GridT* mGrid;
476  AccessorT mAccessor;
477  RayT mRay;
478  RealType mTmax;
479  math::CoordBBox mBBox;
481 
482 };// VolumeRayIntersector
483 
484 
485 //////////////////////////////////////// LinearSearchImpl ////////////////////////////////////////
486 
487 
488 /// @brief Implements linear iterative search for an iso-value of
489 /// the level set along the direction of the ray.
490 ///
491 /// @note Since this class is used internally in
492 /// LevelSetRayIntersector (define above) and LevelSetHDDA (defined below)
493 /// client code should never interact directly with its API. This also
494 /// explains why we are not concerned with the fact that several of
495 /// its methods are unsafe to call unless roots were already detected.
496 ///
497 /// @details It is approximate due to the limited number of iterations
498 /// which can can be defined with a template parameter. However the default value
499 /// has proven surprisingly accurate and fast. In fact more iterations
500 /// are not guaranteed to give significantly better results.
501 ///
502 /// @warning Since the root-searching algorithm is approximate
503 /// (first-order) it is possible to miss intersections if the
504 /// iso-value is too close to the inside or outside of the narrow
505 /// band (typically a distance less than a voxel unit).
506 ///
507 /// @warning Since this class internally stores a ValueAccessor it is NOT thread-safe,
508 /// so make sure to give each thread its own instance. This of course also means that
509 /// the cost of allocating an instance should (if possible) be amortized over
510 /// as many ray intersections as possible.
511 template<typename GridT, int Iterations, typename RealT>
512 class LinearSearchImpl
513 {
514 public:
515  using RayT = math::Ray<RealT>;
517  using ValueT = typename GridT::ValueType;
518  using AccessorT = typename GridT::ConstAccessor;
520 
521  /// @brief Constructor from a grid.
522  /// @throw RunTimeError if the grid is empty.
523  /// @throw ValueError if the isoValue is not inside the narrow-band.
524  LinearSearchImpl(const GridT& grid, const ValueT& isoValue = zeroVal<ValueT>())
525  : mStencil(grid),
526  mIsoValue(isoValue),
527  mMinValue(isoValue - ValueT(2 * grid.voxelSize()[0])),
528  mMaxValue(isoValue + ValueT(2 * grid.voxelSize()[0]))
529  {
530  if ( grid.empty() ) {
531  OPENVDB_THROW(RuntimeError, "LinearSearchImpl does not supports empty grids");
532  }
533  if (mIsoValue<= -grid.background() ||
534  mIsoValue>= grid.background() ){
535  OPENVDB_THROW(ValueError, "The iso-value must be inside the narrow-band!");
536  }
537  grid.tree().root().evalActiveBoundingBox(mBBox, /*visit individual voxels*/false);
538  }
539 
540  /// @brief Return the iso-value used for ray-intersections
541  const ValueT& getIsoValue() const { return mIsoValue; }
542 
543  /// @brief Return @c false if the ray misses the bbox of the grid.
544  /// @param iRay Ray represented in index space.
545  /// @warning Call this method before the ray traversal starts.
546  inline bool setIndexRay(const RayT& iRay)
547  {
548  mRay = iRay;
549  return mRay.clip(mBBox);//did it hit the bbox
550  }
551 
552  /// @brief Return @c false if the ray misses the bbox of the grid.
553  /// @param wRay Ray represented in world space.
554  /// @warning Call this method before the ray traversal starts.
555  inline bool setWorldRay(const RayT& wRay)
556  {
557  mRay = wRay.worldToIndex(mStencil.grid());
558  return mRay.clip(mBBox);//did it hit the bbox
559  }
560 
561  /// @brief Get the intersection point in index space.
562  /// @param xyz The position in index space of the intersection.
563  inline void getIndexPos(VecT& xyz) const { xyz = mRay(mTime); }
564 
565  /// @brief Get the intersection point in world space.
566  /// @param xyz The position in world space of the intersection.
567  inline void getWorldPos(VecT& xyz) const { xyz = mStencil.grid().indexToWorld(mRay(mTime)); }
568 
569  /// @brief Get the intersection point and normal in world space
570  /// @param xyz The position in world space of the intersection.
571  /// @param nml The surface normal in world space of the intersection.
572  inline void getWorldPosAndNml(VecT& xyz, VecT& nml)
573  {
574  this->getIndexPos(xyz);
575  mStencil.moveTo(xyz);
576  nml = mStencil.gradient(xyz);
577  nml.normalize();
578  xyz = mStencil.grid().indexToWorld(xyz);
579  }
580 
581  /// @brief Return the time of intersection along the index ray.
582  inline RealT getIndexTime() const { return mTime; }
583 
584  /// @brief Return the time of intersection along the world ray.
585  inline RealT getWorldTime() const
586  {
587  return mTime*mStencil.grid().transform().baseMap()->applyJacobian(mRay.dir()).length();
588  }
589 
590 private:
591 
592  /// @brief Initiate the local voxel intersection test.
593  /// @warning Make sure to call this method before the local voxel intersection test.
594  inline void init(RealT t0)
595  {
596  mT[0] = t0;
597  mV[0] = static_cast<ValueT>(this->interpValue(t0));
598  }
599 
600  inline void setRange(RealT t0, RealT t1) { mRay.setTimes(t0, t1); }
601 
602  /// @brief Return a const reference to the ray.
603  inline const RayT& ray() const { return mRay; }
604 
605  /// @brief Return true if a node of the specified type exists at ijk.
606  template <typename NodeT>
607  inline bool hasNode(const Coord& ijk)
608  {
609  return mStencil.accessor().template probeConstNode<NodeT>(ijk) != nullptr;
610  }
611 
612  /// @brief Return @c true if an intersection is detected.
613  /// @param ijk Grid coordinate of the node origin or voxel being tested.
614  /// @param time Time along the index ray being tested.
615  /// @warning Only if an intersection is detected is it safe to
616  /// call getIndexPos, getWorldPos and getWorldPosAndNml!
617  inline bool operator()(const Coord& ijk, RealT time)
618  {
619  ValueT V;
620  if (mStencil.accessor().probeValue(ijk, V) &&//within narrow band
621  V>mMinValue && V<mMaxValue) {// and close to iso-value?
622  mT[1] = time;
623  mV[1] = static_cast<ValueT>(this->interpValue(time));
624  if (math::ZeroCrossing(mV[0], mV[1])) {
625  mTime = this->interpTime();
627  for (int n=0; Iterations>0 && n<Iterations; ++n) {//resolved at compile-time
628  V = static_cast<ValueT>(this->interpValue(mTime));
629  const int m = math::ZeroCrossing(mV[0], V) ? 1 : 0;
630  mV[m] = V;
631  mT[m] = mTime;
632  mTime = this->interpTime();
633  }
635  return true;
636  }
637  mT[0] = mT[1];
638  mV[0] = mV[1];
639  }
640  return false;
641  }
642 
643  inline RealT interpTime()
644  {
645  assert( math::isApproxLarger(mT[1], mT[0], RealT(1e-6) ) );
646  return mT[0]+(mT[1]-mT[0])*mV[0]/(mV[0]-mV[1]);
647  }
648 
649  inline RealT interpValue(RealT time)
650  {
651  const VecT pos = mRay(time);
652  mStencil.moveTo(pos);
653  return mStencil.interpolation(pos) - mIsoValue;
654  }
655 
656  template<typename, int> friend struct math::LevelSetHDDA;
657 
658  RayT mRay;
659  StencilT mStencil;
660  RealT mTime;//time of intersection
661  ValueT mV[2];
662  RealT mT[2];
663  const ValueT mIsoValue, mMinValue, mMaxValue;
664  math::CoordBBox mBBox;
665 };// LinearSearchImpl
666 
667 } // namespace tools
668 } // namespace OPENVDB_VERSION_NAME
669 } // namespace openvdb
670 
671 #endif // OPENVDB_TOOLS_RAYINTERSECTOR_HAS_BEEN_INCLUDED
bool setWorldRay(const RayT &wRay)
Return false if the ray misses the bbox of the grid.
Definition: RayIntersector.h:555
RayT RayType
Definition: RayIntersector.h:281
RayT::TimeSpan march()
Definition: RayIntersector.h:389
typename RayT::RealType RealType
Definition: RayIntersector.h:282
Vec3R getIndexPos(RealType time) const
Return the floating-point index position along the current index ray at the specified time...
Definition: RayIntersector.h:433
This class provides the public API for intersecting a ray with a narrow-band level set...
Definition: RayIntersector.h:82
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
Definition: RayIntersector.h:114
RealT getWorldTime() const
Return the time of intersection along the world ray.
Definition: RayIntersector.h:585
bool intersectsWS(const RayType &wRay, Vec3Type &world) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:185
const math::CoordBBox & bbox() const
Return a const reference to the BBOX of the grid.
Definition: RayIntersector.h:452
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
typename RayT::Vec3T Vec3Type
Definition: RayIntersector.h:88
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
VolumeRayIntersector(const VolumeRayIntersector &other)
Shallow copy constructor.
Definition: RayIntersector.h:346
bool intersectsIS(const RayType &iRay) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:118
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:214
Definition: Tree.h:177
void getWorldPosAndNml(VecT &xyz, VecT &nml)
Get the intersection point and normal in world space.
Definition: RayIntersector.h:572
Definition: Morphology.h:80
RealT getIndexTime() const
Return the time of intersection along the index ray.
Definition: RayIntersector.h:582
LevelSetRayIntersector(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor.
Definition: RayIntersector.h:99
Definition: Exceptions.h:63
void getIndexPos(VecT &xyz) const
Get the intersection point in index space.
Definition: RayIntersector.h:563
void dilateActiveValues(TreeOrLeafManagerT &tree, const int iterations=1, const NearestNeighbors nn=NN_FACE, const TilePolicy mode=PRESERVE_TILES, const bool threaded=true)
Topologically dilate all active values (i.e. both voxels and tiles) in a tree using one of three near...
Definition: Morphology.h:1055
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:248
bool setWorldRay(const RayT &wRay)
Return false if the world ray misses the bbox of the grid.
Definition: RayIntersector.h:384
RealType getWorldTime(RealType time) const
Definition: RayIntersector.h:439
bool setIndexRay(const RayT &iRay)
Return false if the ray misses the bbox of the grid.
Definition: RayIntersector.h:546
Definition: Exceptions.h:65
void hits(ListType &list)
Generates a list of hits along the ray.
Definition: RayIntersector.h:426
~VolumeRayIntersector()
Destructor.
Definition: RayIntersector.h:358
GridT GridType
Definition: RayIntersector.h:280
bool intersectsWS(const RayType &wRay, RealType &wTime) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:174
Definition: Types.h:416
VolumeRayIntersector(const GridT &grid, int dilationCount=0)
Grid constructor.
Definition: RayIntersector.h:296
typename GridT::ConstAccessor AccessorT
Definition: RayIntersector.h:518
This class provides the public API for intersecting a ray with a generic (e.g. density) volume...
Definition: RayIntersector.h:277
Vec3R getWorldPos(RealType time) const
Return the floating-point world position along the current index ray at the specified time...
Definition: RayIntersector.h:437
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:139
bool isApproxLarger(const Type &a, const Type &b, const Type &tolerance)
Return true if a is larger than b to within the given tolerance, i.e., if b - a < tolerance...
Definition: Math.h:434
Implements linear iterative search for an iso-value of the level set along the direction of the ray...
Definition: RayIntersector.h:55
Definition: ValueAccessor.h:182
const GridT & grid() const
Return a const reference to the input grid.
Definition: RayIntersector.h:445
LinearSearchImpl(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor from a grid.
Definition: RayIntersector.h:524
typename RayT::RealType RealType
Definition: RayIntersector.h:87
bool intersectsIS(const RayType &iRay, RealType &iTime) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:128
typename GridT::TreeType::RootNodeType RootType
Definition: RayIntersector.h:283
Implementation of morphological dilation and erosion.
RayT RayType
Definition: RayIntersector.h:86
bool intersectsIS(const RayType &iRay, Vec3Type &xyz, RealType &iTime) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:153
Definition: Exceptions.h:13
bool intersectsWS(const RayType &wRay, Vec3Type &world, RealType &wTime) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:199
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
SIMD Intrinsic Headers.
Definition: Platform.h:115
ValueT value
Definition: GridBuilder.h:1287
Delta for small floating-point offsets.
Definition: Math.h:155
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal, RealType &wTime) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:230
Definition: Morphology.h:58
bool clip(const Vec3T &center, RealT radius)
Return true if this ray intersects the specified sphere.
Definition: Ray.h:217
typename GridT::ValueType ValueT
Definition: RayIntersector.h:89
VolumeRayIntersector(const GridT &grid, const math::CoordBBox &bbox)
Grid and BBox constructor.
Definition: RayIntersector.h:325
const TreeT & tree() const
Return a const reference to the (potentially dilated) bool tree used to accelerate the ray marching...
Definition: RayIntersector.h:449
void print(std::ostream &os=std::cout, int verboseLevel=1)
Print bbox, statistics, memory usage and other information.
Definition: RayIntersector.h:458
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
Definition: RayIntersector.h:541
bool intersectsWS(const RayType &wRay) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:164
bool setIndexRay(const RayT &iRay)
Return false if the index ray misses the bbox of the grid.
Definition: RayIntersector.h:365
Helper class that implements Hierarchical Digital Differential Analyzers and is specialized for ray i...
Definition: DDA.h:144
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:116
GridT GridType
Definition: RayIntersector.h:85
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:366
bool ZeroCrossing(const Type &a, const Type &b)
Return true if the interval [a, b] includes zero, i.e., if either a or b is zero or if they have diff...
Definition: Math.h:753
typename GridT::TreeType TreeT
Definition: RayIntersector.h:90
Helper class that implements Hierarchical Digital Differential Analyzers for ray intersections agains...
Definition: DDA.h:187
typename GridT::ValueType ValueT
Definition: RayIntersector.h:517
Ray worldToIndex(const GridType &grid) const
Return a new ray in the index space of the specified grid, assuming the existing ray is represented i...
Definition: Ray.h:171
Digital Differential Analyzers specialized for VDB.
Definition: Ray.h:26
bool march(RealType &t0, RealType &t1)
Return true if the ray intersects active values, i.e. either active voxels or tiles. Only when a hit is detected are t0 and t1 updated with the corresponding entry and exit times along the INDEX ray!
Definition: RayIntersector.h:410
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h.in:116
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:644
void getWorldPos(VecT &xyz) const
Get the intersection point in world space.
Definition: RayIntersector.h:567
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h.in:202