Line data Source code
1 : #pragma once
2 :
3 : #include "elsaDefines.h"
4 : #include "Assertions.h"
5 : #include "Math.hpp"
6 : #include "Luts.hpp"
7 :
8 : #include <cmath>
9 :
10 : namespace elsa
11 : {
12 : namespace bspline
13 : {
14 : /// @brief Evaluate the 1-dimensional B-Spline of degree m, with m + 2 equally spaced knots
15 : ///
16 : /// This is based on the implementation given as in e.g.:
17 : /// Fast B-spline transforms for continuous image representation and interpolation, Unser
18 : /// et. al. (1991), equation 2.2
19 : ///
20 : /// @param x coordinate to evaluate 1-dimensional B-Spline at
21 : /// @param m order of B-Spline
22 : template <typename data_t>
23 : constexpr data_t bspline1d_evaluate(data_t x, index_t m) noexcept
24 993052 : {
25 993052 : const auto xs = x + static_cast<data_t>(m + 1) / data_t{2};
26 993052 : data_t res = 0;
27 :
28 4329466 : for (int n = 0; n <= m + 1; ++n) {
29 3336414 : data_t nf = static_cast<data_t>(n);
30 3336414 : const auto tmp1 = math::heaviside<data_t>(xs - nf, 0);
31 3336414 : const auto tmp2 = std::pow<data_t>(xs - nf, m);
32 3336414 : const auto tmp3 = static_cast<data_t>(math::binom(m + 1, n));
33 3336414 : const auto tmp4 = std::pow<data_t>(-1, n) / static_cast<data_t>(math::factorial(m));
34 3336414 : res += static_cast<data_t>(tmp1 * tmp2 * tmp3 * tmp4);
35 3336414 : }
36 :
37 993052 : return res;
38 993052 : }
39 :
40 : /// @brief Evaluate the 1-dimensional B-Spline Derivative of degree m
41 : ///
42 : /// @param x coordinate to evaluate 1-dimensional B-Spline at
43 : /// @param m order of B-Spline
44 : template <typename data_t>
45 : constexpr data_t bsplineDerivative1d_evaluate(data_t x, index_t m) noexcept
46 328210 : {
47 328210 : return bspline1d_evaluate(x + 0.5, m - 1) - bspline1d_evaluate(x - 0.5, m - 1);
48 328210 : }
49 :
50 : /// @brief Evaluate n-dimensional B-Spline of degree m. As B-Splines are separable, this is
51 : /// just the product of 1-dimensional B-Splines.
52 : ///
53 : /// @param x n-dimensional coordinate to evaluate B-Spline at
54 : /// @param m order of B-Spline
55 : template <typename data_t>
56 : constexpr data_t nd_bspline_evaluate(const Vector_t<data_t>& x, index_t m) noexcept
57 : {
58 : data_t res = bspline1d_evaluate(x[0], m);
59 : for (int i = 1; i < x.size(); ++i) {
60 : res *= bspline1d_evaluate(x[i], m);
61 : }
62 : return res;
63 : }
64 :
65 : /// @brief Evaluate n-dimensional B-Spline at a given coordinate for the first dimension,
66 : /// and at the center (i.e. `0`) at all the other dimensions. This is particular useful as
67 : /// an approximation during the calculation of the line integral
68 : ///
69 : /// @param x 1-dimensional coordinate to evaluate B-Spline at
70 : /// @param m order of B-Spline
71 : /// @param dim dimension of B-Spline
72 : template <typename data_t>
73 : constexpr data_t nd_bspline_centered(data_t x, index_t m, index_t dim) noexcept
74 164648 : {
75 164648 : data_t res = bspline1d_evaluate<data_t>(x, m);
76 222369 : for (int i = 1; i < dim; ++i) {
77 57721 : const auto inc = bspline1d_evaluate<data_t>(0.f, m);
78 57721 : res *= inc;
79 57721 : }
80 164648 : return res;
81 164648 : }
82 :
83 : /// @brief Evaluate the derivative of an n-dimensional B-Spline at a given coordinate for
84 : /// the first dimension, and at the center (i.e. `0`) at all the other dimensions. This is
85 : /// particular useful as an approximation during the calculation of the line integral
86 : ///
87 : /// @param x 1-dimensional coordinate to evaluate B-Spline at
88 : /// @param m order of B-Spline
89 : /// @param dim dimension of B-Spline
90 : template <typename data_t>
91 : constexpr data_t nd_bspline_derivative_centered(data_t x, int m, int dim) noexcept
92 328210 : {
93 328210 : data_t res = bsplineDerivative1d_evaluate(x, m);
94 442410 : for (int i = 1; i < dim; ++i) {
95 114200 : const auto inc = bspline1d_evaluate(0., m);
96 114200 : res *= inc;
97 114200 : }
98 328210 : return res;
99 328210 : }
100 :
101 : static constexpr index_t DEFAULT_ORDER = 2;
102 : } // namespace bspline
103 :
104 : /// @brief Represent a B-Spline basis function of a given dimension and order
105 : template <typename data_t>
106 : class BSpline
107 : {
108 : public:
109 : BSpline(index_t dim, index_t order);
110 :
111 : data_t operator()(Vector_t<data_t> x);
112 :
113 : data_t derivative(data_t s);
114 :
115 : index_t order() const;
116 :
117 : data_t radius() const;
118 :
119 : index_t dim() const;
120 :
121 : private:
122 : /// Dimension of B-Spline
123 : index_t dim_;
124 :
125 : /// Order of B-Spline
126 : index_t order_;
127 : };
128 :
129 : /// @brief Represent a projected B-Spline basis function of a given dimension and order.
130 : /// Projected B-Splines are again B-Spline of n-1 dimensions. Using the fact, that B-Splines
131 : /// are close to symmetrical, we can approximate the projection only based on distance.
132 : template <typename data_t, size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
133 : class ProjectedBSpline
134 : {
135 : public:
136 : constexpr ProjectedBSpline(index_t dim, index_t order)
137 : : dim_(dim),
138 : order_(order),
139 : radius_((order + 1) * 0.5),
140 164404 : lut_([this](data_t s) { return this->operator()<true>(s); }, radius_),
141 : derivative_lut_(
142 164404 : [this](data_t s) { return order_ > 0 ? this->derivative<true>(s) : 0; }, radius_),
143 : normalized_gradient_lut_(
144 164404 : [this](data_t s) { return order_ > 0 ? this->normalized_gradient(s) : 0; },
145 : radius_)
146 1668 : {
147 1668 : }
148 :
149 : template <bool accurate = false>
150 : constexpr data_t operator()(data_t x) const
151 171478 : {
152 171478 : if constexpr (accurate)
153 164648 : return bspline::nd_bspline_centered(x, order_, dim_ - 1);
154 6830 : else
155 6830 : return lut_(std::abs(x));
156 171478 : }
157 :
158 : template <bool accurate = false>
159 : constexpr data_t derivative(data_t x) const
160 164206 : {
161 164206 : if constexpr (accurate)
162 164206 : return bspline::nd_bspline_derivative_centered(x, order_, dim_ - 1);
163 164206 : else
164 164206 : return derivative_lut_(std::abs(x)) * math::sgn(x);
165 164206 : }
166 :
167 : constexpr data_t normalized_gradient(data_t x) const
168 164004 : {
169 : // compute f'(x)/x
170 164004 : if (x == 0)
171 1660 : x = 1e-10;
172 164004 : return bspline::nd_bspline_derivative_centered(x, order_, dim_ - 1) / x;
173 164004 : }
174 :
175 2 : constexpr index_t dim() const { return dim_; }
176 :
177 2 : constexpr index_t order() const { return order_; }
178 :
179 0 : constexpr data_t radius() const { return radius_; }
180 :
181 1662 : constexpr const Lut<data_t, N>& get_lut() const { return lut_; }
182 :
183 2 : constexpr const Lut<data_t, N>& get_derivative_lut() const { return derivative_lut_; }
184 :
185 : constexpr const Lut<data_t, N>& get_normalized_gradient_lut() const
186 0 : {
187 0 : return normalized_gradient_lut_;
188 0 : }
189 :
190 : private:
191 : /// Dimension of B-Spline
192 : index_t dim_;
193 :
194 : /// Order of B-Spline
195 : index_t order_;
196 :
197 : /// Radius of B-Spline
198 : data_t radius_;
199 :
200 : /// LUT for projected B-Spline
201 : Lut<data_t, N> lut_;
202 :
203 : /// LUT for projected derivative
204 : Lut<data_t, N> derivative_lut_;
205 :
206 : /// LUT for projected normalized gradient
207 : Lut<data_t, N> normalized_gradient_lut_;
208 : };
209 : } // namespace elsa
|
