LCOV - code coverage report
Current view: top level - projectors - TraverseAABBJosephsMethod.cpp (source / functions) Hit Total Coverage
Test: test_coverage.info.cleaned Lines: 0 127 0.0 %
Date: 2022-08-04 03:43:28 Functions: 0 14 0.0 %

          Line data    Source code
       1             : #include "TraverseAABBJosephsMethod.h"
       2             : 
       3             : namespace elsa
       4             : {
       5             : 
       6           0 :     TraverseAABBJosephsMethod::TraverseAABBJosephsMethod(const BoundingBox& aabb,
       7           0 :                                                          const RealRay_t& r)
       8           0 :         : _aabb{aabb}
       9             :     {
      10           0 :         initStepDirection(r.direction());
      11             : 
      12           0 :         RealRay_t rt(r.origin(), r.direction());
      13             : 
      14             :         // determinge length of entire intersection with AABB
      15           0 :         real_t intersectionLength = calculateAABBIntersections(rt);
      16             : 
      17           0 :         if (!isInBoundingBox())
      18           0 :             return;
      19             : 
      20           0 :         selectClosestVoxel(_entryPoint);
      21             : 
      22             :         // exit direction stored in _exitDirection
      23           0 :         (_exitPoint - (_exitPoint.array().floor() + static_cast<real_t>(0.5)).matrix())
      24           0 :             .cwiseAbs()
      25           0 :             .maxCoeff(&_exitDirection);
      26             : 
      27           0 :         moveToFirstSamplingPoint(r.direction(), intersectionLength);
      28             : 
      29             :         // last pixel/voxel handled separately, reduce box dimensionality for easy detection
      30             :         int mainDir;
      31           0 :         rt.direction().cwiseAbs().maxCoeff(&mainDir);
      32           0 :         real_t prevBorder = std::floor(_exitPoint(mainDir));
      33           0 :         if (prevBorder == _exitPoint(mainDir) && rt.direction()[mainDir] > 0) {
      34           0 :             prevBorder--;
      35             :         }
      36           0 :         if (rt.direction()[mainDir] < 0) {
      37           0 :             _aabb._min[mainDir] = prevBorder + 1;
      38             :         } else {
      39           0 :             _aabb._max[mainDir] = prevBorder;
      40             :         }
      41           0 :     }
      42             : 
      43           0 :     void TraverseAABBJosephsMethod::updateTraverse()
      44             :     {
      45           0 :         _fractionals += _nextStep;
      46           0 :         for (int i = 0; i < _aabb._dim; i++) {
      47             :             //*0.5 because a change of 1 spacing corresponds to a change of
      48             :             // fractionals from -0.5 to 0.5  0.5 or -0.5 = voxel border is still
      49             :             // ok
      50           0 :             if (std::abs(_fractionals(i)) > 0.5) {
      51           0 :                 _fractionals(i) -= _stepDirection(i);
      52           0 :                 _currentPos(i) += _stepDirection(i);
      53             :                 // --> is the traverse algorithm still in the bounding box?
      54           0 :                 _isInAABB = _isInAABB && isCurrentPositionInAABB(i);
      55             :             }
      56             :         }
      57           0 :         switch (_stage) {
      58           0 :             case FIRST:
      59           0 :                 _nextStep.cwiseAbs().maxCoeff(&_ignoreDirection);
      60           0 :                 if (_isInAABB) {
      61             :                     // now ignore main direction
      62           0 :                     _nextStep /= std::abs(_nextStep(_ignoreDirection));
      63           0 :                     _fractionals(_ignoreDirection) = 0;
      64           0 :                     _intersectionLength = _nextStep.norm();
      65           0 :                     _stage = INTERIOR;
      66             :                 }
      67             :                 [[fallthrough]];
      68             :             case INTERIOR:
      69           0 :                 if (!_isInAABB) {
      70             :                     // revert to exit position and adjust values
      71           0 :                     real_t prevBorder = _nextStep(_ignoreDirection) < 0
      72           0 :                                             ? _aabb._min[_ignoreDirection]
      73           0 :                                             : _aabb._max[_ignoreDirection];
      74           0 :                     _nextStep.normalize();
      75           0 :                     _intersectionLength =
      76           0 :                         (_exitPoint[_ignoreDirection] - prevBorder) / (_nextStep[_ignoreDirection]);
      77             :                     // determine exit direction (the exit coordinate furthest from the center of the
      78             :                     // volume)
      79           0 :                     _ignoreDirection = _exitDirection;
      80             :                     // move to last sampling point
      81           0 :                     RealVector_t currentPosition = _exitPoint - _intersectionLength * _nextStep / 2;
      82           0 :                     selectClosestVoxel(currentPosition);
      83           0 :                     initFractionals(currentPosition);
      84           0 :                     _isInAABB = true;
      85           0 :                     _stage = LAST;
      86           0 :                 }
      87           0 :                 break;
      88           0 :             case LAST:
      89           0 :                 _isInAABB = false;
      90           0 :                 break;
      91           0 :             default:
      92           0 :                 break;
      93             :         }
      94           0 :     }
      95             : 
      96           0 :     const RealVector_t& TraverseAABBJosephsMethod::getFractionals() const { return _fractionals; }
      97           0 :     int TraverseAABBJosephsMethod::getIgnoreDirection() const { return _ignoreDirection; }
      98             : 
      99           0 :     void TraverseAABBJosephsMethod::initFractionals(const RealVector_t& currentPosition)
     100             :     {
     101           0 :         for (int i = 0; i < _aabb._dim; i++)
     102           0 :             _fractionals(i) =
     103           0 :                 static_cast<real_t>(std::abs(currentPosition(i)) - _currentPos(i) - 0.5);
     104           0 :     }
     105             : 
     106           0 :     void TraverseAABBJosephsMethod::moveToFirstSamplingPoint(const RealVector_t& rd,
     107             :                                                              real_t intersectionLength)
     108             :     {
     109             :         // determine main direction
     110           0 :         rd.cwiseAbs().maxCoeff(&_ignoreDirection);
     111             : 
     112             :         // initialize _nextStep as the step for interior pixels
     113           0 :         _nextStep(_ignoreDirection) = _stepDirection(_ignoreDirection);
     114           0 :         for (int i = 0; i < _aabb._dim; ++i) {
     115           0 :             if (i != _ignoreDirection) {
     116             :                 // tDelta(i) is given in relation to tDelta(_ignoreDirection)
     117           0 :                 _nextStep(i) = _nextStep(_ignoreDirection) * rd(i) / rd(_ignoreDirection);
     118             :             }
     119             :         }
     120             : 
     121           0 :         RealVector_t currentPosition = _entryPoint;
     122             :         // move to first sampling point
     123           0 :         if (std::abs(_entryPoint[_ignoreDirection] - _exitPoint[_ignoreDirection])
     124           0 :             <= static_cast<real_t>(1.)) {
     125           0 :             _stage = LAST;
     126             :             // use midpoint of intersection as the interpolation value
     127           0 :             currentPosition += rd * intersectionLength / 2;
     128           0 :             _intersectionLength = intersectionLength;
     129             :         } else {
     130             :             // find distance to next plane orthogonal to main direction
     131           0 :             real_t nextBoundary = std::trunc(_currentPos(_ignoreDirection));
     132           0 :             if (_stepDirection(_ignoreDirection) > 0)
     133           0 :                 nextBoundary += static_cast<real_t>(1.);
     134             : 
     135             :             real_t distToBoundary =
     136           0 :                 (nextBoundary - _entryPoint[_ignoreDirection]) / rd[_ignoreDirection];
     137             : 
     138           0 :             currentPosition += rd * distToBoundary / 2;
     139           0 :             _nextStep = _nextStep / 2 + _nextStep * (distToBoundary / (2 * _nextStep.norm()));
     140           0 :             _intersectionLength = distToBoundary;
     141             : 
     142             :             // determine entry direction (the entry coordinate furthest from the center of the
     143             :             // volume)
     144           0 :             (_entryPoint - (_entryPoint.array().floor() + static_cast<real_t>(0.5)).matrix())
     145           0 :                 .cwiseAbs()
     146           0 :                 .maxCoeff(&_ignoreDirection);
     147             :         }
     148           0 :         selectClosestVoxel(currentPosition);
     149             : 
     150             :         // init fractionals
     151           0 :         initFractionals(currentPosition);
     152           0 :     }
     153             : 
     154           0 :     inline bool TraverseAABBJosephsMethod::isCurrentPositionInAABB(index_t index) const
     155             :     {
     156           0 :         return _currentPos(index) < _aabb._max(index) && _currentPos(index) >= _aabb._min(index);
     157             :     }
     158             : 
     159           0 :     void TraverseAABBJosephsMethod::selectClosestVoxel(const RealVector_t& currentPosition)
     160             :     {
     161           0 :         RealVector_t lowerCorner = currentPosition.array().floor();
     162             :         lowerCorner =
     163           0 :             (((lowerCorner.array() == currentPosition.array()) && (_stepDirection.array() < 0.0))
     164           0 :              || currentPosition.array() == _aabb._max.array())
     165           0 :                 .select(lowerCorner.array() - 1, lowerCorner);
     166             : 
     167             :         // --> If ray is parallel and we are close, choose the next previous/next voxel
     168           0 :         _currentPos = lowerCorner;
     169             : 
     170             :         // check if we are still inside the aabb
     171           0 :         if ((_currentPos.array() >= _aabb._max.array()).any()
     172           0 :             || (_currentPos.array() < _aabb._min.array()).any())
     173           0 :             _isInAABB = false;
     174           0 :     }
     175             : 
     176           0 :     RealVector_t floorFallback(const RealVector_t& v)
     177             :     {
     178           0 :         RealVector_t ret(v.rows());
     179           0 :         for (int i = 0; i < v.rows(); i++)
     180           0 :             ret(i) = std::floor(v(i));
     181           0 :         return ret;
     182           0 :     }
     183             : 
     184           0 :     void TraverseAABBJosephsMethod::initStepDirection(const RealVector_t& rd)
     185             :     {
     186           0 :         _stepDirection = rd.array().sign();
     187           0 :     }
     188             : 
     189           0 :     real_t TraverseAABBJosephsMethod::calculateAABBIntersections(const RealRay_t& ray)
     190             :     {
     191             :         real_t tmin;
     192             : 
     193             :         // --> calculate intersection parameter and if the volume is hit
     194           0 :         auto opt = Intersection::withRay(_aabb, ray);
     195             : 
     196           0 :         if (opt) {
     197           0 :             _isInAABB = true;
     198           0 :             tmin = opt->_tmin;
     199             : 
     200             :             // --> get points at which they intersect
     201           0 :             _entryPoint = ray.pointAt(tmin);
     202             : 
     203             :             // --> because of floating point error it can happen, that values are out of
     204             :             // the bounding box, this can lead to errors
     205             :             _entryPoint =
     206           0 :                 (_entryPoint.array() < _aabb._min.array()).select(_aabb._min, _entryPoint);
     207             :             _entryPoint =
     208           0 :                 (_entryPoint.array() > _aabb._max.array()).select(_aabb._max, _entryPoint);
     209             : 
     210             :             // exit point stored in _exitPoint
     211           0 :             _exitPoint = ray.pointAt(opt->_tmax);
     212           0 :             _exitPoint = (_exitPoint.array() < _aabb._min.array()).select(_aabb._min, _exitPoint);
     213           0 :             _exitPoint = (_exitPoint.array() > _aabb._max.array()).select(_aabb._max, _exitPoint);
     214             : 
     215           0 :             return opt->_tmax - opt->_tmin;
     216             :         } else {
     217           0 :             return 0;
     218             :         }
     219             :     }
     220             : 
     221           0 :     bool TraverseAABBJosephsMethod::isInBoundingBox() const { return _isInAABB; }
     222             : 
     223           0 :     real_t TraverseAABBJosephsMethod::getIntersectionLength() const { return _intersectionLength; }
     224             : 
     225           0 :     IndexVector_t TraverseAABBJosephsMethod::getCurrentVoxel() const
     226             :     {
     227           0 :         return _currentPos.template cast<index_t>();
     228             :     }
     229             : 
     230             : } // namespace elsa

Generated by: LCOV version 1.14