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