STIR 6.4.0
BSplines_weights.inl
Go to the documentation of this file.
1//
2//
3/*
4 Copyright (C) 2005 - 2009-10-27, Hammersmith Imanet Ltd
5 Copyright (C) 2011-07-01 - 2011, Kris Thielemans
6 This file is part of STIR.
7
8 SPDX-License-Identifier: Apache-2.0
9
10 See STIR/LICENSE.txt for details
11*/
21
22#include "stir/error.h"
23
24START_NAMESPACE_STIR
25namespace BSpline
26{
27
28template <class posT>
29class PieceWiseFunction
30{
31public:
32 typedef posT result_type;
33
34 virtual ~PieceWiseFunction() {}
35
36 virtual posT kernel_length_left() const = 0;
37 virtual int kernel_total_length() const = 0;
38 posT kernel_length_right() const { return -this->kernel_length_left() + this->kernel_total_length(); }
39
40 virtual posT function_piece(const posT x, int p) const = 0;
41 virtual posT derivative_piece(const posT x, int p) const = 0;
42 virtual int find_piece(const posT x) const = 0;
43 virtual int find_highest_piece() const = 0;
44 posT function(const posT x) const { return this->function_piece(x, find_piece(x)); }
45 posT derivative(const posT x) const { return this->derivative_piece(x, find_piece(x)); }
46 posT operator()(const posT x) const { return this->function(x); }
47 posT operator()(const posT p, const BSplineType) { return this->function(p); }
48};
49
50template <BSplineType bspline_type, class posT>
51class BSplineFunction : public PieceWiseFunction<posT>
52{
53};
54
55template <class posT>
56class BSplineFunction<near_n, posT> : public PieceWiseFunction<posT>
57{
58public:
59 BSplineFunction() {}
60 posT kernel_length_left() const override { return static_cast<posT>(.5); }
61 BOOST_STATIC_CONSTANT(int, kernel_total_length_value = 1);
62 int kernel_total_length() const override { return kernel_total_length_value; }
63 posT function_piece(const posT x, int p) const override
64 {
65 switch (p)
66 {
67 case 0:
68 return static_cast<posT>(1);
69 // todo probably should add another piece containing a just .5, now returns 0
70 default:
71 return 0;
72 }
73 }
74 posT derivative_piece(const posT, int) const override { return 0; }
75
76 int find_abs_piece(const posT absx) const { return static_cast<int>(absx + .5); }
77 int find_piece(const posT x) const override
78 {
79 const int abs_p = this->find_abs_piece(fabs(x));
80 if (x > 0)
81 return abs_p;
82 else
83 return -abs_p;
84 }
85 int find_highest_piece() const override { return 0; }
86};
87
88template <class posT>
89class BSplineFunction<linear, posT> : public PieceWiseFunction<posT>
90{
91public:
92 BSplineFunction() {}
93 posT kernel_length_left() const override { return static_cast<posT>(1); }
94 BOOST_STATIC_CONSTANT(int, kernel_total_length_value = 2);
95 int kernel_total_length() const override { return kernel_total_length_value; }
96 posT function_piece(const posT x, int p) const override
97 {
98 switch (p)
99 {
100 case 0:
101 return 1 - x;
102 case -1:
103 return 1 + x;
104 default:
105 return 0;
106 }
107 }
108 posT derivative_piece(const posT x, int p) const override
109 {
110 switch (p)
111 {
112 case 0:
113 return -1;
114 case -1:
115 return 1;
116 default:
117 return 0;
118 }
119 }
120 int find_piece(const posT x) const override { return static_cast<int>(floor(x)); }
121
122 int find_highest_piece() const override { return 0; }
123};
124
125template <class posT>
126class BSplineFunction<quadratic, posT> : public PieceWiseFunction<posT>
127{
128public:
129 BSplineFunction() {}
130 posT kernel_length_left() const override { return static_cast<posT>(1.5); }
131 BOOST_STATIC_CONSTANT(int, kernel_total_length_value = 3);
132 int kernel_total_length() const override { return kernel_total_length_value; }
133 posT function_piece(const posT x, int p) const override
134 {
135 switch (std::abs(p))
136 {
137 case 0:
138 return static_cast<posT>(3) / 4 - square(x);
139 case 1: {
140 const posT absx = std::fabs(x);
141 return square(2 * absx - 3) / 8;
142 }
143 default:
144 return 0;
145 }
146 }
147 posT derivative_piece(const posT x, int p) const override
148 {
149 switch (std::abs(p))
150 {
151 case 0:
152 return -2 * x;
153 case 1: {
154 const int sign = x > 0 ? 1 : -1;
155 return x - sign * static_cast<posT>(1.5);
156 }
157 default:
158 return 0;
159 }
160 }
161
162private:
163 int find_abs_piece(const posT absx) const
164 {
165#if 1
166 return static_cast<int>(absx + .5);
167#else
168 // can use this if guaranteed never out of range
169 if (absx <= .5)
170 return 0;
171 else if (absx <= 1.5)
172 return 1;
173 else
174 return 2;
175#endif
176 }
177
178public:
179 int find_piece(const posT x) const override
180 {
181 const int abs_p = this->find_abs_piece(fabs(x));
182 if (x > 0)
183 return abs_p;
184 else
185 return -abs_p;
186 }
187 int find_highest_piece() const override
188 {
189 return 1;
190 }
191
192#if 0
193 posT function(const posT x) const
194 {
195 const posT absx = fabs(x);
196 // note: gcc inlines next function call, even getting rid of switch
197 // so, this is faster than using virtual functions etc.
198 // however, ugly code repetition, so we can't be bothered, as this is hardly called anyway
199 if (absx<=.5)
200 return self_t::function_piece(absx,0);
201 else if (absx<=1.5)
202 return self_t::function_piece(absx,1);
203 else
204 return 0;
205 }
206#endif
207};
208
209template <class posT>
210class BSplineFunction<cubic, posT> : public PieceWiseFunction<posT>
211{
212public:
213 BSplineFunction() {}
214 posT kernel_length_left() const override { return static_cast<posT>(2); }
215 BOOST_STATIC_CONSTANT(int, kernel_total_length_value = 4);
216 int kernel_total_length() const override { return kernel_total_length_value; }
217 posT function_piece(const posT x, int p) const override
218 {
219 const posT absx = std::fabs(x);
220 switch (p)
221 {
222 case 0:
223 case -1:
224 return 2. / 3. + (absx / 2 - 1) * absx * absx;
225 case 1:
226 case -2: {
227 const posT tmp = 2 - absx;
228 return tmp * tmp * tmp / 6;
229 }
230 default:
231 return 0;
232 }
233 }
234 posT derivative_piece(const posT x, int p) const override
235 {
236 switch (p)
237 {
238 case 0:
239 return x * (1.5 * x - 2);
240 case -1:
241 return -x * (1.5 * x + 2);
242 case 1:
243 return -square(x - 2) / 2;
244 case -2:
245 return square(x + 2) / 2;
246 default:
247 return 0;
248 }
249 }
250 int find_piece(const posT x) const override { return static_cast<int>(floor(x)); }
251 int find_highest_piece() const override { return 1; }
252};
253
254template <class posT>
255class BSplineFunction<oMoms, posT> : public PieceWiseFunction<posT>
256{
257public:
258 BSplineFunction() {}
259 posT kernel_length_left() const override { return static_cast<posT>(2); }
260 BOOST_STATIC_CONSTANT(int, kernel_total_length_value = 4);
261 int kernel_total_length() const override { return kernel_total_length_value; }
262 posT function_piece(const posT x, int p) const override
263 {
264 const posT absx = std::fabs(x);
265 switch (p)
266 {
267 case 0:
268 case -1:
269 return 13. / 21. + absx * (1. / 14. + absx * (absx / 2 - 1));
270 case 1:
271 case -2:
272 return 29. / 21. + absx * (-85. / 42. + absx * (1 - absx / 6));
273 default:
274 return 0;
275 }
276 }
277 posT derivative_piece(const posT x, int p) const override
278 {
279 const posT absx = std::fabs(x);
280 const int sign = x > 0 ? 1 : -1;
281 switch (p)
282 {
283 case 0:
284 case -1:
285 return sign * (1 / 14. + absx * (-2 + 3 * absx / 2));
286 case 1:
287 case -2:
288 return sign * (-85 / 42. + absx * (2 - absx / 2));
289 default:
290 return 0;
291 }
292 }
293 int find_piece(const posT x) const override { return static_cast<int>(floor(x)); }
294 int find_highest_piece() const override { return 1; }
295};
296
297static const BSplineFunction<near_n, pos_type> near_n_BSpline_function;
298static const BSplineFunction<linear, pos_type> linear_BSpline_function;
299static const BSplineFunction<quadratic, pos_type> quadratic_BSpline_function;
300static const BSplineFunction<cubic, pos_type> cubic_BSpline_function;
301static const BSplineFunction<oMoms, pos_type> oMoms_BSpline_function;
302
303inline const PieceWiseFunction<pos_type>&
304bspline_function(BSplineType type)
305{
306 switch (type)
307 {
308 case near_n:
309 return near_n_BSpline_function;
310 case linear:
311 return linear_BSpline_function;
312 case quadratic:
313 return quadratic_BSpline_function;
314 case cubic:
315 return cubic_BSpline_function;
316 case oMoms:
317 return oMoms_BSpline_function;
318 default:
319 std::cerr << "quartic,quantic to do\n";
320 exit(EXIT_FAILURE);
321 return cubic_BSpline_function; // WARNING WRONG
322 }
323}
324
325template <typename posT>
326inline posT
327cubic_BSplines_weight(const posT relative_position)
328{
329 const posT abs_relative_position = fabs(relative_position);
330 assert(abs_relative_position >= 0);
331 if (abs_relative_position < 1)
332 return 2. / 3. + (0.5 * abs_relative_position - 1) * abs_relative_position * abs_relative_position;
333 if (abs_relative_position >= 2)
334 return 0;
335 const posT tmp = 2 - abs_relative_position;
336 return tmp * tmp * tmp / 6;
337}
338
339template <typename posT>
340inline posT
341oMoms_weight(const posT relative_position)
342{
343 const posT abs_relative_position = fabs(relative_position);
344 assert(abs_relative_position >= 0);
345 if (abs_relative_position >= 2)
346 return 0;
347 if (abs_relative_position >= 1)
348 return 29. / 21. + abs_relative_position * (-85. / 42. + abs_relative_position * (1. - abs_relative_position / 6));
349 else
350 return 13. / 21. + abs_relative_position * (1. / 14. + abs_relative_position * (0.5 * abs_relative_position - 1.));
351}
352
353template <typename posT>
354inline posT
355cubic_BSplines_1st_der_weight(const posT relative_position)
356{
357 const posT abs_relative_position = fabs(relative_position);
358 if (abs_relative_position >= 2)
359 return 0;
360 int sign = relative_position > 0 ? 1 : -1;
361 if (abs_relative_position >= 1)
362 return -0.5 * sign * (abs_relative_position - 2) * (abs_relative_position - 2);
363 return sign * abs_relative_position * (1.5 * abs_relative_position - 2.);
364}
365
366template <typename posT>
367posT
368BSplines_1st_der_weight(const posT relative_position, const BSplineType spline_type)
369{
370 switch (spline_type)
371 {
372 case cubic:
373 return cubic_BSplines_1st_der_weight(relative_position);
374 default:
375 error("BSplines_1st_der_weight currently only implemented for cubic splines.");
376 return 0;
377 // TODO
378 }
379}
380
381template <typename posT>
382posT
383BSplines_weights(const posT relative_position, const BSplineType spline_type)
384{
385 // double relative_position = rel_position;
386 switch (spline_type)
387 {
388 case cubic:
389 return cubic_BSplines_weight(relative_position);
390 case near_n: {
391 if (fabs(relative_position) < 0.5)
392 return 1;
393 if (fabs(relative_position) == 0.5)
394 return 0.5;
395 else
396 return 0;
397 }
398 case linear: {
399 if (fabs(relative_position) < 1)
400 return 1 - fabs(relative_position);
401 else
402 return 0;
403 }
404 case quadratic:
405 return BSplineFunction<quadratic, posT>().function(relative_position);
406 case quintic:
407 return (pow(std::max(0., -3. + relative_position), 5) - 6 * pow(std::max(0., -2. + relative_position), 5)
408 + 15 * pow(std::max(0., -1. + relative_position), 5)
409 - 20 * pow(std::max(static_cast<posT>(0), relative_position), 5) + 15 * pow(std::max(0., 1. + relative_position), 5)
410 - 6 * pow(std::max(0., 2. + relative_position), 5) + pow(std::max(0., 3. + relative_position), 5))
411 / 120.;
412 case oMoms:
413 return oMoms_weight(relative_position);
414 default:
415 error("Not implemented b-spline type");
416 return -1000000;
417 }
418}
419
420} // namespace BSpline
421
422END_NAMESPACE_STIR
error.h
Declaration of stir::error()
stir::BSpline::BSplines_1st_der_weight
pos_type BSplines_1st_der_weight(const pos_type relative_position, const BSplineType spline_type)
return value of the first derivative of the spline
stir::BSpline::BSplines_weights
pos_type BSplines_weights(const pos_type relative_position, const BSplineType spline_type)
return spline value
stir::error
void error(const char *const s,...)
Print error with format string a la printf and throw exception.
Definition error.cxx:42
stir::square
NUMBER square(const NUMBER &x)
returns the square of a number, templated.
Definition common.h:154
Generated on Wed May 27 2026 17:04:05 for STIR by doxygen 1.13.2