Developer documentation
Version 3.0.3-105-gd3941f44
smooth.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 __image_filter_gaussian_h__
18#define __image_filter_gaussian_h__
19
20#include "memory.h"
21#include "image.h"
22#include "algo/copy.h"
23#include "algo/threaded_copy.h"
24#include "adapter/gaussian1D.h"
25#include "filter/base.h"
26
27namespace MR
28{
29 namespace Filter
30 {
47 class Smooth : public Base
48 { MEMALIGN (Smooth)
49
50 public:
51 template <class HeaderType>
52 Smooth (const HeaderType& in) :
53 Base (in),
54 extent (3, 0),
55 stdev (3, 0.0),
56 stride_order (Stride::order (in)),
57 zero_boundary (false)
58 {
59 for (int i = 0; i < 3; i++)
60 stdev[i] = in.spacing(i);
61 datatype() = DataType::Float32;
62 }
63
64 template <class HeaderType>
65 Smooth (const HeaderType& in, const vector<default_type>& stdev_in):
66 Base (in),
67 extent (3, 0),
68 stdev (3, 0.0),
69 stride_order (Stride::order (in))
70 {
71 set_stdev (stdev_in);
72 datatype() = DataType::Float32;
73 }
74
78 void set_extent (const vector<uint32_t>& new_extent)
79 {
80 if (new_extent.size() != 1 && new_extent.size() != 3)
81 throw Exception ("Please supply a single kernel extent value, or three values (one for each spatial dimension)");
82 for (size_t i = 0; i < new_extent.size(); ++i) {
83 if (!(new_extent[i] & uint32_t(1)))
84 throw Exception ("expected odd number for extent");
85 }
86 if (new_extent.size() == 1)
87 for (size_t i = 0; i < 3; i++)
88 extent[i] = new_extent[0];
89 else
90 extent = new_extent;
91 }
92
93 void set_stdev (default_type stdev_in) {
94 set_stdev (vector<default_type> (3, stdev_in));
95 }
96
98 void set_zero_boundary (bool do_zero_boundary) {
99 zero_boundary = do_zero_boundary;
100 }
101
105 void set_stdev (const vector<default_type>& std_dev)
106 {
107 for (size_t i = 0; i < std_dev.size(); ++i)
108 if (stdev[i] < 0.0)
109 throw Exception ("the Gaussian stdev values cannot be negative");
110 if (std_dev.size() == 1) {
111 for (unsigned int i = 0; i < 3; i++)
112 stdev[i] = std_dev[0];
113 } else {
114 if (std_dev.size() != 3)
115 throw Exception ("Please supply a single standard deviation value, or three values (one for each spatial dimension)");
116 for (unsigned int i = 0; i < 3; i++)
117 stdev[i] = std_dev[i];
118 }
119 }
120
122 template <class InputImageType, class OutputImageType, typename ValueType = float>
123 void operator() (InputImageType& input, OutputImageType& output)
124 {
125 std::shared_ptr <Image<ValueType> > in (make_shared<Image<ValueType> > (Image<ValueType>::scratch (input)));
126 threaded_copy (input, *in);
127 std::shared_ptr <Image<ValueType> > out;
128
129 std::unique_ptr<ProgressBar> progress;
130 if (message.size()) {
131 size_t axes_to_smooth = 0;
132 for (vector<default_type>::const_iterator i = stdev.begin(); i != stdev.end(); ++i)
133 if (*i)
134 ++axes_to_smooth;
135 progress.reset (new ProgressBar (message, axes_to_smooth + 1));
136 }
137
138 for (size_t dim = 0; dim < 3; dim++) {
139 if (stdev[dim] > 0) {
140 DEBUG ("creating scratch image for smoothing image along dimension " + str(dim));
141 out = make_shared<Image<ValueType> > (Image<ValueType>::scratch (input));
142 Adapter::Gaussian1D<Image<ValueType> > gaussian (*in, stdev[dim], dim, extent[dim], zero_boundary);
143 threaded_copy (gaussian, *out, 0, input.ndim(), 2);
144 in = out;
145 if (progress)
146 ++(*progress);
147 }
148 }
149 threaded_copy (*in, output);
150 }
151
153 template <class ImageType>
154 void operator() (ImageType& in_and_output)
155 {
156 std::unique_ptr<ProgressBar> progress;
157 if (message.size()) {
158 size_t axes_to_smooth = 0;
159 for (vector<default_type>::const_iterator i = stdev.begin(); i != stdev.end(); ++i)
160 if (*i)
161 ++axes_to_smooth;
162 progress.reset (new ProgressBar (message, axes_to_smooth + 1));
163 }
164
165 for (size_t dim = 0; dim < 3; dim++) {
166 if (stdev[dim] > 0) {
167 vector<size_t> axes (in_and_output.ndim(), dim);
168 size_t axdim = 1;
169 for (size_t i = 0; i < in_and_output.ndim(); ++i) {
170 if (stride_order[i] == dim)
171 continue;
172 axes[axdim++] = stride_order[i];
173 }
174 DEBUG ("smoothing dimension " + str(dim) + " in place with stride order: " + str(axes));
175 SmoothFunctor1D<ImageType> smooth (in_and_output, stdev[dim], dim, extent[dim], zero_boundary);
176 ThreadedLoop (in_and_output, axes, 1).run (smooth, in_and_output);
177 if (progress)
178 ++(*progress);
179 }
180 }
181 }
182
183 protected:
184 vector<uint32_t> extent;
185 vector<default_type> stdev;
186 const vector<size_t> stride_order;
187 bool zero_boundary;
188
189 template <class ImageType>
190 class SmoothFunctor1D { MEMALIGN (SmoothFunctor1D)
191 public:
192 SmoothFunctor1D (ImageType& image,
193 default_type stdev_in = 1.0,
194 size_t axis_in = 0,
195 size_t extent = 0,
196 bool zero_boundary_in = false):
197 stdev (stdev_in),
198 axis (axis_in),
199 zero_boundary (zero_boundary_in),
200 spacing (image.spacing(axis_in)),
201 buffer_size (image.size(axis_in)) {
202 buffer.resize(buffer_size);
203 if (!extent)
204 radius = std::ceil(2 * stdev / spacing);
205 else if (extent == 1)
206 radius = 0;
207 else
208 radius = (extent - 1) / 2;
209 compute_kernel();
210 }
211
212 using value_type = typename ImageType::value_type;
213
214 void compute_kernel() {
215 if ((radius < 1) || stdev <= 0.0)
216 return;
217 kernel.resize (2 * radius + 1);
218 default_type norm_factor = 0.0;
219 for (ssize_t c = 0; c < kernel.size(); ++c) {
220 kernel[c] = exp(-((c-radius) * (c-radius) * spacing * spacing) / (2 * stdev * stdev));
221 norm_factor += kernel[c];
222 }
223 for (ssize_t c = 0; c < kernel.size(); c++) {
224 kernel[c] /= norm_factor;
225 }
226 }
227
228 // SmoothFunctor1D operator():
229 // the inner loop axis has to be the dimension the smoothing is applied to and
230 // the loop has to start with image.index (smooth_axis) == 0
231 void operator () (ImageType& image) {
232 if (!kernel.size())
233 return;
234
235 const ssize_t pos = image.index (axis);
236
237 // fill buffer for current image line if necessary
238 if (pos == 0) {
239 for (ssize_t k = 0; k < buffer_size; ++k) {
240 image.index(axis) = k;
241 buffer(k) = image.value();
242 }
243 image.index (axis) = pos;
244 }
245
246 if (zero_boundary)
247 if (pos == 0 || pos == image.size(axis) - 1) {
248 image.value() = 0.0;
249 return;
250 }
251
252 const ssize_t from = (pos < radius) ? 0 : pos - radius;
253 const ssize_t to = (pos + radius) >= image.size(axis) ? image.size(axis) - 1 : pos + radius;
254
255 ssize_t c = (pos < radius) ? radius - pos : 0;
256 ssize_t kernel_size = to - from + 1;
257
258 value_type result = kernel.segment(c, kernel_size).dot(buffer.segment(from, kernel_size));
259
260 if (!std::isfinite(result)) {
261 result = 0.0;
262 value_type av_weights = 0.0;
263 for (ssize_t k = from; k <= to; ++k, ++c) {
264 value_type neighbour_value = buffer(k);
265 if (std::isfinite (neighbour_value)) {
266 av_weights += kernel[c];
267 result += neighbour_value * kernel[c];
268 }
269 }
270 result /= av_weights;
271 } else if (kernel_size != kernel.size()) {
272 result /= kernel.segment(c, kernel_size).sum();
273 }
274 image.value() = result;
275 }
276
277 private:
278 const default_type stdev;
279 ssize_t radius;
280 size_t axis;
281 Eigen::VectorXd kernel;
282 const bool zero_boundary;
283 const default_type spacing;
284 ssize_t buffer_size;
285 Eigen::VectorXd buffer;
286 };
287 };
289 }
290}
291
292
293#endif
static constexpr uint8_t Float32
Definition: datatype.h:147
std::string message
Definition: base.h:66
vector< default_type > stdev
Definition: smooth.h:185
vector< uint32_t > extent
Definition: smooth.h:184
bool zero_boundary
Definition: smooth.h:187
const vector< size_t > stride_order
Definition: smooth.h:186
static Image scratch(const Header &template_header, const std::string &label="scratch image")
Definition: image.h:195
#define DEBUG(msg)
Definition: exception.h:75
matrix_type stdev(const matrix_type &measurements, const matrix_type &design)
constexpr I ceil(const T x)
template function with cast to different type
Definition: math.h:86
MR::default_type value_type
Definition: typedefs.h:33
vector< size_t > order(const HeaderType &header, size_t from_axis=0, size_t to_axis=std::numeric_limits< size_t >::max())
sort range of axes with respect to their absolute stride.
Definition: stride.h:159
Definition: base.h:24
double default_type
the default type used throughout MRtrix
Definition: types.h:228
std::string str(const T &value, int precision=0)
Definition: mrtrix.h:247
ThreadedLoopRunOuter< decltype(Loop(vector< size_t >()))> ThreadedLoop(const HeaderType &source, const vector< size_t > &outer_axes, const vector< size_t > &inner_axes)
Multi-threaded loop object.
std::shared_ptr< X > make_shared(Args &&... args)
Definition: types.h:283
void threaded_copy(InputImageType &source, OutputImageType &destination, const vector< size_t > &axes, size_t num_axes_in_thread=1)
Definition: threaded_copy.h:43
T to(const std::string &string)
Definition: mrtrix.h:260