Line data    Source code

       1             : #include "TraverseAABB.h"
       2             : #include "Intersection.h"
       3             : 
       4             : namespace elsa
       5             : {
       6             :     TraverseAABB::TraverseAABB(const BoundingBox& aabb, const RealRay_t& r) : _aabb{aabb}
       7      287989 :     {
       8             :         // compute the first intersection
       9      287989 :         calculateAABBIntersections(r);
      10      287989 :         if (!isInBoundingBox()) // early abort if necessary
      11       22970 :             return;
      12             : 
      13             :         // determine whether we go up/down or left/right
      14      265019 :         initStepDirection(r.direction());
      15             : 
      16             :         // select the initial voxel (closest to the entry point)
      17      265019 :         selectClosestVoxel();
      18             : 
      19             :         // initialize the step sizes for the step parameter
      20      265019 :         initDelta(r.direction());
      21             : 
      22             :         // initialize the maximum step parameters
      23      265019 :         initMax(r.direction());
      24      265019 :     }
      25             : 
      26             :     void TraverseAABB::updateTraverse()
      27     4401733 :     {
      28             :         // --> select the index that has lowest t value
      29     4401733 :         index_t indexToChange;
      30     4401733 :         _tMax.minCoeff(&indexToChange);
      31             : 
      32             :         // --> step into the current direction
      33     4401733 :         _currentPos(indexToChange) += _stepDirection(indexToChange);
      34             : 
      35             :         // --> update the tMax for next iteration
      36     4401733 :         _tMax(indexToChange) += _tDelta(indexToChange);
      37             : 
      38             :         // --> check if we are still in bounding box
      39     4401733 :         _isInAABB = isCurrentPositionInAABB(indexToChange);
      40     4401733 :     }
      41             : 
      42             :     real_t TraverseAABB::updateTraverseAndGetDistance()
      43     4325180 :     {
      44             :         // --> select the index that has lowest t value
      45     4325180 :         index_t indexToChange;
      46     4325180 :         _tMax.minCoeff(&indexToChange);
      47             : 
      48             :         // --> compute the distance
      49     4325180 :         real_t tEntry = _tExit;
      50     4325180 :         _tExit = _tMax(indexToChange);
      51             : 
      52             :         // --> do the update
      53     4325180 :         updateTraverse();
      54             : 
      55     4325180 :         return (_tExit - tEntry);
      56     4325180 :     }
      57             : 
      58     5332470 :     bool TraverseAABB::isInBoundingBox() const { return _isInAABB; }
      59             : 
      60             :     IndexVector_t TraverseAABB::getCurrentVoxel() const
      61     5034707 :     {
      62     5034707 :         return _currentPos.template cast<IndexVector_t::Scalar>();
      63     5034707 :     }
      64             : 
      65             :     void TraverseAABB::calculateAABBIntersections(const RealRay_t& r)
      66      288103 :     {
      67             :         // entry and exit point parameters
      68      288103 :         real_t tmin;
      69             : 
      70             :         // --> calculate intersection parameter and if the volume is hit
      71      288103 :         auto opt = Intersection::withRay(_aabb, r);
      72             : 
      73      288103 :         if (opt) { // hit!
      74      266021 :             _isInAABB = true;
      75      266021 :             tmin = opt->_tmin;
      76             : 
      77             :             // --> get points at which they intersect
      78      266021 :             _entryPoint = r.pointAt(tmin);
      79             : 
      80             :             // --> because of floating point error it can happen, that values are out of
      81             :             // the bounding box, this can lead to errors
      82      266021 :             _entryPoint =
      83      266021 :                 (_entryPoint.array() < _aabb._min.array()).select(_aabb._min, _entryPoint);
      84      266021 :             _entryPoint =
      85      266021 :                 (_entryPoint.array() > _aabb._max.array()).select(_aabb._max, _entryPoint);
      86      266021 :         }
      87      288103 :     }
      88             : 
      89             :     void TraverseAABB::initStepDirection(const RealVector_t& rd)
      90      263174 :     {
      91      263174 :         _stepDirection = rd.array().sign();
      92      263174 :     }
      93             : 
      94             :     void TraverseAABB::selectClosestVoxel()
      95      263346 :     {
      96      263346 :         RealVector_t lowerCorner = _entryPoint.array().floor();
      97      263346 :         lowerCorner =
      98      263346 :             ((lowerCorner.array() == _entryPoint.array()) && (_stepDirection.array() < 0.0))
      99      263346 :                 .select(lowerCorner.array() - 1, lowerCorner);
     100             : 
     101      263346 :         _currentPos = lowerCorner;
     102             : 
     103             :         // check if we are still inside the aabb
     104      263346 :         if ((_currentPos.array() >= _aabb._max.array()).any()
     105      265259 :             || (_currentPos.array() < _aabb._min.array()).any())
     106          12 :             _isInAABB = false;
     107      263346 :     }
     108             : 
     109             :     void TraverseAABB::initDelta(const RealVector_t& rd)
     110      266093 :     {
     111      266093 :         RealVector_t tdelta = _stepDirection.template cast<real_t>().cwiseQuotient(rd);
     112             : 
     113      266093 :         _tDelta = (Eigen::abs(rd.array()) > _EPS.array()).select(tdelta, _MAX);
     114      266093 :     }
     115             : 
     116             :     void TraverseAABB::initMax(const RealVector_t& rd)
     117      266652 :     {
     118      266652 :         RealVector_t tmax =
     119      266652 :             (((rd.array() > 0.0f).select(_currentPos.array() + 1., _currentPos)).matrix()
     120      266652 :              - _entryPoint)
     121      266652 :                 .cwiseQuotient(rd)
     122      266652 :                 .matrix();
     123             : 
     124      266652 :         _tMax = (Eigen::abs(rd.array()) > _EPS.array()).select(tmax, _MAX);
     125      266652 :     }
     126             : 
     127             :     bool TraverseAABB::isCurrentPositionInAABB(index_t index) const
     128     4409776 :     {
     129     4409776 :         return _currentPos(index) < _aabb._max(index) && _currentPos(index) >= _aabb._min(index);
     130     4409776 :     }
     131             : 
     132             : } // namespace elsa