Merge branch 'feature/cpu-projectors' into 'master'

Feature/cpu projectors

See merge request !2

elsa/core/DataContainer.h

@@ -6,6 +6,7 @@
 #include "DataHandler.h"
 
 #include <memory>
+#include <type_traits>
 
 namespace elsa
 {
@@ -107,6 +108,21 @@ namespace elsa
         /// return an element by n-dimensional coordinate as read-only (not bounds-checked!)
         const data_t& operator()(IndexVector_t coordinate) const;
 
+        template <typename idx0_t, typename... idx_t, 
+                    typename = std::enable_if_t<std::is_integral_v<idx0_t> && (... && std::is_integral_v<idx_t>)>>
+        data_t& operator()(idx0_t idx0, idx_t... indices) {
+            IndexVector_t coordinate(sizeof...(indices)+1);
+            ((coordinate<<idx0) , ... , indices);
+            return operator()(coordinate);
+        }
+
+        template <typename idx0_t, typename... idx_t, 
+                    typename = std::enable_if_t<std::is_integral_v<idx0_t> && (... && std::is_integral_v<idx_t>)>>
+        const data_t& operator()(idx0_t idx0, idx_t... indices) const{
+            IndexVector_t coordinate(sizeof...(indices)+1);
+            ((coordinate<<idx0) , ... , indices);
+            return operator()(coordinate);
+        }
 
         /// return the dot product of this signal with the one from container other
         data_t dot(const DataContainer<data_t>& other) const;

elsa/core/tests/test_DataContainer.cpp

@@ -140,14 +140,17 @@ SCENARIO("Element-wise access of DataContainers") {
             THEN("it works as expected when using indices/coordinates") {
                 REQUIRE(dc[index] == 0.0_a);
                 REQUIRE(dc(coord) == 0.0_a);
+                REQUIRE(dc(17, 4) == 0.0_a);
 
                 dc[index] = 2.2;
                 REQUIRE(dc[index] == 2.2_a);
                 REQUIRE(dc(coord) == 2.2_a);
+                REQUIRE(dc(17, 4) == 2.2_a);
 
                 dc(coord) = 3.3;
                 REQUIRE(dc[index] == 3.3_a);
                 REQUIRE(dc(coord) == 3.3_a);
+                REQUIRE(dc(17, 4) == 3.3_a);
             }
         }
     }

elsa/projectors/CMakeLists.txt

@@ -12,7 +12,10 @@ set(MODULE_HEADERS
         BoundingBox.h
         Intersection.h
         TraverseAABB.h
-        BinaryMethod.h)
+        TraverseAABBJosephsMethod.h
+        BinaryMethod.h
+        SiddonsMethod.h
+        JosephsMethod.h)
 
 # list all the code files of the module
 set(MODULE_SOURCES
@@ -20,7 +23,10 @@ set(MODULE_SOURCES
         BoundingBox.cpp
         Intersection.cpp
         TraverseAABB.cpp
-        BinaryMethod.cpp)
+        TraverseAABBJosephsMethod.cpp
+        BinaryMethod.cpp
+        SiddonsMethod.cpp
+        JosephsMethod.cpp)
 
 
 # build the module library

elsa/projectors/JosephsMethod.cpp

