Developer documentation
Version 3.0.3-105-gd3941f44
sinc.h
Go to the documentation of this file.
1/* Copyright (c) 2008-2022 the MRtrix3 contributors.
2 *
3 * This Source Code Form is subject to the terms of the Mozilla Public
4 * License, v. 2.0. If a copy of the MPL was not distributed with this
5 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
6 *
7 * Covered Software is provided under this License on an "as is"
8 * basis, without warranty of any kind, either expressed, implied, or
9 * statutory, including, without limitation, warranties that the
10 * Covered Software is free of defects, merchantable, fit for a
11 * particular purpose or non-infringing.
12 * See the Mozilla Public License v. 2.0 for more details.
13 *
14 * For more details, see http://www.mrtrix.org/.
15 */
16
17#ifndef __math_sinc_h__
18#define __math_sinc_h__
19
20#include "math/math.h"
21
22namespace MR
23{
24 namespace Math
25 {
26
27 template <typename T = float> class Sinc
29 public:
30 using value_type = T;
31
32 Sinc (const size_t w) :
33 window_size (w),
34 max_offset_from_kernel_centre ((w-1) / 2),
35 indices (w),
36 weights (w),
37 current_pos (NAN)
38 {
39 assert (w % 2);
40 }
41
42 template <class ImageType>
43 void set (const ImageType& image, const size_t axis, const value_type position) {
44
45 if (position == current_pos)
46 return;
47
48 const int kernel_centre = std::round (position);
49 value_type sum_weights = 0.0;
50
51 for (size_t i = 0; i != window_size; ++i) {
52
53 const int voxel = kernel_centre - max_offset_from_kernel_centre + i;
54 if (voxel < 0)
55 indices[i] = -voxel - 1;
56 else if (voxel >= image.size (axis))
57 indices[i] = (2 * int(image.size (axis))) - voxel - 1;
58 else
59 indices[i] = voxel;
60
61 const value_type offset = position - (value_type)voxel;
62
63 const value_type sinc = offset ? std::sin (Math::pi * offset) / (Math::pi * offset) : 1.0;
64
65 //const value_type hann_cos_term = Math::pi * offset / (value_type(max_offset_from_kernel_centre) + 0.5);
66 //const value_type hann_factor = (abs (hann_cos_term) < Math::pi) ? 0.5 * (1.0 + std::cos (hann_cos_term)) : 0.0;
67 //const value_type this_weight = hann_factor * sinc;
68
69 const value_type lanczos_sinc_term = abs (Math::pi * offset / (double(max_offset_from_kernel_centre) + 0.5));
70 value_type lanczos_factor = 0.0;
71 if (lanczos_sinc_term < Math::pi) {
72 if (lanczos_sinc_term)
73 lanczos_factor = std::sin (lanczos_sinc_term) / lanczos_sinc_term;
74 else
75 lanczos_factor = 1.0;
76 }
77 const value_type this_weight = lanczos_factor * sinc;
78
79 weights[i] = this_weight;
80 sum_weights += this_weight;
81
82 }
83
84 const value_type normalisation = 1.0 / sum_weights;
85 for (size_t i = 0; i != window_size; ++i)
86 weights[i] *= normalisation;
87
88 current_pos = position;
89
90 }
91
92 size_t index (const size_t i) const { return indices[i]; }
93
94 template <class ImageType>
95 value_type value (ImageType& image, const size_t axis) const {
96 assert (current_pos != NAN);
97 const size_t init_pos = image.index(axis);
98 value_type sum = 0.0;
99 for (size_t i = 0; i != window_size; ++i) {
100 image.index(axis) = indices[i];
101 sum += image.value() * weights[i];
102 }
103 image.index(axis) = init_pos;
104 return sum;
105 }
106
107 template <class Cont>
108 value_type value (Cont& data) const {
109 assert (data.size() == window_size);
110 assert (current_pos != NAN);
111 value_type sum = 0.0;
112 for (size_t i = 0; i != window_size; ++i)
113 sum += data[i] * weights[i];
114 return sum;
115 }
116
118 assert (current_pos != NAN);
119 value_type sum = 0.0;
120 for (size_t i = 0; i != window_size; ++i)
121 sum += data[i] * weights[i];
122 return sum;
123 }
124
125 private:
126 const size_t window_size, max_offset_from_kernel_centre;
127 vector<size_t> indices;
128 vector<value_type> weights;
129 value_type current_pos;
130
131 };
132
133 }
134}
135
136#endif
T value_type
Definition: sinc.h:30
Sinc(const size_t w)
Definition: sinc.h:32
void set(const ImageType &image, const size_t axis, const value_type position)
Definition: sinc.h:43
value_type value(value_type *data) const
Definition: sinc.h:117
value_type value(ImageType &image, const size_t axis) const
Definition: sinc.h:95
size_t index(const size_t i) const
Definition: sinc.h:92
value_type value(Cont &data) const
Definition: sinc.h:108
index_type offset
Definition: loop.h:33
constexpr I round(const T x)
Definition: math.h:64
constexpr double pi
Definition: math.h:40
#define NOMEMALIGN
Definition: memory.h:22
Definition: base.h:24
constexpr std::enable_if< std::is_arithmetic< X >::value &&std::is_unsigned< X >::value, X >::type abs(X x)
Definition: types.h:297
int axis