Line data    Source code

       1             : #include "doctest/doctest.h"
       2             : 
       3             : #include "BSplines.h"
       4             : 
       5             : #include "spdlog/fmt/fmt.h"
       6             : #include "spdlog/fmt/ostr.h"
       7             : #include "spdlog/fmt/bundled/ranges.h"
       8             : #include <utility>
       9             : 
      10             : using namespace elsa;
      11             : 
      12             : TEST_SUITE_BEGIN("projectors::BSplines");
      13             : 
      14             : Eigen::IOFormat vecfmt(10, 0, ", ", ", ", "", "", "[", "]");
      15             : 
      16             : // TEST_CASE_TEMPLATE("BSplines: 1d BSpline evaluation", data_t, float, double)
      17             : // TEST_CASE_TEMPLATE("BSplines: 1d BSpline evaluation", data_t, float)
      18             : // {
      19             : //     constexpr auto size = 11;
      20             : //
      21             : //     const auto low = -2;
      22             : //     const auto high = 2;
      23             : //
      24             : //     const auto linspace = Vector_t<data_t>::LinSpaced(size, low, high);
      25             : //     std::array<data_t, size> spline;
      26             : //
      27             : //     auto degree = 3;
      28             : //     for (int i = 0; i < size; ++i) {
      29             : //         for (int j = 0; j < size; ++j) {
      30             : //             RealVector_t coord({{linspace[i], linspace[j]}});
      31             : //             // fmt::print("bspline({}, {}) = {}\n", coord.format(vecfmt), degree,
      32             : //             //            bspline::nd_bspline_evaluate(coord, degree));
      33             : //         }
      34             : //     }
      35             : //
      36             : //     // fmt::print("{}\n", linspace.format(vecfmt));
      37             : //     // fmt::print("{}\n", spline);
      38             : // }
      39             : 
      40             : /// Quick and easy simpsons rule for numerical integration
      41             : template <typename data_t, typename Func>
      42             : data_t simpsons(Func f, data_t a, data_t b, int N = 50)
      43             : {
      44             :     const auto h = (b - a) / N;
      45             : 
      46             :     data_t sum_odds = 0.0;
      47             :     for (int i = 1; i < N; i += 2) {
      48             :         sum_odds += f(a + i * h);
      49             :     }
      50             : 
      51             :     data_t sum_evens = 0.0;
      52             :     for (int i = 2; i < N; i += 2) {
      53             :         sum_evens += f(a + i * h);
      54             :     }
      55             : 
      56             :     return (f(a) + f(b) + 2 * sum_evens + 4 * sum_odds) * h / 3;
      57             : }
      58             : 
      59             : /// Integrate from distance s from center, to a
      60             : /// see Fig 3.6 (p. 53) and Listing 3.1 (p. 58) in Three-Dimensional Digital Tomosynthesis by
      61             : /// Levakhina
      62             : template <typename data_t>
      63             : data_t bspline_line_integral(data_t s, index_t m, index_t dim)
      64             : {
      65             :     auto x = [=](auto t) {
      66             :         data_t res = bspline::bspline1d_evaluate<data_t>(std::sqrt(t * t + s * s), m);
      67             :         for (int i = 1; i < dim; ++i) {
      68             :             res *= bspline::bspline1d_evaluate<data_t>(0., m);
      69             :         }
      70             :         return res;
      71             :     };
      72             : 
      73             :     return 2 * simpsons<data_t>(x, 0, 3);
      74             : }
      75             : 
      76             : TEST_CASE_TEMPLATE("BSplines: 1d line integal", data_t, float)
      77           1 : {
      78           1 :     constexpr auto size = 21;
      79             : 
      80           1 :     const auto low = -2;
      81           1 :     const auto high = 2;
      82             : 
      83           1 :     const auto linspace = Vector_t<data_t>::LinSpaced(size, low, high);
      84             : 
      85           1 :     const int degree = 2;
      86             : 
      87             :     // const data_t distance = 1.f;
      88             :     //
      89             :     // fmt::print("1D bspline at distance {:4.2f}: {:8.5f}\n", distance,
      90             :     //            bspline::bspline1d_evaluate<data_t>(distance, degree));
      91             :     // fmt::print("2D bspline at distance {:4.2f}: {:8.5f}\n", distance,
      92             :     //            bspline::bspline1d_evaluate<data_t>(0, degree)
      93             :     //                * bspline::bspline1d_evaluate<data_t>(distance, degree));
      94             :     // fmt::print("2D line integral: {:8.5f}\n", bspline_line_integral<data_t>(distance, degree,
      95             :     // 2)); fmt::print("3D line integral: {:8.5f}\n", bspline_line_integral<data_t>(distance,
      96             :     // degree, 3)); MESSAGE("helloo");
      97             : 
      98             :     // BSpline<data_t> bspline_1d(1, degree);
      99             :     // BSpline<data_t> bspline_2d(2, degree);
     100             :     //
     101             :     // for (int i = 0; i < 101; ++i) {
     102             :     //     const data_t x = (i / 25.) - 2.;
     103             :     //
     104             :     //     CAPTURE(x);
     105             :     //     CAPTURE(bspline::bspline1d_evaluate(x, degree));
     106             :     //     CAPTURE(bspline::bspline1d_evaluate(0., degree));
     107             :     //
     108             :     //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree),
     109             :     //              doctest::Approx(bspline_1d(Vector_t<data_t>({{x}}))));
     110             :     //
     111             :     //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(0.,
     112             :     //     degree),
     113             :     //              doctest::Approx(bspline_2d(Vector_t<data_t>({{x, 0}}))));
     114             :     //     CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(x, degree),
     115             :     //              doctest::Approx(bspline_2d(Vector_t<data_t>({{x, x}}))));
     116             :     // }
     117             : 
     118           1 :     ProjectedBSpline<data_t> proj_bspline_2d(2, degree);
     119           1 :     ProjectedBSpline<data_t> proj_bspline_3d(3, degree);
     120             : 
     121           1 :     CHECK_EQ(proj_bspline_2d.order(), degree);
     122           1 :     CHECK_EQ(proj_bspline_3d.order(), degree);
     123           1 :     CHECK_EQ(proj_bspline_2d.dim(), 2);
     124           1 :     CHECK_EQ(proj_bspline_3d.dim(), 3);
     125             : 
     126          22 :     for (int i = 0; i < 21; ++i) {
     127          21 :         const data_t x = (i / 5.) - 2.;
     128             : 
     129          21 :         CAPTURE(x);
     130          21 :         CAPTURE(bspline::bspline1d_evaluate(x, degree));
     131          21 :         CAPTURE(bspline::bspline1d_evaluate(0., degree));
     132             : 
     133          21 :         CHECK_EQ(bspline::bspline1d_evaluate(x, degree), doctest::Approx(proj_bspline_2d(x)));
     134             : 
     135          21 :         CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(0., degree),
     136          21 :                  doctest::Approx(proj_bspline_3d(x)));
     137          21 :     }
     138           1 : }