Line data Source code
1 : #pragma once
2 :
3 : #include "Error.h"
4 : #include "elsaDefines.h"
5 : #include "Bessel.h"
6 : #include "Luts.hpp"
7 : #include "Math.hpp"
8 :
9 : namespace elsa
10 : {
11 : namespace blobs
12 : {
13 : // The devinition given by Lewitt (1992):
14 : // \f$
15 : // b_{m, alpha, a}(r) = \frac{I_m(\alpha * \sqrt{1 - (r/a)^2)}{I_m(\alpha)} * (\sqrt{1 -
16 : // (r/a)^2})^m \f$ for any \f$ 0 <= r <= a \f$.
17 : template <typename data_t>
18 : constexpr data_t blob_evaluate(data_t r, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
19 : index_t m) noexcept
20 9104328 : {
21 9104328 : const auto w = static_cast<data_t>(1) - std::pow(r / a, static_cast<data_t>(2));
22 9104328 : if (w >= 0) {
23 7233828 : const data_t Im1 = static_cast<data_t>(math::bessi(m, alpha));
24 7233828 : const data_t arg = std::sqrt(w);
25 7233828 : const data_t Im2 = static_cast<data_t>(math::bessi(m, alpha * arg));
26 :
27 7233828 : return (Im2 / Im1) * static_cast<data_t>(std::pow(arg, m));
28 7233828 : }
29 1870500 : return 0;
30 1870500 : }
31 :
32 : /// @brief Compute line integral of blob through a straight line
33 : /// The exact formulations are quite easily derived using the fact, that the line
34 : /// integral for a single blob is given by:
35 : ///
36 : /// @f[
37 : /// \frac{a}{I_m(\alpha)} \sqrt{\frac{2 \pi}{\alpha}} \left(\sqrt{1 -
38 : /// \left(\frac{s}{a}\right)^2}\right)^{m + 0.5} I_{m + 0.5}\left(\alpha \sqrt{1 -
39 : /// \left(\frac{s}{a}\right)^2}\right)
40 : /// @f]
41 : ///
42 : /// Then for a given order substitute I_{m + 0.5}(x) with the elementary function
43 : /// representations, as they can be found in Abramowitz and Stegun's Handbook of
44 : /// mathematical functions (1972):
45 : /// - \f$ I_{0.5}(x) = \sqrt{\frac{2}{\pi x}} \sinh(x)\$f
46 : /// - \f$ I_{1.5}(x) = \sqrt{\frac{2}{\pi x}} \left( \cosh(x) - \frac{\sinh(x)}{x} \right)
47 : /// \$f
48 : /// - \f$ I_{2.5}(x) = \sqrt{\frac{2}{\pi x}} \left(\left(\frac{3}{x^2} +
49 : /// \frac{1}{x}\right)\sinh(x) - \frac{3}{x^2} \cosh(x)\right)\$f
50 : ///
51 : /// Which will result in the below formulations. In theory using the recurrent relations,
52 : /// this could be extended to further orders, but is deemed unnecessary for now.
53 : ///
54 : /// TODO: What about alpha = 0? Check if you find anything in the literature for that.
55 : ///
56 : /// Ref:
57 : /// Distance-Driven Projection and Backprojection for Spherically Symmetric Basis Functions
58 : /// in CT - Levakhina
59 : /// Spherically symmetric volume elements as basis functions for image reconstructions in
60 : /// computed laminography - P. Trampert
61 : /// Semi-Discrete Iteration Methods in X-Ray Tomography - Jonas Vogelgesang
62 : ///
63 : /// @param distance distance of blob center to straight line, in literature often referred
64 : /// to as `r`
65 : /// @param radius radius of blob, often referred to as `a`
66 : /// @param alpha smoothness factor of blob, expected to be larger than 0
67 : /// @param order order of Bessel function, often referred to as `m`
68 : template <typename data_t>
69 : constexpr data_t blob_projected(data_t s, SelfType_t<data_t> a, SelfType_t<data_t> alpha,
70 : index_t m)
71 410936 : {
72 : // expect alpha > 0
73 410936 : using namespace elsa::math;
74 :
75 410936 : const data_t sda = s / a;
76 410936 : const data_t sdas = std::pow(sda, 2.f);
77 410936 : const data_t w = 1.f - sdas;
78 :
79 410936 : if (w > 1.0e-10) {
80 406826 : const auto arg = alpha * std::sqrt(w);
81 406826 : if (m == 0) {
82 2376 : return (2.f * a / alpha) * std::sinh(arg)
83 2376 : / static_cast<data_t>(math::bessi0(alpha));
84 404450 : } else if (m == 1) {
85 5560 : return (2.f * a / alpha) * std::sqrt(w)
86 5560 : * (std::cosh(arg) - std::sinh(arg) / arg)
87 5560 : / static_cast<data_t>(math::bessi1(alpha));
88 :
89 398890 : } else if (m == 2) {
90 398890 : return (2.f * a / alpha) * w
91 398890 : * ((3.f / (arg * arg) + 1.f) * std::sinh(arg)
92 398890 : - (3.f / arg) * std::cosh(arg))
93 398890 : / static_cast<data_t>(math::bessi2(alpha));
94 398890 : } else {
95 0 : throw Error("m out of range in blob_projected()");
96 0 : }
97 4110 : }
98 4110 : return 0.0f;
99 4110 : }
100 :
101 : template <typename data_t>
102 : constexpr data_t blob_projected(data_t s)
103 : {
104 : return blob_projected(s, 2.f, 10.83f, 2);
105 : }
106 :
107 : /// @brief Compute line integral of blob derivative through a straight line
108 : /// The exact formulations are quite easily derived using the fact, that the line
109 : /// integral for the derivative of a single blob is given by:
110 : ///
111 : /// @f[
112 : /// - \frac{\sqrt{2 \pi \alpha}}{I_m(\alpha)} \frac{s}{a} \left(\sqrt{1 -
113 : /// \left(\frac{s}{a}\right)^2}\right)^{m - 0.5} I_{m - 0.5}\left(\alpha \sqrt{1 -
114 : /// \left(\frac{s}{a}\right)^2}\right)
115 : /// @f]
116 : ///
117 : /// By substitution with the above formulas (blob_projected) one can derive the below
118 : /// formulations. In theory using the recurrent relations, this could be extended to
119 : /// further orders, but is deemed unnecessary for now.
120 : ///
121 : /// Ref:
122 : /// Investigation of discrete imaging models and iterative image reconstruction
123 : /// in differential X-ray phase-contrast tomography - Qiaofeng Xu (Appendix B)
124 : ///
125 : /// @param distance distance of blob center to straight line, in literature often referred
126 : /// to as `r`
127 : /// @param radius radius of blob, often referred to as `a`
128 : /// @param alpha smoothness factor of blob, expected to be larger than 0
129 : /// @param order order of Bessel function, often referred to as `m`
130 : template <typename data_t>
131 : constexpr data_t blob_derivative_projected(data_t s, SelfType_t<data_t> a,
132 : SelfType_t<data_t> alpha, int m)
133 411666 : {
134 : // expect alpha > 0
135 411666 : using namespace elsa::math;
136 :
137 411666 : const data_t sda = s / a;
138 411666 : const data_t sdas = std::pow(sda, 2);
139 411666 : const data_t w = 1.0 - sdas;
140 :
141 411666 : if (w > 1.0e-10) {
142 407550 : const auto arg = alpha * std::sqrt(w);
143 407550 : if (m == 1) {
144 7160 : return (-2.0 * s / a) * std::sinh(arg) / bessi1(alpha);
145 :
146 400390 : } else if (m == 2) {
147 400390 : return (-2.0 * s / alpha / a) * (std::cosh(arg) * arg - std::sinh(arg))
148 400390 : / bessi2(alpha);
149 400390 : } else {
150 0 : throw Error("m out of range in blob_derivative_projected()");
151 0 : }
152 4116 : }
153 4116 : return 0.0;
154 4116 : }
155 :
156 : /// @brief For 3D objects we need to compute the directional gradient given by
157 : /// @f$ \frac{-g(\left\lVert \vec x\right\rVert)}{\left\lVert \vec x\right\rVert} \vec x^T
158 : /// @f$
159 : /// where @f$g@f$ is the derivative computed below
160 : ///
161 : /// As we divide by the argument we have a potential divide by zero. This can be solved by
162 : /// moving the division inside g as it cancels out
163 : ///
164 : /// Compute line integral of blob derivative through a straight line
165 : /// The exact formulations are quite easily derived using the fact, that the line
166 : /// integral for the derivative of a single blob is given by:
167 : ///
168 : /// @f[
169 : /// - \frac{\sqrt{2 \pi \alpha}}{I_m(\alpha)} \frac{s}{a} \left(\sqrt{1 -
170 : /// \left(\frac{s}{a}\right)^2}\right)^{m - 0.5} I_{m - 0.5}\left(\alpha \sqrt{1 -
171 : /// \left(\frac{s}{a}\right)^2}\right)
172 : /// @f]
173 : ///
174 : /// By substitution with the above formulas (blob_projected) one can derive the below
175 : /// formulations. In theory using the recurrent relations, this could be extended to
176 : /// further orders, but is deemed unnecessary for now.
177 : ///
178 : /// Ref:
179 : /// Investigation of discrete imaging models and iterative image reconstruction
180 : /// in differential X-ray phase-contrast tomography - Qiaofeng Xu (Appendix B)
181 : ///
182 : /// @param distance distance of blob center to straight line, in literature often referred
183 : /// to as `r`
184 : /// @param radius radius of blob, often referred to as `a`
185 : /// @param alpha smoothness factor of blob, expected to be larger than 0
186 : /// @param order order of Bessel function, often referred to as `m`
187 : template <typename data_t>
188 : constexpr data_t blob_normalized_derivative_projected(data_t s, SelfType_t<data_t> a,
189 : SelfType_t<data_t> alpha, int m)
190 411634 : {
191 : // expect alpha > 0
192 411634 : using namespace elsa::math;
193 :
194 411634 : const data_t sda = s / a;
195 411634 : const data_t sdas = std::pow(sda, 2);
196 411634 : const data_t w = 1.0 - sdas;
197 :
198 411634 : if (w > 1.0e-10) {
199 407518 : const auto arg = alpha * std::sqrt(w);
200 407518 : if (m == 1) {
201 7144 : return (-2.0 / a) * std::sinh(arg) / bessi1(alpha);
202 :
203 400374 : } else if (m == 2) {
204 400374 : return (-2.0 / alpha / a) * (std::cosh(arg) * arg - std::sinh(arg))
205 400374 : / bessi2(alpha);
206 400374 : } else {
207 0 : throw Error("m out of range in blob_normalized_derivative_projected()");
208 0 : }
209 4116 : }
210 4116 : return 0.0;
211 4116 : }
212 :
213 : template <typename data_t>
214 : constexpr data_t blob_derivative_projected(data_t s)
215 : {
216 : return blob_derivative_projected(s, 2.f, 10.83f, 2);
217 : }
218 :
219 : static constexpr double DEFAULT_RADIUS = 2.0;
220 : static constexpr double DEFAULT_ALPHA = 10.83;
221 : static constexpr index_t DEFAULT_ORDER = 2;
222 :
223 : } // namespace blobs
224 :
225 : template <typename data_t>
226 : class Blob
227 : {
228 : public:
229 : constexpr Blob(data_t radius, SelfType_t<data_t> alpha, index_t order)
230 : : radius_(radius), alpha_(alpha), order_(order)
231 96 : {
232 96 : }
233 : constexpr data_t operator()(data_t s)
234 9696 : {
235 9696 : return blobs::blob_evaluate(s, radius_, alpha_, order_);
236 9696 : }
237 :
238 : constexpr data_t radius() const { return radius_; }
239 :
240 : constexpr data_t alpha() const { return alpha_; }
241 :
242 : constexpr index_t order() const { return order_; }
243 :
244 : private:
245 : data_t radius_;
246 : data_t alpha_;
247 : index_t order_;
248 : };
249 :
250 : template <typename data_t, size_t N = DEFAULT_PROJECTOR_LUT_SIZE>
251 : class ProjectedBlob
252 : {
253 : public:
254 : constexpr ProjectedBlob(data_t radius, SelfType_t<data_t> alpha, index_t order)
255 : : radius_(radius),
256 : alpha_(alpha),
257 : order_(order),
258 410834 : lut_([this](data_t s) { return this->operator()<true>(s); }, radius_),
259 : derivative_lut_(
260 410834 : [this](data_t s) { return order_ > 0 ? this->derivative<true>(s) : 0; }, radius_),
261 : normalized_gradient_lut_(
262 410834 : [this](data_t s) { return order_ > 0 ? this->normalized_gradient(s) : 0; },
263 : radius_)
264 4108 : {
265 4108 : }
266 :
267 : template <bool accurate = false>
268 : constexpr data_t operator()(data_t s) const
269 420235 : {
270 420235 : if constexpr (accurate)
271 410936 : return blobs::blob_projected(s, radius_, alpha_, order_);
272 9299 : else
273 9299 : return lut_(std::abs(s));
274 420235 : }
275 :
276 : template <bool accurate = false>
277 : constexpr data_t derivative(data_t s) const
278 411762 : {
279 411762 : if constexpr (accurate)
280 411666 : return blobs::blob_derivative_projected(s, radius_, alpha_, order_);
281 96 : else
282 96 : return derivative_lut_(std::abs(s)) * math::sgn(s);
283 411762 : }
284 :
285 : constexpr data_t normalized_gradient(data_t s) const
286 411634 : {
287 411634 : return blobs::blob_normalized_derivative_projected(s, radius_, alpha_, order_);
288 411634 : }
289 :
290 1594 : constexpr data_t radius() const { return radius_; }
291 :
292 0 : constexpr data_t alpha() const { return alpha_; }
293 :
294 0 : constexpr index_t order() const { return order_; }
295 :
296 4256083 : constexpr const Lut<data_t, N>& get_lut() const { return lut_; }
297 :
298 4203132 : constexpr const Lut<data_t, N>& get_derivative_lut() const { return derivative_lut_; }
299 :
300 : constexpr const Lut<data_t, N>& get_normalized_gradient_lut() const
301 0 : {
302 0 : return normalized_gradient_lut_;
303 0 : }
304 :
305 : private:
306 : data_t radius_;
307 : data_t alpha_;
308 : index_t order_;
309 :
310 : /// LUT for projected Blob
311 : Lut<data_t, N> lut_;
312 :
313 : /// LUT for projected derivative
314 : Lut<data_t, N> derivative_lut_;
315 :
316 : /// LUT for projected normalized gradient
317 : Lut<data_t, N> normalized_gradient_lut_;
318 : };
319 :
320 : } // namespace elsa
|