+#include "JosephsMethod.h"
+#include "Timer.h"
+#include "TraverseAABBJosephsMethod.h"
+
+#include <stdexcept>
+
+namespace elsa
+{
+    template <typename data_t>
+    JosephsMethod<data_t>::JosephsMethod(const DataDescriptor& domainDescriptor, const DataDescriptor& rangeDescriptor,
+                                         const std::vector<Geometry>& geometryList, Interpolation interpolation)
+        : LinearOperator<data_t>(domainDescriptor, rangeDescriptor), _geometryList{geometryList},
+          _boundingBox{domainDescriptor.getNumberOfCoefficientsPerDimension()}, _interpolation{interpolation}
+    {
+        auto dim = _domainDescriptor->getNumberOfDimensions();
+        if (dim!=2 && dim != 3) {
+            throw std::invalid_argument("JosephsMethod:only supporting 2d/3d operations");
+        }
+
+        if (dim != _rangeDescriptor->getNumberOfDimensions()) {
+            throw std::invalid_argument("JosephsMethod: domain and range dimension need to match");
+        }
+
+        if (_geometryList.empty()) {
+            throw std::invalid_argument("JosephsMethod: geometry list was empty");
+        }
+    }
+
+    template <typename data_t>
+    void JosephsMethod<data_t>::_apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax)
+    {
+        Timer<> timeguard("JosephsMethod", "apply");
+        traverseVolume<false>(x, Ax);
+    }
+
+    template <typename data_t>
+    void JosephsMethod<data_t>::_applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty)
+    {
+        Timer<> timeguard("JosephsMethod", "applyAdjoint");
+        traverseVolume<true>(y, Aty);
+    }
+
+    template <typename data_t>
+    JosephsMethod<data_t>* JosephsMethod<data_t>::cloneImpl() const{
+        return new JosephsMethod(*_domainDescriptor, *_rangeDescriptor, _geometryList);
+    }
+
+    template <typename data_t>
+    bool JosephsMethod<data_t>::isEqual(const LinearOperator<data_t>& other) const
+    {
+        if (!LinearOperator<data_t>::isEqual(other))
+            return false;
+
+        auto otherJM = dynamic_cast<const JosephsMethod*>(&other);
+        if (!otherJM)
+            return false;
+
+        if (_geometryList != otherJM->_geometryList || _interpolation != otherJM->_interpolation)
+            return false;
+
+        return true;
+    }
+
+    template <typename data_t>
+    template <bool adjoint>
+    void JosephsMethod<data_t>::traverseVolume(const DataContainer<data_t>& vector, DataContainer<data_t>& result) const
+    {
+        if(adjoint) result = 0;
+
+        const index_t sizeOfRange = _rangeDescriptor->getNumberOfCoefficients();
+        const auto rangeDim = _rangeDescriptor->getNumberOfDimensions();
+
+        // iterate over all rays
+    #pragma omp parallel for
+        for (index_t ir = 0; ir < sizeOfRange; ir++) {
+            Ray ray = computeRayToDetector(ir,rangeDim);
+
+            // --> setup traversal algorithm
+            TraverseAABBJosephsMethod traverse(_boundingBox, ray);
+            
+            if(!adjoint) result[ir] = 0;
+
+            // Make steps through the volume
+            while (traverse.isInBoundingBox()) {
+            
+                IndexVector_t currentVoxel = traverse.getCurrentVoxel();
+                float intersection = traverse.getIntersectionLength();
+
+                // to avoid code duplicates for apply and applyAdjoint
+                index_t from;
+                index_t to;
+                if (adjoint) {
+                    to = _domainDescriptor->getIndexFromCoordinate(currentVoxel);
+                    from = ir;
+                } else {
+                    to = ir;
+                    from = _domainDescriptor->getIndexFromCoordinate(currentVoxel);
+                }
+                
+                switch (_interpolation) {
+                case Interpolation::LINEAR:
+                    LINEAR(vector, result, traverse.getFractionals(), adjoint, rangeDim,
+                        currentVoxel, intersection, from, to, traverse.getIgnoreDirection());
+                    break;
+                case Interpolation::NN:
+                    if (adjoint) {
+                        #pragma omp atomic
+                        result[to] += intersection * vector[from];
+                    }
+                    else {
+                        result[to] += intersection * vector[from];
+                    }
+                    break;
+                }
+
+                // update Traverse
+                traverse.updateTraverse();
+            }
+        }
+    }
+
+    template <typename data_t>
+    typename JosephsMethod<data_t>::Ray JosephsMethod<data_t>::computeRayToDetector(index_t detectorIndex, index_t dimension) const
+    {
+        auto detectorCoord = _rangeDescriptor->getCoordinateFromIndex(detectorIndex);
+
+        //center of detector pixel is 0.5 units away from the corresponding detector coordinates
+        auto geometry = _geometryList.at(detectorCoord(dimension - 1));
+        auto [ro, rd] = geometry.computeRayTo(detectorCoord.block(0, 0, dimension-1, 1).template cast<real_t>().array() + 0.5);
+
+        return Ray(ro, rd);
+    }
+
+    template <typename data_t>
+    void JosephsMethod<data_t>::LINEAR(const DataContainer<data_t>& vector, DataContainer<data_t>& result,
+        const RealVector_t& fractionals, bool adjoint, int domainDim,
+        const IndexVector_t& currentVoxel, float intersection, index_t from, index_t to,
+        int mainDirection) const
+    {
+        float weight = intersection;
+        IndexVector_t interpol = currentVoxel;
+
+        //handle diagonal if 3D
+        if(domainDim==3) {
+            for (int i=0;i<domainDim;i++) {
+                if(i != mainDirection) {
+                    weight *= fabs(fractionals(i));
+                    interpol(i) += (fractionals(i) < 0.0) ? -1 : 1;
+                    // mirror values at border if outside the volume
+                    if(interpol(i)<_boundingBox._min(i) || interpol(i)>_boundingBox._max(i))
+                        interpol(i) = _boundingBox._min(i);
+                    else if(interpol(i) == _boundingBox._max(i))
+                        interpol(i) = _boundingBox._max(i)-1;
+                }
+            }
+            if (adjoint) {
+                #pragma omp atomic
+                result(interpol) += weight * vector[from];
+            } else {
+                result[to] += weight * vector(interpol);
+            }
+        }
+
+        // handle current voxel
+        weight = intersection*(1-fractionals.array().abs()).prod()/(1-fabs(fractionals(mainDirection)));
+        if (adjoint) {
+            #pragma omp atomic
+            result[to] += weight * vector[from];
+        }
+        else {
+            result[to] += weight * vector[from];
+        }
+
+        // handle neighbors not along the main direction
+        for (int i = 0; i < domainDim; i++) {
+            if (i != mainDirection) {
+                float weightn = weight * fabs(fractionals(i))/(1-fabs(fractionals(i)));
+                interpol = currentVoxel;
+                interpol(i) += (fractionals(i) < 0.0) ? -1 : 1;
+                
+                // mirror values at border if outside the volume
+                if(interpol(i)<_boundingBox._min(i) || interpol(i)>_boundingBox._max(i))
+                        interpol(i) = _boundingBox._min(i);
+                else if(interpol(i) == _boundingBox._max(i))
+                    interpol(i) = _boundingBox._max(i)-1;
+                
+                if (adjoint) {
+                    #pragma omp atomic
+                    result(interpol) += weightn * vector[from];
+                } else {
+                    result[to] += weightn * vector(interpol);
+                }
+            }
+        }
+    }
+
+    // ------------------------------------------
+    // explicit template instantiation
+    template class JosephsMethod<float>;
+    template class JosephsMethod<double>;
+
+} // namespace elsa

