Add rudimentary support for SliceTraversal

docs/modules/projectors.rst

@@ -47,6 +47,13 @@ TraverseAABBJosephsMethod
 -------------------------
 
 .. doxygenclass:: elsa::TraverseAABBJosephsMethod
+    
+SliceTraversal
+------------
+
+.. doxygenclass:: elsa::TransformToTraversal
+
+.. doxygenclass:: elsa::SliceTraversal
 
 SubsetSampler
 -------------------------

elsa/projectors/CMakeLists.txt

@@ -4,6 +4,7 @@ set(MODULE_HEADERS
     Intersection.h
     TraverseAABB.h
     TraverseAABBJosephsMethod.h
+    SliceTraversal.h
     BinaryMethod.h
     SiddonsMethod.h
     JosephsMethod.h
@@ -16,6 +17,7 @@ set(MODULE_SOURCES
     Intersection.cpp
     TraverseAABB.cpp
     TraverseAABBJosephsMethod.cpp
+    SliceTraversal.cpp
     BinaryMethod.cpp
     SiddonsMethod.cpp
     JosephsMethod.cpp

elsa/projectors/Intersection.cpp

 #include "Intersection.h"
+#include "spdlog/fmt/fmt.h"
+#include "spdlog/fmt/ostr.h"
 
 namespace elsa
 {
-    std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
-                                                            const RealRay_t& r)
+    std::tuple<real_t, real_t, real_t, real_t> intersect(const BoundingBox& aabb,
+                                                         const RealRay_t& r)
     {
         real_t invDir = 1 / r.direction()(0);
 
@@ -12,6 +14,8 @@ namespace elsa
 
         real_t tmin = invDir >= 0 ? t1 : t2;
         real_t tmax = invDir >= 0 ? t2 : t1;
+        const auto txmin = tmin;
+        const auto txmax = tmax;
 
         for (int i = 1; i < aabb._min.rows(); ++i) {
             invDir = 1 / r.direction()(i);
@@ -23,10 +27,85 @@ namespace elsa
             tmax = minNum(tmax, invDir >= 0 ? t2 : t1);
         }
 
+        return {tmin, tmax, txmin, txmax};
+    }
+    std::optional<IntersectionResult> Intersection::withRay(const BoundingBox& aabb,
+                                                            const RealRay_t& r)
+    {
+        const auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
+
         if (tmax == 0.0 && tmin == 0.0)
             return std::nullopt;
         if (tmax >= maxNum(tmin, 0.0f)) // hit
             return std::make_optional<IntersectionResult>(tmin, tmax);
         return std::nullopt;
     }
+
+    std::optional<IntersectionResult> Intersection::xPlanesWithRay(BoundingBox aabb,
+                                                                   const RealRay_t& r)
+    {
+        // Slightly increase the size of the bounding box, such that center of voxels are actually
+        // inside the bounding box, parallel rays which directly go through the center of the
+        // voxel, will be recognized with this cheapo hack
+        aabb._min[0] += 0.5f;
+        aabb._max[0] -= 0.5f;
+
+        auto [tmin, tmax, txmin, txmax] = intersect(aabb, r);
+
+        auto dist_to_integer = [](auto real) {
+            const auto aabbmin = static_cast<real_t>(static_cast<int>(std::round(real)));
+            return std::abs(real - aabbmin);
+        };
+
+        auto fractional = [](auto real) {
+            real_t trash = 0;
+            return std::modf(real, &trash);
+        };
+
+        auto advance_to_next_voxel = [&](auto aabb, auto& tmin) -> real_t {
+            // the x-axis coord for voxel centers
+            const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
+
+            RealVector_t entry = r.pointAt(tmin);
+
+            // Calculate the distance from entry.x to the x coord of the next voxel
+            const auto tmp = std::abs(fractional(entry[0] - xvoxelcenter));
+            const auto frac = entry[0] < xvoxelcenter ? tmp : 1 - tmp;
+
+            // Calculate distance in t, but don't advance if it's too small
+            return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
+        };
+
+        auto retreat_to_last_voxel = [&](auto aabb, auto& tmax) -> real_t {
+            const auto xvoxelcenter = dist_to_integer(aabb._min[0]);
+
+            RealVector_t exit = r.pointAt(tmax);
+
+            // calculate the distance from entry.x to the x coord of the previous voxel
+            const auto tmp = std::abs(fractional(exit[0] - xvoxelcenter));
+            const auto frac = exit[0] < xvoxelcenter ? 1 - tmp : tmp;
+
+            // Calculate distance in t, but don't advance if it's too small
+            return tmp > 0.00001 && frac > 0.00001 ? frac / r.direction()[0] : 0;
+        };
+
+        // Sometimes we hit the y axis first, not the x axis, but we don't want
+        // that, we always want the intersection with the x-plane
+        if (txmin < tmin) {
+            tmin += advance_to_next_voxel(aabb, tmin);
+        }
+
+        // This can also happen if we leave at the top of bottom
+        if (txmax > tmax) {
+            tmax -= retreat_to_last_voxel(aabb, tmax);
+        }
+
+        if (tmax == 0.0 && tmin == 0.0)
+            return std::nullopt;
+
+        if (tmax - maxNum(tmin, 0.0f) > -0.000001) // hit
+            return std::make_optional<IntersectionResult>(tmin, tmax);
+
+        return std::nullopt;
+    }
 } // namespace elsa

elsa/projectors/Intersection.h

@@ -3,6 +3,8 @@
 #include "elsaDefines.h"
 #include "BoundingBox.h"
 
+#include "spdlog/fmt/fmt.h"
+
 #include <Eigen/Geometry>
 #include <limits>
 #include <optional>
@@ -57,6 +59,9 @@ namespace elsa
          */
         static std::optional<IntersectionResult> withRay(const BoundingBox& aabb,
                                                          const RealRay_t& r);
+
+        static std::optional<IntersectionResult> xPlanesWithRay(BoundingBox aabb,
+                                                                const RealRay_t& r);
     };
 
     /**
@@ -92,3 +97,18 @@ namespace elsa
     }
 
 } // namespace elsa
+
+template <>
+struct fmt::formatter<elsa::IntersectionResult> {
+    template <typename ParseContext>
+    constexpr auto parse(ParseContext& ctx)
+    {
+        return ctx.begin();
+    }
+
+    template <typename FormatContext>
+    auto format(const elsa::IntersectionResult& hit, FormatContext& ctx)
+    {
+        return fmt::format_to(ctx.out(), "{{ tmin: {}, tmax: {} }}", hit._tmin, hit._tmax);
+    }
+};

elsa/projectors/SliceTraversal.cpp

+#include "SliceTraversal.h"
+#include "elsaDefines.h"
+#include "spdlog/fmt/fmt.h"
+#include "spdlog/fmt/ostr.h"
+
+namespace elsa
+{
+    RealMatrix_t TransformToTraversal::create2DRotation(const real_t leadingCoeff,
+                                                        const index_t leadingAxisIndex)
+    {
+        auto createRotationMatrix = [](real_t radian) {
+            real_t c = std::cos(radian);
+            real_t s = std::sin(radian);
+            return RealMatrix_t({{c, -s}, {s, c}});
+        };
+
+        // TODO: Does a closed form solution exist?
+        // Use conversion to polar coordinates
+        if (leadingAxisIndex == 0 && leadingCoeff >= 0) {
+            // Already identity, so do nothing
+            return createRotationMatrix(0);
+        } else if (leadingAxisIndex == 0 && leadingCoeff < 0) {
+            // create a 2D 180° rotation matrix
+            return createRotationMatrix(pi_t);
+        } else if (leadingAxisIndex == 1 && leadingCoeff >= 0) {
+            // create a 2D 270° rotation matrix counter clockwise
+            return createRotationMatrix(3 * pi_t / 2);
+        } else if (leadingAxisIndex == 1 && leadingCoeff <= 0) {
+            // create a 2D 90° rotation matrix counter clockwise
+            return createRotationMatrix(pi_t / 2);
+        }
+        return createRotationMatrix(0);
+    }
+
+    RealMatrix_t TransformToTraversal::createRotation(RealRay_t ray)
+    {
+        // Get the leading axis, by absolute value
+        index_t leadingAxisIndex = 0;
+        ray.direction().array().abs().maxCoeff(&leadingAxisIndex);
+
+        // Get the signed leading coefficient
+        const auto leadingCoeff = ray.direction()(leadingAxisIndex);
+
+        return create2DRotation(leadingCoeff, leadingAxisIndex);
+    }
+
+    RealMatrix_t TransformToTraversal::createTransformation(const RealRay_t& ray,
+                                                            const RealVector_t& centerOfRotation)
+    {
+        const auto dim = ray.dim();
+
+        // Setup translation
+        RealMatrix_t translation(dim + 1, dim + 1);
+        translation.setIdentity();
+        translation.block(0, dim, dim, 1) = -centerOfRotation;
+
+        // Setup rotation
+        RealMatrix_t rotation(dim + 1, dim + 1);
+        rotation.setIdentity();
+        rotation.block(0, 0, dim, dim) = createRotation(ray);
+
+        return rotation * translation;
+    }
+
+    TransformToTraversal::TransformToTraversal(const RealRay_t& ray,
+                                               const RealVector_t& centerOfRotation)
+        : transformation_(createTransformation(ray, centerOfRotation)),
+          translation_(-centerOfRotation)
+    {
+    }
+
+    RealRay_t TransformToTraversal::toTraversalCoordinates(RealRay_t ray) const
+    {
+        ray.origin() = *this * Point(ray.origin());
+        ray.direction() = *this * Vec(ray.direction());
+        return ray;
+    }
+
+    BoundingBox TransformToTraversal::toTraversalCoordinates(BoundingBox aabb) const
+    {
+        // Only translate, as the rotations are always 90, 180, 270 degrees,
+        // TODO: this might be wrong, idk, what about non square things
+        aabb._min = *this * Point(aabb._min);
+        aabb._max = *this * Point(aabb._max);
+        aabb.recomputeBounds();
+        return aabb;
+    }
+
+    RealMatrix_t TransformToTraversal::transformation() const { return transformation_; }
+
+    RealMatrix_t TransformToTraversal::invTransformation() const
+    {
+        return transformation().inverse();
+    }
+
+    RealMatrix_t TransformToTraversal::rotation() const
+    {
+        const auto dim = transformation_.rows();
+        return transformation_.block(0, 0, dim - 1, dim - 1);
+    }
+
+    RealMatrix_t TransformToTraversal::invRotation() const { return rotation().transpose(); }
+
+    RealVector_t TransformToTraversal::translation() const { return translation_; }
+
+    RealMatrix_t TransformToTraversal::linear() const { return rotation(); }
+
+    RealVector_t TransformToTraversal::operator*(const Point<real_t>& point) const
+    {
+        return (transformation() * point.point_.homogeneous()).hnormalized();
+    }
+
+    RealVector_t TransformToTraversal::operator*(const Vec<real_t>& vec) const
+    {
+        return linear() * vec.vec_;
+    }
+
+    SliceTraversal::SliceTraversal(BoundingBox aabb, RealRay_t ray)
+        : transformation_(ray, aabb._max.array() / 2), ray_(ray)
+    {
+        Eigen::IOFormat vecfmt(10, 0, ", ", ", ", "", "", "[", "]");
+        // {
+        //     auto tmphit = Intersection::xPlanesWithRay(aabb, ray);
+        //
+        //     fmt::print("aabb: {}\n", aabb);
+        //     fmt::print("ro: {}\n", ray.origin().format(vecfmt));
+        //     fmt::print("rd: {}\n", ray.direction().format(vecfmt));
+        //     if (tmphit) {
+        //         fmt::print("hit: {}\n", *tmphit);
+        //         fmt::print("diff: {}\n", tmphit->_tmax - tmphit->_tmin);
+        //         fmt::print("firstVoxel: {}\n", ray.pointAt(tmphit->_tmin).format(vecfmt));
+        //         fmt::print("lastVoxel: {}\n", ray.pointAt(tmphit->_tmax).format(vecfmt));
+        //     }
+        // }
+
+        // Transform ray and aabb to traversal coordinate space
+        ray = transformation_.toTraversalCoordinates(ray);
+        aabb = transformation_.toTraversalCoordinates(aabb);
+        // fmt::print("aabb: {}\n", aabb);
+
+        // TODO: We only want to shift in x direction by 0.5, not in y
+        auto hit = Intersection::xPlanesWithRay(aabb, ray);
+
+        // Default init is sufficient if we don't hit,  TODO: test that!
+        if (hit) {
+            const auto firstVoxel = ray.pointAt(hit->_tmin);
+            const auto lastVoxel = ray.pointAt(hit->_tmax);
+
+            // fmt::print("ro: {}\n", ray.origin().format(vecfmt));
+            // fmt::print("rd: {}\n", ray.direction().format(vecfmt));
+            // fmt::print("hit: {}\n", *hit);
+            // fmt::print("diff: {}\n", hit->_tmax - hit->_tmin);
+            // fmt::print("firstVoxel: {}\n", firstVoxel.format(vecfmt));
+            // fmt::print("lastVoxel: {}\n", lastVoxel.format(vecfmt));
+
+            endIndex_ = std::max<index_t>(
+                1, static_cast<index_t>(std::ceil(lastVoxel[0] - firstVoxel[0] + 0.5f)));
+
+            tDelta_ = 1 / ray.direction()[0];
+
+            t_ = hit->_tmin;
+        }
+    }
+
+    SliceTraversal::Iter SliceTraversal::begin() const { return {startIndex_, ray_, tDelta_, t_}; }
+
+    SliceTraversal::Iter SliceTraversal::end() const { return {endIndex_, ray_, tDelta_, t_}; }
+
+    index_t SliceTraversal::startIndex() const { return startIndex_; }
+
+    index_t SliceTraversal::endIndex() const { return endIndex_; }
+
+    SliceTraversal::Iter::value_type SliceTraversal::Iter::operator*() const
+    {
+        const RealVector_t curPos = ray_.pointAt(t_);
+        const IndexVector_t curVoxel = curPos.template cast<index_t>();
+        return {curPos, curVoxel, t_};
+    }
+
+    SliceTraversal::Iter& SliceTraversal::Iter::operator++()
+    {
+        ++pos_;
+        t_ += tDelta_;
+        return *this;
+    }
+
+    SliceTraversal::Iter SliceTraversal::Iter::operator++(int)
+    {
+        auto copy = *this;
+        ++(*this);
+        return copy;
+    }
+
+    bool operator==(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
+    {
+        return lhs.pos_ == rhs.pos_;
+    }
+
+    bool operator!=(const SliceTraversal::Iter& lhs, const SliceTraversal::Iter& rhs)
+    {
+        return !(lhs == rhs);
+    }
+
+} // namespace elsa

elsa/projectors/SliceTraversal.h

+#pragma once
+
+#include "Intersection.h"
+#include "elsaDefines.h"
+#include "BoundingBox.h"
+#include "Logger.h"
+#include "Assertions.h"
+
+#include <Eigen/Geometry>
+#include <cassert>
+#include <optional>
+#include <cmath>
+#include <utility>
+#include <variant>
+
+// TODO: remove
+#include <iostream>
+
+namespace elsa
+{
+    /// Strong type to distinguish transformation of points vs transformation of vectors
+    template <class data_t>
+    struct Point {
+        explicit Point(Vector_t<data_t> point) : point_(std::move(point)) {}
+
+        Vector_t<data_t> point_;
+    };
+
+    /// Strong type to distinguish transformation of points vs transformation of vectors
+    template <class data_t>
+    struct Vec {
+        explicit Vec(Vector_t<data_t> vec) : vec_(std::move(vec)) {}
+
+        Vector_t<data_t> vec_;
+    };
+
+    /**
+     * @brief Represent a transformation, which will transform any point into the
+     * traversal coordinate system.
+     *
+     * In the traversal coordinate system, the reference direction - form which it is constructed -
+     * is transformed, such that the leading coefficient (i.e. the coefficient of largest absolute
+     * value), is in the positive x-direction (i.e. the first component)
+     *
+     * The algorithm to determine the rotation in a 2D case is quite simple: Given the vector
+     * \f$\begin{bmatrix} x & y \end{bmatrix}\f$,
+     * first determine the leading axis, by finding the index of the component with maximum
+     * absolute value, i.e. \f$\text{axis} = \text{maxAxis} (\lvert x \rvert, \lvert y \rvert) \f$
+     * (in this case, the max returns "x", or "y" for our use cases instead of the actual value).
+     * Assume `axis` is either `"x"` or `"y"` and `leadingCoeff` stores corresponding value
+     * to the leading axis (but it's signed!). Then in pseudocode the algorithm determines the
+     * rotation the following way:
+     *
+     * ```python
+     * if axis == "x" and maxCoeff >= 0:
+     *     return rotate_by(0)
+     * elif axis == "x" and maxCoeff <= 0:
+     *     return rotate_by(180)
+     * elif axis == "y" and maxCoeff >= 0:
+     *     return rotate_by(90)
+     * elif axis == "y" and maxCoeff <= 0:
+     *     return rotate_by(270)
+     * ```
+     *
+     * To extent it to 3D one further decision has to be made, do we only rotate, such that the
+     * leading direction is in the x direction, or do we rotate, such that y is the second largest
+     * direction and z is the smallest.
+     *
+     * TODO: It might be nice to move this to a separate file and module, to wrap some of Eigens
+     * transformation stuff, such that it's usable for dynamic cases, as we have here. This
+     * might be nice to have for the Geometry class as well, but not for now.
+     */
+    class TransformToTraversal
+    {
+    private:
+        static RealMatrix_t create2DRotation(const real_t leadingCoeff,
+                                             const index_t leadingAxisIndex);
+
+        static RealMatrix_t createRotation(RealRay_t ray);
+
+        static RealMatrix_t createTransformation(const RealRay_t& ray,
+                                                 const RealVector_t& centerOfRotation);
+
+        RealMatrix_t transformation_;
+        RealVector_t translation_;
+
+    public:
+        TransformToTraversal(const RealRay_t& ray, const RealVector_t& centerOfRotation);
+
+        RealRay_t toTraversalCoordinates(RealRay_t ray) const;
+
+        BoundingBox toTraversalCoordinates(BoundingBox aabb) const;
+
+        RealMatrix_t transformation() const;
+
+        RealMatrix_t invTransformation() const;
+
+        RealMatrix_t rotation() const;
+
+        RealMatrix_t invRotation() const;
+
+        RealVector_t translation() const;
+
+        RealMatrix_t linear() const;
+
+        RealVector_t operator*(const Point<real_t>& point) const;
+
+        RealVector_t operator*(const Vec<real_t>& vec) const;
+    };
+
+    /**
+     * @brief Traverse a volume along the direction of a given ray. Each step it will advance one
+     * slice further in the direction of the ray. The traversal visits voxels at the center planes,
+     * i.e. it always evaluates at the center of voxel in the plane of leading direction.
+     *
+     * This is a slightly modified version of the algorithm of Amanatides & Woo's "A Fast Voxel
+     * Traversal Algorithm". The reference visits each and every single voxel along the ray.
+     * However, for this case we don't need that.
+     *
+     * Given a bounding box of the volume and a ray to traverse the volume, one can simply call:
+     *
+     * ```cpp
+     * for(auto [pos, voxel, t] = SliceTraversal(aabb, ray)) {
+     *     // pos: exact position on center plane of voxel
+     *     // voxel: current voxel
+     *     // t: position on the ray (i.e. pos = ray.origin() + t * ray.direction())
+     * }
+     * ```
+     *
+     * The number of visited voxels is known at construction, therefore moving along the ray is
+     * basically decrementing a count. However, to return useful information, more bookkeeping is
+     * needed.
+     *
+     * Dereferencing the iterator will return a small struct, which has information about the
+     * current position in the volume, the current voxel and the t parameter of the ray.
+     *
+     * Note: When using voxel returned from iterating over the volume, it can happen that voxel
+     * on the boundary of the volume are returned. This can trigger wrong behaviour for certain
+     * applications and therefore should be handled with care.
+     *
+     * A note about the implementation: To ease the construction and handling, the incoming ray and
+     * bounding box are transformed to a traversal internal coordinate space, in which the leading
+     * direction lies in the positive x-axis. The 'world' space is partitioned in 4 quadrants based
+     * on the direction. The first ranges for all direction from  with a polar angle in the range of
+     * [-45°, 45°], the second (45°, 135°), the third [135°, 225°], and the fourth (225°, 315°).
+     * Note that -45° = 315°.
+     *
+     * TODO:
+     * - Make the iterator random access
+     *
+     * @author
+     * - David Frank - initial code
+     */
+    class SliceTraversal
+    {
+    private:
+        TransformToTraversal transformation_;
+
+        /// TODO: Do we want to consider points that spawn inside the volume?
+        index_t startIndex_{0};
+
+        index_t endIndex_{0};
+
+        /// t value for the first intersection of the ray and the bounding box
+        real_t t_{0.0};
+
+        /// t value to increment each iteration
+        real_t tDelta_{0.0};
+
+        RealRay_t ray_;
+
+    public:
+        /// The Dereference type of the iterator
+        /// TODO: With C++20 support for proxy iterators is given, maybe this could change then
+        struct IterValue {
+            RealVector_t curPosition_;
+            IndexVector_t curVoxel_;
+            real_t t_;
+        };
+
+        /// Traversal iterator, models forward iterator, maybe this should actually be an input
+        /// iterator due to the non reference type of the dereference type. IDK, as we use it, this
+        /// works, but in the future this might should be different.
+        struct Iter {
+        public:
+            using iterator_category = std::forward_iterator_tag;
+            using difference_type = std::ptrdiff_t;
+            using value_type = IterValue;
+            using pointer = value_type*;
+            using reference = value_type&;
+
+        private:
+            index_t pos_;
+            RealRay_t ray_;
+            real_t tDelta_;
+            real_t t_;
+
+        public:
+            /// Construct iterator
+            ///
+            /// @param pos index position of the traversal, used mainly for comparing two iterators
+            /// @param ray traversed ray used to compute exact position on dereference
+            /// @param deltat increment of t each increment
+            /// @param t position along the ray
+            Iter(index_t pos, RealRay_t ray, real_t deltat, real_t t)
+                : pos_(pos), ray_(ray), tDelta_(deltat), t_(t)
+            {
+            }
+
+            /// Dereference iterator
+            value_type operator*() const;
+
+            /// Pre increment
+            Iter& operator++();
+
+            /// Post increment
+            Iter operator++(int);
+
+            /// Equal to operator with iterator
+            friend bool operator==(const Iter& lhs, const Iter& rhs);
+
+            /// Not equal to operator with iterator
+            friend bool operator!=(const Iter& lhs, const Iter& rhs);
+        };
+
+        // Delete default construction
+        SliceTraversal() = delete;
+
+        // Constr