Add BSpline projector

elsa/core/Utilities/Math.hpp

@@ -4,6 +4,48 @@
 
 namespace elsa
 {
+    namespace math
+    {
+        /// Compute factorial \f$n!\f$ recursively
+        constexpr inline index_t factorial(index_t n) noexcept
+        {
+            return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
+        }
+
+        /// Compute binomial coefficient
+        constexpr inline index_t binom(index_t n, index_t k) noexcept
+        {
+            return (k > n)
+                       ? 0
+                       : (k == 0 || k == n) ? 1
+                                            : (k == 1 || k == n - 1)
+                                                  ? n
+                                                  : (k + k < n) ? (binom(n - 1, k - 1) * n) / k
+                                                                : (binom(n - 1, k) * n) / (n - k);
+        }
+
+        /// Compute Heaviside-function
+        /// \f[
+        /// x \mapsto
+        /// \begin{cases}
+        /// 0: & x < 0 \\
+        /// c: & x = 0 \\
+        /// 1: & x > 0
+        /// \end{cases}
+        /// \f]
+        template <typename data_t>
+        constexpr data_t heaviside(data_t x1, data_t c)
+        {
+            if (x1 == 0) {
+                return c;
+            } else if (x1 < 0) {
+                return 0;
+            } else {
+                return 1;
+            }
+        }
+    } // namespace math
+
     /// proposed in Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution
     /// Equations. AMS, 2001
     template <typename data_t>

elsa/core/tests/CMakeLists.txt

@@ -29,3 +29,4 @@ ELSA_DOCTEST(DataContainer)
 ELSA_DOCTEST(DataContainerFormatter)
 ELSA_DOCTEST(CartesianIndices)
 ELSA_DOCTEST(Bessel)
+ELSA_DOCTEST(Math)

elsa/core/tests/test_Math.cpp

+#include "doctest/doctest.h"
+
+#include "Math.hpp"
+
+TEST_SUITE_BEGIN("Math");
+
+using namespace elsa;
+
+TEST_CASE("Math::factorial")
+{
+    CHECK_EQ(1, math::factorial(0));
+    CHECK_EQ(1, math::factorial(1));
+
+    auto fac = 1;
+    for (int i = 1; i < 10; ++i) {
+        fac *= i;
+        CHECK_EQ(fac, math::factorial(i));
+    }
+}
+
+TEST_CASE("Math::binom")
+{
+    GIVEN("n == 10")
+    {
+        const index_t n = 10;
+
+        WHEN("k larger than n") { CHECK_EQ(math::binom(n, 15), 0); }
+
+        WHEN("k == 0 or k == n")
+        {
+            CHECK_EQ(math::binom(n, 0), 1);
+            CHECK_EQ(math::binom(n, 10), 1);
+        }
+
+        CHECK_EQ(math::binom(n, 1), 10);
+        CHECK_EQ(math::binom(n, 2), 45);
+        CHECK_EQ(math::binom(n, 3), 120);
+        CHECK_EQ(math::binom(n, 4), 210);
+        CHECK_EQ(math::binom(n, 5), 252);
+        CHECK_EQ(math::binom(n, 6), 210);
+        CHECK_EQ(math::binom(n, 7), 120);
+        CHECK_EQ(math::binom(n, 8), 45);
+        CHECK_EQ(math::binom(n, 9), 10);
+    }
+}
+
+TEST_CASE("Math::heaviside")
+{
+    constexpr index_t size = 200;
+    constexpr real_t c = 0.5;
+    const auto linspace = Vector_t<real_t>::LinSpaced(size, -2, 2);
+
+    for (std::size_t i = 0; i < size; ++i) {
+        auto res = math::heaviside(linspace[i], c);
+        if (linspace[i] == 0.) {
+            CHECK_EQ(res, c);
+        } else if (linspace[i] < 0) {
+            CHECK_EQ(res, 0);
+        } else if (linspace[i] > 0) {
+            CHECK_EQ(res, 1);
+        }
+    }
+}
+
+TEST_SUITE_END();

elsa/projectors/BSplines.cpp

+#include "BSplines.h"
+
+namespace elsa
+{
+    template <typename data_t>
+    BSpline<data_t>::BSpline(index_t dim, index_t order) : dim_(dim), order_(order)
+    {
+    }
+
+    template <typename data_t>
+    data_t BSpline<data_t>::operator()(Vector_t<data_t> x)
+    {
+        ELSA_VERIFY(dim_ == x.size());
+        return bspline::nd_bspline_evaluate(x, order_);
+    }
+
+    template <typename data_t>
+    index_t BSpline<data_t>::order() const
+    {
+        return order_;
+    }
+
+    template <typename data_t>
+    index_t BSpline<data_t>::dim() const
+    {
+        return dim_;
+    }
+
+    template <typename data_t>
+
+    ProjectedBSpline<data_t>::ProjectedBSpline(index_t dim, index_t order)
+        : dim_(dim), order_(order)
+    {
+    }
+
+    template <typename data_t>
+    data_t ProjectedBSpline<data_t>::operator()(data_t x)
+    {
+        return bspline::nd_bspline_centered(x, order_, dim_ - 1);
+    }
+
+    template <typename data_t>
+    index_t ProjectedBSpline<data_t>::order() const
+    {
+        return order_;
+    }
+
+    template <typename data_t>
+    index_t ProjectedBSpline<data_t>::dim() const
+    {
+        return dim_;
+    }
+
+    // ------------------------------------------
+    // explicit template instantiation
+    template class BSpline<float>;
+    template class BSpline<double>;
+
+    template class ProjectedBSpline<float>;
+    template class ProjectedBSpline<double>;
+} // namespace elsa

elsa/projectors/BSplines.h

+#pragma once
+
+#include "elsaDefines.h"
+#include "Assertions.h"
+#include "Math.hpp"
+
+#include <cmath>
+#include "spdlog/fmt/bundled/core.h"
+#include "spdlog/fmt/fmt.h"
+
+namespace elsa
+{
+    namespace bspline
+    {
+        /// @brief Evaluate the 1-dimensional B-Spline of degree m, with m + 2 equally spaced knots
+        ///
+        /// This is based on the implementation given as in e.g.:
+        /// Fast B-spline transforms for continuous image representation and interpolation, Unser
+        /// et. al. (1991), equation 2.2
+        ///
+        /// @param x coordinate to evaluate 1-dimensional B-Spline at
+        /// @param m order of B-Spline
+        template <typename data_t>
+        constexpr data_t bspline1d_evaluate(data_t x, index_t m) noexcept
+        {
+            const auto xs = x + (m + 1) / data_t{2};
+            data_t res = 0;
+
+            for (int n = 0; n <= m + 1; ++n) {
+                const auto tmp1 = math::heaviside<data_t>(xs - n, 0);
+                const auto tmp2 = std::pow<data_t>(xs - n, m);
+                const auto tmp3 = math::binom(m + 1, n);
+                const auto tmp4 = std::pow<data_t>(-1, n) / math::factorial(m);
+                res += tmp1 * tmp2 * tmp3 * tmp4;
+            }
+
+            return res;
+        }
+
+        /// @brief Evaluate n-dimensional B-Spline of degree m. As B-Splines are separable, this is
+        /// just the product of 1-dimensional B-Splines.
+        ///
+        /// @param x n-dimensional coordinate to evaluate B-Spline at
+        /// @param m order of B-Spline
+        template <typename data_t>
+        constexpr data_t nd_bspline_evaluate(const Vector_t<data_t>& x, index_t m) noexcept
+        {
+            data_t res = bspline1d_evaluate(x[0], m);
+            for (int i = 1; i < x.size(); ++i) {
+                res *= bspline1d_evaluate(x[i], m);
+            }
+            return res;
+        }
+
+        /// @brief Evaluate n-dimensional B-Spline at a given coordinate for the first dimension,
+        /// and at the center (i.e. `0`) at all the other dimensions. This is particular useful as
+        /// an approximation during the calculation of the line integral
+        ///
+        /// @param x 1-dimensional coordinate to evaluate B-Spline at
+        /// @param m order of B-Spline
+        /// @param dim dimension of B-Spline
+        template <typename data_t>
+        constexpr data_t nd_bspline_centered(data_t x, int m, int dim) noexcept
+        {
+            data_t res = bspline1d_evaluate<data_t>(x, m);
+            for (int i = 1; i < dim; ++i) {
+                const auto inc = bspline1d_evaluate<data_t>(0., m);
+                res *= inc;
+            }
+            return res;
+        }
+    } // namespace bspline
+
+    /// @brief Represent a B-Spline basis function of a given dimension and order
+    template <typename data_t>
+    class BSpline
+    {
+    public:
+        BSpline(index_t dim, index_t order);
+
+        data_t operator()(Vector_t<data_t> x);
+
+        index_t order() const;
+
+        index_t dim() const;
+
+    private:
+        /// Dimension of B-Spline
+        index_t dim_;
+
+        /// Order of B-Spline
+        index_t order_;
+    };
+
+    /// @brief Represent a projected B-Spline basis function of a given dimension and order.
+    /// Projected B-Splines are again B-Spline of n-1 dimensions. Using the fact, that B-Splines
+    /// are close to symmetrical, we can approximate the projection only based on distance.
+    template <typename data_t>
+    class ProjectedBSpline
+    {
+    public:
+        ProjectedBSpline(index_t dim, index_t order);
+
+        data_t operator()(data_t s);
+
+        index_t order() const;
+
+        index_t dim() const;
+
+    private:
+        /// Dimension of B-Spline
+        index_t dim_;
+
+        /// Order of B-Spline
+        index_t order_;
+    };
+} // namespace elsa

elsa/projectors/CMakeLists.txt

@@ -11,6 +11,8 @@ set(MODULE_HEADERS
     JosephsMethod.h
     SubsetSampler.h
     Blobs.h
+    BSplines.h
+    Luts.hpp
 )
 
 # list all the code files of the module
@@ -26,6 +28,7 @@ set(MODULE_SOURCES
     JosephsMethod.cpp
     SubsetSampler.cpp
     Blobs.cpp
+    BSplines.cpp
 )
 
 list(APPEND MODULE_PUBLIC_DEPS elsa_core elsa_logging)

elsa/projectors/LutProjector.cpp

@@ -34,8 +34,69 @@ namespace elsa
     {
     }
 
+    template <typename data_t>
+    BSplineProjector<data_t>::BSplineProjector(data_t degree,
+                                               const VolumeDescriptor& domainDescriptor,
+                                               const DetectorDescriptor& rangeDescriptor)
+        : LutProjector<data_t, BSplineProjector<data_t>>(domainDescriptor, rangeDescriptor),
+          lut_(domainDescriptor.getNumberOfDimensions(), degree)
+    {
+        // sanity checks
+        auto dim = domainDescriptor.getNumberOfDimensions();
+        if (dim < 2 || dim > 3) {
+            throw InvalidArgumentError("BSplineProjector: only supporting 2d/3d operations");
+        }
+
+        if (dim != rangeDescriptor.getNumberOfDimensions()) {
+            throw InvalidArgumentError(
+                "BSplineProjector: domain and range dimension need to match");
+        }
+
+        if (_detectorDescriptor.getNumberOfGeometryPoses() == 0) {
+            throw InvalidArgumentError("BSplineProjector: rangeDescriptor without any geometry");
+        }
+    }
+
+    template <typename data_t>
+    BSplineProjector<data_t>::BSplineProjector(const VolumeDescriptor& domainDescriptor,
+                                               const DetectorDescriptor& rangeDescriptor)
+        : BSplineProjector(2, domainDescriptor, rangeDescriptor)
+    {
+    }
+
+    template <typename data_t>
+    data_t BSplineProjector<data_t>::weight(data_t distance) const
+    {
+        return lut_(distance);
+    }
+
+    template <typename data_t>
+    index_t BSplineProjector<data_t>::support() const
+    {
+        return static_cast<index_t>(2);
+    }
+
+    template <typename data_t>
+    BSplineProjector<data_t>* BSplineProjector<data_t>::cloneImpl() const
+    {
+        return new BSplineProjector(_volumeDescriptor, _detectorDescriptor);
+    }
+
+    template <typename data_t>
+    bool BSplineProjector<data_t>::isEqual(const LinearOperator<data_t>& other) const
+    {
+        if (!Base::isEqual(other))
+            return false;
+
+        auto otherOp = downcast_safe<BSplineProjector>(&other);
+        return static_cast<bool>(otherOp);
+    }
+
     // ------------------------------------------
     // explicit template instantiation
     template class BlobProjector<float>;
     template class BlobProjector<double>;
+
+    template class BSplineProjector<float>;
+    template class BSplineProjector<double>;
 } // namespace elsa