elsa/projectors/JosephsMethod.h

+#pragma once
+
+#include "LinearOperator.h"
+#include "Geometry.h"
+#include "BoundingBox.h"
+
+#include <vector>
+#include <utility>
+
+#include <Eigen/Geometry>
+
+namespace elsa
+{
+    /**
+     * \brief Operator representing the discretized X-ray transform in 2d/3d using Joseph's method.
+     *
+     * \author Christoph Hahn - initial implementation
+     * \author Maximilian Hornung - modularization
+     * \author Nikola Dinev - fixes
+     * 
+     * \tparam data_t data type for the domain and range of the operator, defaulting to real_t
+     * 
+     * The volume is traversed along the rays as specified by the Geometry. For interior voxels
+     * the sampling point is located in the middle of the two planes orthogonal to the main
+     * direction of the ray. For boundary voxels the sampling point is located at the center of the
+     * ray intersection with the voxel.
+     * 
+     * The geometry is represented as a list of projection matrices (see class Geometry), one for each
+     * acquisition pose.
+     * 
+     * Two modes of interpolation are available:
+     * NN (NearestNeighbours) takes the value of the pixel/voxel containing the point
+     * LINEAR performs linear interpolation for the nearest 2 pixels (in 2D)
+     * or the nearest 4 voxels (in 3D). 
+     * 
+     * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
+     * This projector is matched.
+     */
+    template <typename data_t = real_t>
+    class JosephsMethod : public LinearOperator<data_t> {
+    public:
+        /// Available interpolation modes
+        enum class Interpolation { NN, LINEAR };
+
+        /**
+         * \brief Constructor for Joseph's traversal method.
+         *
+         * \param[in] domainDescriptor describing the domain of the operator (the volume)
+         * \param[in] rangeDescriptor describing the range of the operator (the sinogram)
+         * \param[in] geometryList vector containing the geometries for the acquisition poses
+         *
+         * The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]),
+         * the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
+         */
+        JosephsMethod(const DataDescriptor& domainDescriptor, const DataDescriptor& rangeDescriptor,
+                const std::vector<Geometry>& geometryList, Interpolation interpolation = Interpolation::LINEAR);
+
+        /// default destructor
+        ~JosephsMethod() = default;
+
+    protected:
+        /// apply the binary method (i.e. forward projection)
+        void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) override;
+
+        /// apply the adjoint of the  binary method (i.e. backward projection)
+        void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) override;
+
+        /// implement the polymorphic clone operation
+        JosephsMethod<data_t>* cloneImpl() const override;
+
+        /// implement the polymorphic comparison operation
+        bool isEqual(const LinearOperator<data_t>& other) const override;
+
+    private:
+        /// the bounding box of the volume
+        BoundingBox _boundingBox;
+
+        /// the geometry list
+        std::vector<Geometry> _geometryList;
+
+        /// the interpolation mode
+        Interpolation _interpolation;
+
+        /// the traversal routine (for both apply/applyAdjoint)
+        template <bool adjoint>
+        void traverseVolume(const DataContainer<data_t>& vector, DataContainer<data_t>& result) const;
+
+        /// convenience typedef for ray
+        using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
+
+        /**
+         * \brief computes the ray to the middle of the detector element
+         *
+         * \param[in] detectorIndex the index of the detector element
+         * \param[in] dimension the dimension of the detector (1 or 2)
+         *
+         * \returns the ray
+         */
+        Ray computeRayToDetector(index_t detectorIndex, index_t dimension) const;
+
+
+        /**
+         * \brief  Linear interpolation, works in any dimension
+         *
+         * \param vector the input DataContainer
+         * \param result DataContainer for results
+         * \param fractionals the fractional numbers used in the interpolation
+         * \param adjoint true for backward projection, false for forward
+         * \param domainDim number of dimensions
+         * \param currentVoxel coordinates of voxel for interpolation
+         * \param intersection weighting for the interpolated values depending on the incidence angle
+         * \param from index of the current vector position
+         * \param to index of the current result position
+         */
+        void LINEAR(const DataContainer<data_t>& vector, DataContainer<data_t>& result, const RealVector_t& fractionals, bool adjoint,
+                int domainDim, const IndexVector_t& currentVoxel, float intersection, index_t from, index_t to, int mainDirection) const;    
+
+        /// lift from base class
+        using LinearOperator<data_t>::_domainDescriptor;
+        using LinearOperator<data_t>::_rangeDescriptor;
+    };
+} // namespace elsa

