Line data    Source code

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