elsa/projectors/LutProjector.h

@@ -233,4 +233,32 @@ namespace elsa
         using Base::_detectorDescriptor;
         using Base::_volumeDescriptor;
     };
+
+    template <typename data_t = real_t>
+    class BSplineProjector : public LutProjector<data_t, BSplineProjector<data_t>>
+    {
+    public:
+        BSplineProjector(data_t degree, const VolumeDescriptor& domainDescriptor,
+                         const DetectorDescriptor& rangeDescriptor);
+
+        BSplineProjector(const VolumeDescriptor& domainDescriptor,
+                         const DetectorDescriptor& rangeDescriptor);
+
+        data_t weight(data_t distance) const;
+
+        index_t support() const;
+
+        /// implement the polymorphic clone operation
+        BSplineProjector<data_t>* cloneImpl() const override;
+
+        /// implement the polymorphic comparison operation
+        bool isEqual(const LinearOperator<data_t>& other) const override;
+
+    private:
+        ProjectedBSplineLut<data_t, 100> lut_;
+
+        using Base = LutProjector<data_t, BSplineProjector<data_t>>;
+        using Base::_detectorDescriptor;
+        using Base::_volumeDescriptor;
+    };
 } // namespace elsa

elsa/projectors/Luts.hpp

 #pragma once
 
 #include "Blobs.h"