elsa/projectors/SiddonsMethod.cpp

+#include "SiddonsMethod.h"
+#include "Timer.h"
+#include "TraverseAABB.h"
+
+#include <stdexcept>
+
+namespace elsa
+{
+    template <typename data_t>
+    SiddonsMethod<data_t>::SiddonsMethod(const DataDescriptor& domainDescriptor, const DataDescriptor& rangeDescriptor,
+                               const std::vector<Geometry>& geometryList)
+            : LinearOperator<data_t>(domainDescriptor, rangeDescriptor), _geometryList{geometryList},
+              _boundingBox(domainDescriptor.getNumberOfCoefficientsPerDimension())
+    {
+        auto dim = _domainDescriptor->getNumberOfDimensions();
+        if (dim != _rangeDescriptor->getNumberOfDimensions()) {
+            throw std::logic_error("SiddonsMethod: domain and range dimension need to match");
+        }
+
+        if (dim!=2 && dim != 3) {
+            throw std::logic_error("SiddonsMethod: only supporting 2d/3d operations");
+        }
+
+        if (_geometryList.empty()) {
+            throw std::logic_error("SiddonsMethod: geometry list was empty");
+        }
+    }
+
+    template <typename data_t>
+    void SiddonsMethod<data_t>::_apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax)
+    {
+        Timer t("SiddonsMethod", "apply");
+        traverseVolume<false>(x, Ax);
+    }
+
+    template <typename data_t>
+    void SiddonsMethod<data_t>::_applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty)
+    {
+        Timer t("SiddonsMethod", "applyAdjoint");
+        traverseVolume<true>(y, Aty);
+    }
+
+    template <typename data_t>
+    SiddonsMethod<data_t>* SiddonsMethod<data_t>::cloneImpl() const
+    {
+        return new SiddonsMethod(*_domainDescriptor, *_rangeDescriptor, _geometryList);
+    }
+
+    template <typename data_t>
+    bool SiddonsMethod<data_t>::isEqual(const LinearOperator<data_t>& other) const
+    {
+        if (!LinearOperator<data_t>::isEqual(other))
+            return false;
+
+        auto otherSM = dynamic_cast<const SiddonsMethod*>(&other);
+        if (!otherSM)
+            return false;
+
+        if (_geometryList != otherSM->_geometryList)
+            return false;
+
+        return true;
+    }
+
+    template<typename data_t>
+    template<bool adjoint>
+    void SiddonsMethod<data_t>::traverseVolume(const DataContainer<data_t>& vector, DataContainer<data_t>& result) const
+    {
+        index_t maxIterations{0};
+        if (adjoint) {
+            maxIterations = vector.getSize();
+            result = 0; // initialize volume to 0, because we are not going to hit every voxel!
+        } else
+            maxIterations = result.getSize();
+
+        const auto rangeDim = _rangeDescriptor->getNumberOfDimensions();
+
+        // --> loop either over every voxel that should  updated or every detector
+        // cell that should be calculated
+#pragma omp parallel for
+        for (size_t rangeIndex = 0; rangeIndex < maxIterations; ++rangeIndex)
+        {
+            // --> get the current ray to the detector center
+            auto ray = computeRayToDetector(rangeIndex, rangeDim);
+
+            // --> setup traversal algorithm
+            TraverseAABB traverse(_boundingBox, ray);
+
+            if(!adjoint) 
+                result[rangeIndex] = 0;
+
+            // --> initial index to access the data vector
+            auto dataIndexForCurrentVoxel = _domainDescriptor->getIndexFromCoordinate(traverse.getCurrentVoxel());
+
+            while (traverse.isInBoundingBox())
+            {
+                
+                auto weight = traverse.updateTraverseAndGetDistance();
+                // --> update result depending on the operation performed
+                if (adjoint)
+#pragma omp atomic
+                    result[dataIndexForCurrentVoxel] += vector[rangeIndex] * weight;
+                else
+                    result[rangeIndex] += vector[dataIndexForCurrentVoxel] * weight;
+                
+                dataIndexForCurrentVoxel = _domainDescriptor->getIndexFromCoordinate(traverse.getCurrentVoxel());
+            }
+        }
+    }
+
+    template <typename data_t>
+    typename SiddonsMethod<data_t>::Ray SiddonsMethod<data_t>::computeRayToDetector(index_t detectorIndex, index_t dimension) const
+    {
+        auto detectorCoord = _rangeDescriptor->getCoordinateFromIndex(detectorIndex);
+
+        //center of detector pixel is 0.5 units away from the corresponding detector coordinates
+        auto geometry = _geometryList.at(detectorCoord(dimension - 1));
+        auto [ro, rd] = geometry.computeRayTo(detectorCoord.block(0, 0, dimension-1, 1).template cast<real_t>().array() + 0.5);
+
+        return Ray(ro, rd);
+    }
+
+    // ------------------------------------------
+    // explicit template instantiation
+    template class SiddonsMethod<float>;
+    template class SiddonsMethod<double>;
+}

