OpenVDB  7.0.0
RayIntersector.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
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 
59 
60 
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");
93  static_assert(std::is_floating_point<ValueT>::value,
94  "level set grids must have scalar, floating-point value types");
95 
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 
114  const ValueT& getIsoValue() const { return mTester.getIsoValue(); }
115 
118  bool intersectsIS(const RayType& iRay) const
119  {
120  if (!mTester.setIndexRay(iRay)) return false;//missed bbox
122  }
123 
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 
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 
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 
164  bool intersectsWS(const RayType& wRay) const
165  {
166  if (!mTester.setWorldRay(wRay)) return false;//missed bbox
168  }
169 
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 
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 
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 
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 
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 
247 
248 
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 
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
311  tools::dilateVoxels(*mTree, dilationCount);
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 
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 
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 
358  ~VolumeRayIntersector() { if (mIsMaster) delete mTree; }
359 
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 
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 
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 
425  template <typename ListType>
426  inline void hits(ListType& list)
427  {
428  mHDDA.hits(mRay, mAccessor, list);
429  }
430 
433  inline Vec3R getIndexPos(RealType time) const { return mRay(time); }
434 
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 
445  const GridT& grid() const { return *mGrid; }
446 
449  const TreeT& tree() const { return *mTree; }
450 
452  const math::CoordBBox& bbox() const { return mBBox; }
453 
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 
486 
487 
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 
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 
541  const ValueT& getIsoValue() const { return mIsoValue; }
542 
546  inline bool setIndexRay(const RayT& iRay)
547  {
548  mRay = iRay;
549  return mRay.clip(mBBox);//did it hit the bbox
550  }
551 
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 
563  inline void getIndexPos(VecT& xyz) const { xyz = mRay(mTime); }
564 
567  inline void getWorldPos(VecT& xyz) const { xyz = mStencil.grid().indexToWorld(mRay(mTime)); }
568 
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 
582  inline RealT getIndexTime() const { return mTime; }
583 
585  inline RealT getWorldTime() const
586  {
587  return mTime*mStencil.grid().transform().baseMap()->applyJacobian(mRay.dir()).length();
588  }
589 
590 private:
591 
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 
603  inline const RayT& ray() const { return mRay; }
604 
606  template <typename NodeT>
607  inline bool hasNode(const Coord& ijk)
608  {
609  return mStencil.accessor().template probeConstNode<NodeT>(ijk) != nullptr;
610  }
611 
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 intersectsWS(const RayType &wRay) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:164
This class provides the public API for intersecting a ray with a narrow-band level set...
Definition: RayIntersector.h:82
Delta for small floating-point offsets.
Definition: Math.h:97
bool intersectsWS(const RayType &wRay, RealType &wTime) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:174
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:700
Definition: ValueAccessor.h:193
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
RayT::TimeSpan march()
Definition: RayIntersector.h:389
const math::CoordBBox & bbox() const
Return a const reference to the BBOX of the grid.
Definition: RayIntersector.h:452
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
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
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
This class provides the public API for intersecting a ray with a generic (e.g. density) volume...
Definition: RayIntersector.h:277
const TreeT & tree() const
Return a const reference to the (potentially dilated) bool tree used to accelerate the ray marching...
Definition: RayIntersector.h:449
typename RayT::RealType RealType
Definition: RayIntersector.h:282
typename GridT::TreeType TreeT
Definition: RayIntersector.h:90
typename GridT::ConstAccessor AccessorT
Definition: RayIntersector.h:518
typename RayT::RealType RealType
Definition: RayIntersector.h:87
typename GridT::TreeType::RootNodeType RootType
Definition: RayIntersector.h:283
void getIndexPos(VecT &xyz) const
Get the intersection point in index space.
Definition: RayIntersector.h:563
void hits(ListType &list)
Generates a list of hits along the ray.
Definition: RayIntersector.h:426
RayT RayType
Definition: RayIntersector.h:281
Helper class that implements Hierarchical Digital Differential Analyzers for ray intersections agains...
Definition: DDA.h:188
Definition: Exceptions.h:65
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
RealT getWorldTime() const
Return the time of intersection along the world ray.
Definition: RayIntersector.h:585
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:102
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
Definition: RayIntersector.h:541
VolumeRayIntersector(const VolumeRayIntersector &other)
Shallow copy constructor.
Definition: RayIntersector.h:346
typename GridT::ValueType ValueT
Definition: RayIntersector.h:89
Definition: Ray.h:26
void dilateVoxels(TreeType &tree, int iterations=1, NearestNeighbors nn=NN_FACE)
Topologically dilate all leaf-level active voxels in a tree using one of three nearest neighbor conne...
Definition: Morphology.h:826
bool setWorldRay(const RayT &wRay)
Return false if the ray misses the bbox of the grid.
Definition: RayIntersector.h:555
Vec3R getIndexPos(RealType time) const
Return the floating-point index position along the current index ray at the specified time...
Definition: RayIntersector.h:433
RayT RayType
Definition: RayIntersector.h:86
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
bool intersectsWS(const RayType &wRay, Vec3Type &world) const
Return true if the world-space ray intersects the level set.
Definition: RayIntersector.h:185
Implementation of morphological dilation and erosion.
VolumeRayIntersector(const GridT &grid, const math::CoordBBox &bbox)
Grid and BBox constructor.
Definition: RayIntersector.h:325
Definition: Exceptions.h:13
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_BEGIN
SIMD Intrinsic Headers.
Definition: Platform.h:114
LinearSearchImpl(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor from a grid.
Definition: RayIntersector.h:524
Implements linear iterative search for an iso-value of the level set along the direction of the ray...
Definition: RayIntersector.h:55
void print(std::ostream &os=std::cout, int verboseLevel=1)
Print bbox, statistics, memory usage and other information.
Definition: RayIntersector.h:458
bool clip(const Vec3T &center, RealT radius)
Return true if this ray intersects the specified sphere.
Definition: Ray.h:217
void getWorldPos(VecT &xyz) const
Get the intersection point in world space.
Definition: RayIntersector.h:567
Tag dispatch class that distinguishes topology copy constructors from deep copy constructors.
Definition: Types.h:681
Definition: Exceptions.h:63
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:379
#define OPENVDB_NO_UNREACHABLE_CODE_WARNING_END
Definition: Platform.h:115
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
const GridT & grid() const
Return a const reference to the input grid.
Definition: RayIntersector.h:445
Axis-aligned bounding box of signed integer coordinates.
Definition: Coord.h:248
GridT GridType
Definition: RayIntersector.h:280
VolumeRayIntersector(const GridT &grid, int dilationCount=0)
Grid constructor.
Definition: RayIntersector.h:296
typename GridT::ValueType ValueT
Definition: RayIntersector.h:517
bool setIndexRay(const RayT &iRay)
Return false if the ray misses the bbox of the grid.
Definition: RayIntersector.h:546
bool intersectsIS(const RayType &iRay, Vec3Type &xyz) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:139
bool setIndexRay(const RayT &iRay)
Return false if the index ray misses the bbox of the grid.
Definition: RayIntersector.h:365
Definition: Tree.h:176
Defines various finite difference stencils by means of the "curiously recurring template pattern" on ...
typename RayT::Vec3T Vec3Type
Definition: RayIntersector.h:88
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:154
bool normalize(T eps=T(1.0e-7))
this = normalized this
Definition: Vec3.h:360
Digital Differential Analyzers specialized for VDB.
RealType getWorldTime(RealType time) const
Definition: RayIntersector.h:439
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
A Ray class.
bool intersectsIS(const RayType &iRay) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:118
Vec3R getWorldPos(RealType time) const
Return the floating-point world position along the current index ray at the specified time...
Definition: RayIntersector.h:437
~VolumeRayIntersector()
Destructor.
Definition: RayIntersector.h:358
LevelSetRayIntersector(const GridT &grid, const ValueT &isoValue=zeroVal< ValueT >())
Constructor.
Definition: RayIntersector.h:99
Definition: Types.h:454
bool setWorldRay(const RayT &wRay)
Return false if the world ray misses the bbox of the grid.
Definition: RayIntersector.h:384
void getWorldPosAndNml(VecT &xyz, VecT &nml)
Get the intersection point and normal in world space.
Definition: RayIntersector.h:572
bool intersectsIS(const RayType &iRay, RealType &iTime) const
Return true if the index-space ray intersects the level set.
Definition: RayIntersector.h:128
const ValueT & getIsoValue() const
Return the iso-value used for ray-intersections.
Definition: RayIntersector.h:114
RealT getIndexTime() const
Return the time of intersection along the index ray.
Definition: RayIntersector.h:582
Helper class that implements Hierarchical Digital Differential Analyzers and is specialized for ray i...
Definition: DDA.h:145
GridT GridType
Definition: RayIntersector.h:85