+#include "BSplines.h"
 #include "Logger.h"
 #include "Timer.h"
 
@@ -28,6 +29,24 @@ namespace elsa
             return lut;
         }
 
+        template <typename data_t, index_t N>
+        constexpr std::array<data_t, N> bspline_lut(ProjectedBSpline<data_t> bspline)
+        {
+            Logger::get("bspline_lut")->debug("Calculating lut");
+
+            std::array<data_t, N> lut;
+
+            auto t = static_cast<data_t>(0);
+            const auto step = 2. / N;
+
+            for (std::size_t i = 0; i < N; ++i) {
+                lut[i] = bspline(t);
+                t += step;
+            }
+
+            return lut;
+        }
+
         template <typename data_t>
         data_t lerp(data_t a, SelfType_t<data_t> b, SelfType_t<data_t> t)
         {
@@ -119,4 +138,25 @@ namespace elsa
         ProjectedBlob<data_t> blob_;
         Lut<data_t, N> lut_;
     };
+
+    template <typename data_t, index_t N>
+    class ProjectedBSplineLut
+    {
+    public:
+        constexpr ProjectedBSplineLut(int dim, int degree)
+            : bspline_(dim, degree), lut_(detail::bspline_lut<data_t, N>(bspline_))
+        {
+        }
+
+        constexpr data_t order() const { return bspline_.order(); }
+
+        constexpr data_t operator()(data_t distance) const
+        {
+            return lut_((std::abs(distance) / 2.) * N);
+        }
+
+    private:
+        ProjectedBSpline<data_t> bspline_;
+        Lut<data_t, N> lut_;
+    };
 } // namespace elsa