elsa/projectors/SiddonsMethod.h

+#pragma once
+
+#include "LinearOperator.h"
+#include "Geometry.h"
+#include "BoundingBox.h"
+
+#include <vector>
+#include <utility>
+
+#include <Eigen/Geometry>
+
+namespace elsa
+{
+    /**
+     * \brief Operator representing the discretized X-ray transform in 2d/3d using Siddon's method.
+     *
+     * \author David Frank - initial code
+     * \author Nikola Dinev - modularization, fixes
+     *
+     * \tparam data_t data type for the domain and range of the operator, defaulting to real_t
+     *
+     * The volume is traversed along the rays as specified by the Geometry. Each ray is traversed in a
+     * continguous fashion (i.e. along long voxel borders, not diagonally) and each traversed voxel is
+     * counted as a hit with weight according to the length of the path of the ray through the voxel.
+     *
+     * The geometry is represented as a list of projection matrices (see class Geometry), one for each
+     * acquisition pose.
+     *
+     * Forward projection is accomplished using apply(), backward projection using applyAdjoint().
+     * This projector is matched.
+     */
+     template <typename data_t = real_t>
+     class SiddonsMethod : public LinearOperator<data_t> {
+     public:
+         /**
+          * \brief Constructor for Siddon's method traversal.
+          *
+          * \param[in] domainDescriptor describing the domain of the operator (the volume)
+          * \param[in] rangeDescriptor describing the range of the operator (the sinogram)
+          * \param[in] geometryList vector containing the geometries for the acquisition poses
+          *
+          * The domain is expected to be 2 or 3 dimensional (volSizeX, volSizeY, [volSizeZ]),
+          * the range is expected to be matching the domain (detSizeX, [detSizeY], acqPoses).
+          */
+         SiddonsMethod(const DataDescriptor& domainDescriptor, const DataDescriptor& rangeDescriptor,
+                      const std::vector<Geometry>& geometryList);
+
+         /// default destructor
+         ~SiddonsMethod() override = default;
+
+     protected:
+         /// apply Siddon's method (i.e. forward projection)
+         void _apply(const DataContainer<data_t>& x, DataContainer<data_t>& Ax) override;
+
+         /// apply the adjoint of Siddon's method (i.e. backward projection)
+         void _applyAdjoint(const DataContainer<data_t>& y, DataContainer<data_t>& Aty) override;
+
+         /// implement the polymorphic clone operation
+         SiddonsMethod<data_t>* cloneImpl() const override;
+
+         /// implement the polymorphic comparison operation
+         bool isEqual(const LinearOperator<data_t>& other) const override;
+
+     private:
+         /// the bounding box of the volume
+         BoundingBox _boundingBox;
+
+         /// the geometry list
+         std::vector<Geometry> _geometryList;
+
+         /// the traversal routine (for both apply/applyAdjoint)
+         template <bool adjoint>
+         void traverseVolume(const DataContainer<data_t>& vector, DataContainer<data_t>& result) const;
+
+         /// convenience typedef for ray
+         using Ray = Eigen::ParametrizedLine<real_t, Eigen::Dynamic>;
+
+         /**
+          * \brief computes the ray to the middle of the detector element
+          *
+          * \param[in] detectorIndex the index of the detector element
+          * \param[in] dimension the dimension of the detector (1 or 2)
+          *
+          * \returns the ray
+          */
+         Ray computeRayToDetector(index_t detectorIndex, index_t dimension) const;
+
+         /// lift from base class
+         using LinearOperator<data_t>::_domainDescriptor;
+         using LinearOperator<data_t>::_rangeDescriptor;
+     };
+
+} // namespace elsa

