Line data Source code
1 : #pragma once
2 :
3 : #include "elsaDefines.h"
4 : #include "TypeCasts.hpp"
5 :
6 : namespace elsa
7 : {
8 :
9 : namespace math
10 : {
11 :
12 : template <typename data_t>
13 : constexpr inline data_t sq(data_t x) noexcept
14 68267 : {
15 68267 : return x * x;
16 68267 : }
17 :
18 : /// Compute factorial \f$n!\f$ recursively
19 : constexpr inline index_t factorial(index_t n) noexcept
20 4834305 : {
21 4834305 : return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
22 4834305 : }
23 :
24 : /// Compute binomial coefficient
25 : constexpr inline index_t binom(index_t n, index_t k) noexcept
26 3357473 : {
27 3357473 : return (k > n) ? 0
28 3357473 : : (k == 0 || k == n) ? 1
29 3357472 : : (k == 1 || k == n - 1) ? n
30 1371366 : : (k + k < n) ? (binom(n - 1, k - 1) * n) / k
31 21038 : : (binom(n - 1, k) * n) / (n - k);
32 3357473 : }
33 :
34 : /**
35 : * Compute Heaviside-function
36 : * \f[
37 : * x \mapsto
38 : * \begin{cases}
39 : * 0: & x < 0 \\
40 : * c: & x = 0 \\
41 : * 1: & x > 0
42 : * \end{cases}
43 : * \f]
44 : */
45 : template <typename data_t>
46 : constexpr data_t heaviside(data_t x1, data_t c)
47 3336614 : {
48 3336614 : if (x1 == 0) {
49 18576 : return c;
50 3318038 : } else if (x1 < 0) {
51 1106574 : return 0;
52 2211464 : } else {
53 2211464 : return 1;
54 2211464 : }
55 3336614 : }
56 :
57 : template <typename T>
58 : constexpr inline int sgn(T val)
59 4834468 : {
60 4834468 : return (T(0) < val) - (val < T(0));
61 4834468 : }
62 :
63 : template <typename data_t>
64 : data_t lerp(data_t a, SelfType_t<data_t> b, SelfType_t<data_t> t)
65 : {
66 : if ((a <= 0 && b >= 0) || (a >= 0 && b <= 0))
67 : return t * b + (1 - t) * a;
68 :
69 : if (t == 1)
70 : return b;
71 :
72 : const data_t x = a + t * (b - a);
73 :
74 : if ((t > 1) == (b > a))
75 : return b < x ? x : b;
76 : else
77 : return x < b ? x : b;
78 : }
79 : } // namespace math
80 :
81 : /// proposed in Y. Meyer, Oscillating Patterns in Image Processing and Nonlinear Evolution
82 : /// Equations. AMS, 2001
83 : template <typename data_t>
84 : data_t meyerFunction(data_t x)
85 219570 : {
86 219570 : if (x < 0.f) {
87 116764 : return 0;
88 116764 : } else if (0.f <= x && x <= 1.f) {
89 102806 : return 35 * std::pow(x, 4.f) - 84 * std::pow(x, 5.f) + 70 * std::pow(x, 6.f)
90 102806 : - 20 * std::pow(x, 7.f);
91 102806 : } else {
92 0 : return 1;
93 0 : }
94 219570 : }
95 :
96 : namespace shearlet
97 : {
98 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
99 : /// tutorial, 2014
100 : template <typename data_t>
101 : data_t b(data_t w)
102 312936 : {
103 312936 : if (1 <= std::abs(w) && std::abs(w) <= 2) {
104 41196 : return std::sin(pi<data_t> / 2.f * meyerFunction(std::abs(w) - 1));
105 271740 : } else if (2 < std::abs(w) && std::abs(w) <= 4) {
106 21906 : return std::cos(pi<data_t> / 2.f * meyerFunction(1.f / 2.f * std::abs(w) - 1.f));
107 249834 : } else {
108 249834 : return 0;
109 249834 : }
110 312936 : }
111 :
112 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
113 : /// tutorial, 2014
114 : template <typename data_t>
115 : data_t phi(data_t w)
116 4608 : {
117 4608 : if (std::abs(w) <= 1.f / 2) {
118 6 : return 1;
119 4602 : } else if (1.f / 2 < std::abs(w) && std::abs(w) < 1) {
120 0 : return std::cos(pi<data_t> / 2.f * meyerFunction(2.f * std::abs(w) - 1));
121 4602 : } else {
122 4602 : return 0;
123 4602 : }
124 4608 : }
125 :
126 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
127 : /// tutorial, 2014
128 : template <typename data_t>
129 : data_t phiHat(data_t w, data_t h)
130 4608 : {
131 4608 : if (std::abs(h) <= std::abs(w)) {
132 2458 : return phi(w);
133 2458 : } else {
134 2150 : return phi(h);
135 2150 : }
136 4608 : }
137 :
138 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
139 : /// tutorial, 2014
140 : template <typename data_t>
141 : data_t psiHat1(data_t w)
142 156468 : {
143 156468 : return std::sqrt(std::pow(b(2.f * w), 2.f) + std::pow(b(w), 2.f));
144 156468 : }
145 :
146 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
147 : /// tutorial, 2014
148 : template <typename data_t>
149 : data_t psiHat2(data_t w)
150 156468 : {
151 156468 : if (w <= 0) {
152 80192 : return std::sqrt(meyerFunction(1 + w));
153 80192 : } else {
154 76276 : return std::sqrt(meyerFunction(1 - w));
155 76276 : }
156 156468 : }
157 :
158 : /// defined in Sören Häuser and Gabriele Steidl, Fast Finite Shearlet Transform: a
159 : /// tutorial, 2014
160 : template <typename data_t>
161 : data_t psiHat(data_t w, data_t h)
162 156672 : {
163 156672 : if (w == 0) {
164 204 : return 0;
165 156468 : } else {
166 156468 : return psiHat1(w) * psiHat2(h / w);
167 156468 : }
168 156672 : }
169 : } // namespace shearlet
170 :
171 : namespace axdt
172 : {
173 :
174 : template <typename data_t>
175 : inline data_t legendreAtZero(size_t l)
176 : {
177 : if (l == 0) {
178 : return 1;
179 : }
180 : if (l % 2 == 1) {
181 : return 0;
182 : }
183 :
184 : data_t ratio = ((l / 2) % 2 == 0) ? 1 : -1;
185 :
186 : while (l > 0) {
187 : ratio *= l - 1;
188 : ratio /= l;
189 : l -= 2;
190 : }
191 : return ratio;
192 : }
193 :
194 : template <typename data_t>
195 : data_t ratio_of_factorials(index_t x, index_t y)
196 19000 : {
197 19000 : data_t ratio{1};
198 19000 : if (x == y) {
199 0 : return ratio;
200 19000 : } else if (x > y) {
201 35968 : for (index_t i = y + 1; i <= x; ++i) {
202 26151 : ratio *= as<data_t>(i);
203 26151 : }
204 9817 : return ratio;
205 9817 : } else {
206 9183 : return data_t{1} / ratio_of_factorials<data_t>(y, x);
207 9183 : }
208 19000 : }
209 :
210 : template <typename data_t>
211 : data_t double_factorial(index_t x)
212 2554 : {
213 2554 : data_t df = 1.0;
214 5834 : while (x > 1) {
215 3280 : df *= as<data_t>(x);
216 3280 : x -= 2;
217 3280 : }
218 2554 : return df;
219 2554 : }
220 :
221 : template <typename data_t>
222 : Eigen::Matrix<data_t, Eigen::Dynamic, Eigen::Dynamic> assoc_legendre_pol(index_t l,
223 : data_t x)
224 677 : {
225 691 : if (x < -1 || x > 1) {
226 0 : throw std::invalid_argument("math::axdt::assoc_legendre_pol: Can only evaluate "
227 0 : "polynomials at x in the interval [-1,1].");
228 0 : }
229 :
230 677 : Eigen::Matrix<data_t, Eigen::Dynamic, Eigen::Dynamic> result =
231 677 : Eigen::Matrix<data_t, Eigen::Dynamic, Eigen::Dynamic>::Zero(l + 1, l + 1);
232 :
233 677 : result(0, 0) = as<data_t>(1);
234 :
235 2971 : for (int i = 1; i <= l; ++i) {
236 5124 : for (int j = 0; j <= i - 2; ++j) {
237 2830 : result(i, j) = ((2 * as<data_t>(i) - 1) * x * result(i - 1, j)
238 2830 : - as<data_t>(i + j - 1) * result(i - 2, j))
239 2830 : / as<data_t>(i - j);
240 2830 : }
241 2294 : result(i, i - 1) = x * as<data_t>(2 * i - 1) * result(i - 1, i - 1);
242 2294 : result(i, i) = as<data_t>(pow(-1, i) * double_factorial<data_t>(2 * i - 1)
243 2294 : * pow(1 - x * x, i / 2.0));
244 2294 : }
245 :
246 677 : return result;
247 677 : }
248 :
249 : template <typename data_t>
250 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> SH_basis_real(index_t l, data_t theta, data_t phi)
251 671 : {
252 :
253 671 : Eigen::Matrix<data_t, Eigen::Dynamic, 1> result((l + 1) * (l + 1));
254 :
255 671 : Eigen::Matrix<data_t, Eigen::Dynamic, Eigen::Dynamic> legendre_pol =
256 671 : assoc_legendre_pol(l, as<data_t>(std::cos(theta)));
257 :
258 4110 : for (int i = 0; i <= l; ++i) {
259 3439 : auto c = as<data_t>(sqrt((2 * i + 1) / (4.0 * pi_t)));
260 15090 : for (int j = -i; j <= i; ++j) {
261 11651 : if (j < 0) {
262 5299 : result(i * i + j + i) =
263 5299 : as<data_t>(sqrt(2) * pow(-1, j) * c
264 5299 : * sqrt(ratio_of_factorials<data_t>(i + j, i - j))
265 5299 : * legendre_pol(i, -j) * sin(as<data_t>(-j) * phi));
266 6352 : } else if (j == 0) {
267 3460 : result(i * i + j + i) = c * legendre_pol(i, j);
268 3460 : } else {
269 2892 : result(i * i + j + i) =
270 2892 : as<data_t>(sqrt(2) * pow(-1, j) * c
271 2892 : * sqrt(ratio_of_factorials<data_t>(i - j, i + j))
272 2892 : * legendre_pol(i, j) * cos(as<data_t>(j) * phi));
273 2892 : }
274 11651 : }
275 3439 : }
276 :
277 671 : return result;
278 671 : }
279 : } // namespace axdt
280 :
281 : /// @brief Compute the sign of the given value. Will return -1, for negative values, 1 for
282 : /// positive ones and 0 otherwise
283 : template <typename T, typename Ret = int>
284 : Ret sign(T val)
285 : {
286 : return static_cast<Ret>((T{0} < val) - (val < T{0}));
287 : }
288 :
289 : } // namespace elsa
|