elsa/projectors/tests/CMakeLists.txt

@@ -20,6 +20,7 @@ ELSA_DOCTEST(TraverseAABB)
 ELSA_DOCTEST(SliceTraversal)
 ELSA_DOCTEST(LutProjector)
 ELSA_DOCTEST(Blobs)
+ELSA_DOCTEST(BSplines)
 ELSA_DOCTEST(Luts)
 ELSA_DOCTEST(BinaryMethod)
 ELSA_DOCTEST(SiddonsMethod)

elsa/projectors/tests/test_BSplines.cpp

+#include "doctest/doctest.h"
+
+#include "BSplines.h"
+
+#include "spdlog/fmt/fmt.h"
+#include "spdlog/fmt/ostr.h"
+#include "spdlog/fmt/bundled/ranges.h"
+#include <utility>
+
+using namespace elsa;
+
+TEST_SUITE_BEGIN("projectors::BSplines");
+
+Eigen::IOFormat vecfmt(10, 0, ", ", ", ", "", "", "[", "]");
+
+// TEST_CASE_TEMPLATE("BSplines: 1d BSpline evaluation", data_t, float, double)
+// TEST_CASE_TEMPLATE("BSplines: 1d BSpline evaluation", data_t, float)
+// {
+//     constexpr auto size = 11;
+//
+//     const auto low = -2;
+//     const auto high = 2;
+//
+//     const auto linspace = Vector_t<data_t>::LinSpaced(size, low, high);
+//     std::array<data_t, size> spline;
+//
+//     auto degree = 3;
+//     for (int i = 0; i < size; ++i) {
+//         for (int j = 0; j < size; ++j) {
+//             RealVector_t coord({{linspace[i], linspace[j]}});
+//             // fmt::print("bspline({}, {}) = {}\n", coord.format(vecfmt), degree,
+//             //            bspline::nd_bspline_evaluate(coord, degree));
+//         }
+//     }
+//
+//     // fmt::print("{}\n", linspace.format(vecfmt));
+//     // fmt::print("{}\n", spline);
+// }
+
+/// Quick and easy simpsons rule for numerical integration
+template <typename data_t, typename Func>
+data_t simpsons(Func f, data_t a, data_t b, int N = 50)
+{
+    const auto h = (b - a) / N;
+
+    data_t sum_odds = 0.0;
+    for (int i = 1; i < N; i += 2) {
+        sum_odds += f(a + i * h);
+    }
+
+    data_t sum_evens = 0.0;
+    for (int i = 2; i < N; i += 2) {
+        sum_evens += f(a + i * h);
+    }
+
+    return (f(a) + f(b) + 2 * sum_evens + 4 * sum_odds) * h / 3;
+}
+
+/// Integrate from distance s from center, to a
+/// see Fig 3.6 (p. 53) and Listing 3.1 (p. 58) in Three-Dimensional Digital Tomosynthesis by
+/// Levakhina
+template <typename data_t>
+data_t bspline_line_integral(data_t s, index_t m, index_t dim)
+{
+    auto x = [=](auto t) {
+        data_t res = bspline::bspline1d_evaluate<data_t>(std::sqrt(t * t + s * s), m);
+        for (int i = 1; i < dim; ++i) {
+            res *= bspline::bspline1d_evaluate<data_t>(0., m);
+        }
+        return res;
+    };
+
+    return 2 * simpsons<data_t>(x, 0, 3);
+}
+
+TEST_CASE_TEMPLATE("BSplines: 1d line integal", data_t, float)
+{
+    constexpr auto size = 21;
+
+    const auto low = -2;
+    const auto high = 2;
+
+    const auto linspace = Vector_t<data_t>::LinSpaced(size, low, high);
+
+    const int degree = 2;
+
+    // const data_t distance = 1.f;
+    //
+    // fmt::print("1D bspline at distance {:4.2f}: {:8.5f}\n", distance,
+    //            bspline::bspline1d_evaluate<data_t>(distance, degree));
+    // fmt::print("2D bspline at distance {:4.2f}: {:8.5f}\n", distance,
+    //            bspline::bspline1d_evaluate<data_t>(0, degree)
+    //                * bspline::bspline1d_evaluate<data_t>(distance, degree));
+    // fmt::print("2D line integral: {:8.5f}\n", bspline_line_integral<data_t>(distance, degree,
+    // 2)); fmt::print("3D line integral: {:8.5f}\n", bspline_line_integral<data_t>(distance,
+    // degree, 3)); MESSAGE("helloo");
+
+    // BSpline<data_t> bspline_1d(1, degree);
+    // BSpline<data_t> bspline_2d(2, degree);
+    //
+    // for (int i = 0; i < 101; ++i) {
+    //     const data_t x = (i / 25.) - 2.;
+    //
+    //     CAPTURE(x);
+    //     CAPTURE(bspline::bspline1d_evaluate(x, degree));
+    //     CAPTURE(bspline::bspline1d_evaluate(0., degree));
+    //
+    //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree),
+    //              doctest::Approx(bspline_1d(Vector_t<data_t>({{x}}))));
+    //
+    //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(0.,
+    //     degree),
+    //              doctest::Approx(bspline_2d(Vector_t<data_t>({{x, 0}}))));
+    //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(x, degree),
+    //              doctest::Approx(bspline_2d(Vector_t<data_t>({{x, x}}))));
+    // }
+
+    ProjectedBSpline<data_t> proj_bspline_2d(2, degree);
+    ProjectedBSpline<data_t> proj_bspline_3d(3, degree);
+
+    CHECK_EQ(proj_bspline_2d.order(), degree);
+    CHECK_EQ(proj_bspline_3d.order(), degree);
+    CHECK_EQ(proj_bspline_2d.dim(), 2);