test_BSplines.cpp 4.54 KB
Newer Older
David Frank's avatar
David Frank committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
#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);
    CHECK_EQ(proj_bspline_3d.dim(), 3);

    for (int i = 0; i < 21; ++i) {
        const data_t x = (i / 5.) - 2.;

        CAPTURE(x);
        CAPTURE(bspline::bspline1d_evaluate(x, degree));
        CAPTURE(bspline::bspline1d_evaluate(0., degree));

        CHECK_EQ(bspline::bspline1d_evaluate(x, degree), doctest::Approx(proj_bspline_2d(x)));

        CHECK_EQ(bspline::bspline1d_evaluate(x, degree) * bspline::bspline1d_evaluate(0., degree),
                 doctest::Approx(proj_bspline_3d(x)));
    }
}