HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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() ) {
103  OPENVDB_THROW(RuntimeError,
104  "LevelSetRayIntersector only supports uniform voxels!");
105  }
106  if (grid.getGridClass() != GRID_LEVEL_SET) {
107  OPENVDB_THROW(RuntimeError,
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() ) {
303  OPENVDB_THROW(RuntimeError,
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.
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() ) {
333  OPENVDB_THROW(RuntimeError,
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 
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:
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 setIndexRay(const RayT &iRay)
Return false if the ray misses the bbox of the grid.
LinearSearchImpl(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor from a grid.
void getWorldPos(VecT &xyz) const
Get the intersection point in world space.
void setTimes(RealT t0=math::Delta< RealT >::value(), RealT t1=std::numeric_limits< RealT >::max())
Definition: Ray.h:76
bool clip(const Vec3T &center, RealT radius)
Return true if this ray intersects the specified sphere.
Definition: Ray.h:217
bool intersectsWS(const RayType &wRay) const
Return true if the world-space ray intersects the level set.
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
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:362
bool setIndexRay(const RayT &iRay)
Return false if the index ray misses the bbox of the grid.
This class provides the public API for intersecting a ray with a generic (e.g. density) volume...
RealT getWorldTime() const
Return the time of intersection along the world ray.
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
GT_API const UT_StringHolder time
const math::CoordBBox & bbox() const
Return a const reference to the BBOX of the grid.
GLsizei const GLfloat * value
Definition: glcorearb.h:824
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
bool intersectsIS(const RayType &iRay) const
Return true if the index-space ray intersects the level set.
static bool test(TesterT &tester)
Definition: DDA.h:150
GLuint GLsizei GLsizei * length
Definition: glcorearb.h:795
RootNodeType & root()
Return this tree's root node.
Definition: Tree.h:281
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:239
const GridT & grid() const
Return a const reference to the input grid.
The Value Accessor Implementation and API methods. The majoirty of the API matches the API of a compa...
Definition: ValueAccessor.h:68
tree::Tree< typename RootType::template ValueConverter< bool >::Type > TreeT
bool setWorldRay(const RayT &wRay)
Return false if the ray misses the bbox of the grid.
math::Vec3< ValueType > gradient(const math::Vec3< ValueType > &xyz) const
Return the gradient in world space of the trilinear interpolation kernel.
Definition: Stencils.h:372
const AccessorType & accessor() const
Return a const reference to the ValueAccessor associated with this Stencil.
Definition: Stencils.h:207
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:1056
Vec3R getWorldPos(RealType time) const
Return the floating-point world position along the current index ray at the specified time...
VolumeRayIntersector(const GridT &grid, const math::CoordBBox &bbox)
Grid and BBox constructor.
GLdouble n
Definition: glcorearb.h:2008
bool probeValue(const Coord &xyz, ValueType &value) const
Return the active state of the value at a given coordinate as well as its value.
bool intersectsWS(const RayType &wRay, Vec3Type &world) const
Return true if the world-space ray intersects the level set.
void print(std::ostream &os=std::cout, int verboseLevel=1) const override
Print statistics, memory usage and other information about this tree.
Definition: Tree.h:1928
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
const TreeT & tree() const
Return a const reference to the (potentially dilated) bool tree used to accelerate the ray marching...
Delta for small floating-point offsets.
Definition: Math.h:155
void getIndexPos(VecT &xyz) const
Get the intersection point in index space.
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:250
bool intersectsIS(const RayType &iRay, Vec3Type &xyz, RealType &iTime) const
Return true if the index-space ray intersects the level set.
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal, RealType &wTime) const
Return true if the world-space ray intersects the level set.
bool intersectsIS(const RayType &iRay, RealType &iTime) const
Return true if the index-space ray intersects the level set.
void moveTo(const Coord &ijk)
Initialize the stencil buffer with the values of voxel (i, j, k) and its neighbors.
Definition: Stencils.h:47
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
Definition: Platform.h:146
GLdouble t
Definition: glad.h:2397
bool intersectsWS(const RayType &wRay, Vec3Type &world, Vec3Type &normal) const
Return true if the world-space ray intersects the level set.
void hits(ListType &list)
Generates a list of hits along the ray.
Implements linear iterative search for an iso-value of the level set along the direction of the ray...
void print(std::ostream &os=std::cout, int verboseLevel=1)
Print bbox, statistics, memory usage and other information.
ValueType interpolation(const math::Vec3< ValueType > &xyz) const
Return the trilinear interpolation at the normalized position.
Definition: Stencils.h:338
Helper class that implements Hierarchical Digital Differential Analyzers and is specialized for ray i...
Definition: DDA.h:144
bool intersectsWS(const RayType &wRay, Vec3Type &world, RealType &wTime) const
Return true if the world-space ray intersects the level set.
const GridType & grid() const
Return a const reference to the grid from which this stencil was constructed.
Definition: Stencils.h:203
const Vec3T & dir() const
Definition: Ray.h:99
This class provides the public API for intersecting a ray with a narrow-band level set...
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!
RealT getIndexTime() const
Return the time of intersection along the index ray.
Implementation of morphological dilation and erosion.
typename GridT::TreeType::RootNodeType RootType
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
bool intersectsWS(const RayType &wRay, RealType &wTime) const
Return true if the world-space ray intersects the level set.
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:147
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:683
void getWorldPosAndNml(VecT &xyz, VecT &nml)
Get the intersection point and normal in world space.
Vec3R getIndexPos(RealType time) const
Return the floating-point index position along the current index ray at the specified time...
Digital Differential Analyzers specialized for VDB.
VolumeRayIntersector(const VolumeRayIntersector &other)
Shallow copy constructor.
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:119
VolumeRayIntersector(const GridT &grid, int dilationCount=0)
Grid constructor.
LevelSetRayIntersector(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor.
bool setWorldRay(const RayT &wRay)
Return false if the world ray misses the bbox of the grid.
A Ray class.
Helper class that implements Hierarchical Digital Differential Analyzers for ray intersections agains...
Definition: DDA.h:187
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:74
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