elsa/projectors/TraverseAABBJosephsMethod.cpp

+#include "TraverseAABBJosephsMethod.h"
+
+namespace elsa
+{
+
+TraverseAABBJosephsMethod::TraverseAABBJosephsMethod(const BoundingBox& aabb, const Ray& r)
+    : _aabb{aabb}
+{
+    initStepDirection(r.direction());
+    
+    Ray rt(r.origin(),r.direction());
+
+    // determinge length of entire intersection with AABB
+    real_t intersectionLength = calculateAABBIntersections(rt);
+    
+    if(!isInBoundingBox()) return;
+
+    selectClosestVoxel(rt.direction(),_entryPoint);
+
+    // exit direction stored in _exitDirection
+    (_exitPoint - (_aabb._max - _aabb._min)/2).cwiseAbs().maxCoeff(&_exitDirection);
+    
+    moveToFirstSamplingPoint(r.direction(),intersectionLength);
+
+    // last pixel/voxel handled separately, reduce box dimensionality for easy detection
+    int mainDir;
+    rt.direction().cwiseAbs().maxCoeff(&mainDir);
+    real_t prevBorder = std::floor(_exitPoint(mainDir));
+    if (prevBorder==_exitPoint(mainDir) && rt.direction()[mainDir]>0) {
+        prevBorder--;
+    }
+    if(rt.direction()[mainDir]<0) {
+        _aabb._min[mainDir] = prevBorder+1;
+    }
+    else {
+        _aabb._max[mainDir] = prevBorder;
+    }
+}
+
+void TraverseAABBJosephsMethod::updateTraverse() {
+    _fractionals += _nextStep;
+    for (int i = 0; i < _aabb._dim; i++) {
+        //*0.5 because a change of 1 spacing corresponds to a change of
+        // fractionals from -0.5 to 0.5  0.5 or -0.5 = voxel border is still
+        // ok
+        if (fabs(_fractionals(i)) > 0.5) {
+            _fractionals(i) -= _stepDirection(i);
+            _currentPos(i) += _stepDirection(i);
+            // --> is the traverse algorithm still in the bounding box?
+            _isInAABB = _isInAABB && isCurrentPositionInAABB(i);
+        }
+    }
+    switch (_stage)
+    {
+    case FIRST:
+        if (_isInAABB) {
+            // now ignore main direction
+            _nextStep.cwiseAbs().maxCoeff(&_ignoreDirection);
+            _nextStep /= fabs(_nextStep(_ignoreDirection));
+            _fractionals(_ignoreDirection) = 0;
+            _intersectionLength = _nextStep.norm();
+            _stage = INTERIOR;
+        }
+    case INTERIOR:
+        if(!_isInAABB) {
+            // revert to exit position and adjust values
+            real_t prevBorder = _nextStep(_ignoreDirection)<0 ? _aabb._min[_ignoreDirection] : _aabb._max[_ignoreDirection];
+            _nextStep.normalize();
+            _intersectionLength = (_exitPoint[_ignoreDirection] - prevBorder)/(_nextStep[_ignoreDirection]);
+            // determine exit direction (the exit coordinate furthest from the center of the volume)
+            _ignoreDirection = _exitDirection;
+            // move to last sampling point
+            RealVector_t currentPosition = _exitPoint - _intersectionLength*_nextStep/2;
+            initFractionals(currentPosition);
+            selectClosestVoxel(_nextStep,currentPosition);
+            _isInAABB = true;
+            _stage=LAST;
+        }
+        break;
+    case LAST:
+