STIR 6.4.0
BSplinesDetail.inl
Go to the documentation of this file.
1/*
2 Copyright (C) 2005- 2009, Hammersmith Imanet Ltd
3 Copyright (C) 2013, University College London
4 This file is part of STIR.
5
6 SPDX-License-Identifier: Apache-2.0
7
8 See STIR/LICENSE.txt for details
9*/
18#include "stir/assign.h"
19
20#ifndef __stir_numerics_BSplinesDetail_inl_
21# define __stir_numerics_BSplinesDetail_inl_
22START_NAMESPACE_STIR
23
24namespace BSpline
25{
27namespace detail
28{
29template <typename constantsT>
30static inline void
31set_BSpline_values(constantsT& z1, constantsT& z2, constantsT& lambda, const BSplineType spline_type)
32{
33 switch (spline_type)
34 {
35 case near_n:
36 z1 = static_cast<constantsT>(0);
37 z2 = static_cast<constantsT>(0);
38 break;
39 case linear:
40 z1 = static_cast<constantsT>(0);
41 z2 = static_cast<constantsT>(0);
42 break;
43 case quadratic:
44 z1 = static_cast<constantsT>(sqrt(8.) - 3);
45 z2 = static_cast<constantsT>(0);
46 break;
47 case cubic:
48 z1 = static_cast<constantsT>(sqrt(3.) - 2);
49 z2 = static_cast<constantsT>(0);
50 break;
51 case quartic:
52 z1 = static_cast<constantsT>(sqrt(664. - sqrt(438976.)) + sqrt(304.) - 19.);
53 z2 = static_cast<constantsT>(sqrt(664. - sqrt(438976.)) - sqrt(304.) - 19.);
54 break;
55 case quintic:
56 z1 = static_cast<constantsT>(0.5 * (sqrt(270. - sqrt(70980.)) + sqrt(105.) - 13.));
57 z2 = static_cast<constantsT>(0.5 * (sqrt(270. - sqrt(70980.)) - sqrt(105.) - 13.));
58 break;
59 case oMoms:
60 z1 = static_cast<constantsT>((sqrt(105.) - 13.) / 8.);
61 z2 = static_cast<constantsT>(0);
62 break;
63 }
64 lambda = static_cast<constantsT>((1. - z1) * (1. - (1. / z1)));
65 if (z2 != static_cast<constantsT>(0))
66 lambda *= static_cast<constantsT>((1. - z2) * (1. - (1. / z2)));
67}
68
69// 1d specialisation
70template <typename out_elemT, typename in_elemT, typename constantsT>
71void
72set_coef(Array<1, out_elemT>& coeffs,
73 const Array<1, in_elemT>& input,
74 const BasicCoordinate<1, constantsT>& z1s,
75 const BasicCoordinate<1, constantsT>& z2s,
76 const BasicCoordinate<1, constantsT>& lambdas)
77{
78 BSplines_coef(coeffs.begin(), coeffs.end(), input.begin(), input.end(), z1s[1], z2s[1], lambdas[1]);
79}
80
81template <int num_dimensions, typename out_elemT, typename in_elemT, typename constantsT>
82void
83set_coef(Array<num_dimensions, out_elemT>& coeffs,
84 const Array<num_dimensions, in_elemT>& input,
85 const BasicCoordinate<num_dimensions, constantsT>& z1s,
86 const BasicCoordinate<num_dimensions, constantsT>& z2s,
87 const BasicCoordinate<num_dimensions, constantsT>& lambdas)
88{
89 Array<num_dimensions, out_elemT> temp(input.get_index_range());
90 BSplines_coef(temp.begin(), temp.end(), input.begin(), input.end(), z1s[1], z2s[1], lambdas[1]);
91
92 for (int i = coeffs.get_min_index(); i <= coeffs.get_max_index(); ++i)
93 {
94 set_coef(coeffs[i], temp[i], cut_first_dimension(z1s), cut_first_dimension(z2s), cut_first_dimension(lambdas));
95 }
96}
97
98# if 0
99 template <typename pos_type>
100 struct BW
101 {
102 typedef pos_type result_type;
103 pos_type operator()(const pos_type p, const BSplineType type)
104 {
105 return BSplines_weights(p, type);
106 }
107 };
108
109 template <typename pos_type>
110 struct Bder
111 {
112 typedef pos_type result_type;
113 pos_type operator()(const pos_type p, const BSplineType type)
114 {
115 return BSplines_1st_der_weight(p, type);
116 }
117 };
118# else
119template <typename pos_type>
120struct BW
121{
122 typedef pos_type result_type;
123 pos_type operator()(const pos_type pos, int piece, const PieceWiseFunction<pos_type>& f)
124 {
125 return f.function_piece(pos, piece);
126 }
127};
128
129template <typename pos_type>
130struct Bder
131{
132 typedef pos_type result_type;
133 pos_type operator()(const pos_type pos, int piece, const PieceWiseFunction<pos_type>& f)
134 {
135 return f.derivative_piece(pos, piece);
136 }
137};
138# endif
139
140// TODO later
141/*
142 get_value could be used normally to say get_value(coeffs,index) == coeffs[index],
143 but would allows to use e.g. periodic boundary conditions or extrapolation
144*/
145template <int num_dimensions, int num_dimensions2, typename T, typename FunctionT, typename SplineFunctionT>
146inline typename SplineFunctionT::result_type
147spline_convolution(const Array<num_dimensions, T>& coeffs,
148 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
149 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types,
150 FunctionT f,
151 SplineFunctionT g)
152{
153 const int current_dimension = num_dimensions2 - num_dimensions + 1;
154 const PieceWiseFunction<pos_type>& bspline = bspline_function(spline_types[current_dimension]);
155
156 typename SplineFunctionT::result_type value;
157 assign(value, 0);
158 const int kmin = static_cast<int>(std::ceil(relative_positions[current_dimension] - bspline.kernel_length_right()));
159 const int kmax = kmin + bspline.kernel_total_length() - 1;
160 int k = kmin;
161 pos_type current_pos = relative_positions[current_dimension] - k;
162# define NNN
163# ifdef NNN
164 // TODO doesn't work yet when relative_position is an integer: kmin then becomes highest_piece+1
165 // int p=bspline.find_highest_piece();
166 // assert(p == bspline.find_piece(current_pos));
167 int p = bspline.find_piece(current_pos);
168# define DECR_P , --p
169# else
170# define DECR_P
171# endif
172 for (; k <= kmax; ++k, --current_pos DECR_P)
173 {
174 int index;
175 if (k < coeffs.get_min_index())
176 index = 2 * coeffs.get_min_index() - k;
177 else if (k > coeffs.get_max_index())
178 index = 2 * coeffs.get_max_index() - k;
179 else
180 index = k;
181 assert(coeffs.get_min_index() <= index && index <= coeffs.get_max_index());
182 value += static_cast<typename SplineFunctionT::result_type>(g(coeffs[index], relative_positions, spline_types) *
183# ifdef NNN
184 f(current_pos, p, bspline)
185# else
186 f(current_pos, spline_types[current_dimension])
187# endif
188 );
189 }
190 return value;
191}
192
193template <int num_dimensions2, typename T, typename FunctionT>
194inline T
195spline_convolution(const Array<1, T>& coeffs,
196 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
197 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types,
198 FunctionT f)
199{
200 const int current_dimension = num_dimensions2;
201 const PieceWiseFunction<pos_type>& bspline = bspline_function(spline_types[current_dimension]);
202 T value;
203 assign(value, 0);
204 // x-1.5<k<x+1.5
205 const int kmin = static_cast<int>(std::ceil(relative_positions[current_dimension] - bspline.kernel_length_right()));
206 const int kmax = kmin + bspline.kernel_total_length() - 1;
207 int k = kmin;
208 pos_type current_pos = relative_positions[current_dimension] - k;
209# ifdef NNN
210 // TODO doesn't work yet when relative_position is an integer: kmin then becomes highest_piece+1
211 // int p=bspline.find_highest_piece();
212 // assert(p == bspline.find_piece(current_pos));
213 int p = bspline.find_piece(current_pos);
214# define DECR_P , --p
215# else
216# define DECR_P
217# endif
218 for (; k <= kmax; ++k, --current_pos DECR_P)
219 {
220 int index;
221 if (k < coeffs.get_min_index())
222 index = 2 * coeffs.get_min_index() - k;
223 else if (k > coeffs.get_max_index())
224 index = 2 * coeffs.get_max_index() - k;
225 else
226 index = k;
227 assert(coeffs.get_min_index() <= index && index <= coeffs.get_max_index());
228 value += static_cast<T>(coeffs[index] *
229# ifdef NNN
230 f(current_pos, p, bspline)
231 // bspline.function_piece(current_pos, p)
232# else
233 f(current_pos, spline_types[current_dimension])
234# endif
235 );
236 }
237 return value;
238}
239
240template <int num_dimensions, int num_dimensions2, typename T>
241struct compute_BSplines_value
242{
243 typedef T result_type;
244 T operator()(const Array<num_dimensions, T>& coeffs,
245 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
246 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types) const
247 {
248 return spline_convolution(coeffs,
249 relative_positions,
250 spline_types,
251 BW<pos_type>(),
252 // BSplineFunction<quadratic,pos_type>(),
253 compute_BSplines_value<num_dimensions - 1, num_dimensions2, T>());
254 }
255};
256
257template <typename T, int num_dimensions2>
258struct compute_BSplines_value<1, num_dimensions2, T>
259{
260 typedef T result_type;
261 T operator()(const Array<1, T>& coeffs,
262 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
263 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types) const
264 {
265 return spline_convolution(coeffs, relative_positions, spline_types, BW<pos_type>()
266 // BSplineFunction<quadratic,pos_type>()
267 );
268 }
269};
270
271template <int num_dimensions, int num_dimensions2, typename T>
272struct compute_BSplines_gradient
273{
274 typedef BasicCoordinate<num_dimensions, T> result_type;
275
276 BasicCoordinate<num_dimensions, T> operator()(const Array<num_dimensions, T>& coeffs,
277 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
278 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types) const
279 {
280 const T first_value = spline_convolution(coeffs,
281 relative_positions,
282 spline_types,
283 Bder<pos_type>(),
284 // BSplineFunction<quadratic,pos_type>(),
285 compute_BSplines_value<num_dimensions - 1, num_dimensions2, T>());
286 const BasicCoordinate<num_dimensions - 1, T> rest_value
287 = spline_convolution(coeffs,
288 relative_positions,
289 spline_types,
290 BW<pos_type>(),
291 // BSplineFunction<quadratic,pos_type>(),
292 compute_BSplines_gradient<num_dimensions - 1, num_dimensions2, T>());
293 return join(first_value, rest_value);
294 }
295};
296
297template <int num_dimensions2, typename T>
298struct compute_BSplines_gradient<1, num_dimensions2, T>
299{
300 typedef BasicCoordinate<1, T> result_type;
301
302 BasicCoordinate<1, T> operator()(const Array<1, T>& coeffs,
303 const BasicCoordinate<num_dimensions2, pos_type>& relative_positions,
304 const BasicCoordinate<num_dimensions2, BSplineType>& spline_types) const
305 {
306 BasicCoordinate<1, T> result;
307 result[1] = spline_convolution(coeffs, relative_positions, spline_types, Bder<pos_type>()
308 // BSplineFunction<quadratic,pos_type>()
309 );
310 return result;
311 }
312};
313
314} // end of namespace detail
315} // end of namespace BSpline
316
317END_NAMESPACE_STIR
318#endif
assign.h
defines the stir::assign function to assign values to different data types
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::pos_type
double pos_type
The type used for relative positions between the grid points.
Definition BSplines.h:32
stir::BSpline::BSplines_weights
pos_type BSplines_weights(const pos_type relative_position, const BSplineType spline_type)
return spline value
stir::join
BasicCoordinate< num_dimensions+1, coordT > join(const coordT &a, const BasicCoordinate< num_dimensions, coordT > &c)
make a longer BasicCoordinate, by prepending c with the single coordT
Definition BasicCoordinate.inl:447
stir::cut_first_dimension
BasicCoordinate< num_dimensions - 1, coordT > cut_first_dimension(const BasicCoordinate< num_dimensions, coordT > &c)
make a shorter BasicCoordinate, by cutting the first element from c
Definition BasicCoordinate.inl:477
stir::BSpline::BSplines_coef
void BSplines_coef(RandIterOut c_begin_iterator, RandIterOut c_end_iterator, IterT input_begin_iterator, IterT input_end_iterator, const constantsT z1, const constantsT z2, const constantsT lambda)
compute BSpline coefficients that gives the BSpline that interpolates the given data
Definition BSplines_coef.inl:56
Generated on Wed May 27 2026 17:04:05 for STIR by doxygen 1.13.2