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