17#ifndef __image_filter_gaussian_h__
18#define __image_filter_gaussian_h__
25#include "filter/base.h"
47 class Smooth :
public Base
51 template <
class HeaderType>
52 Smooth (
const HeaderType& in) :
59 for (
int i = 0; i < 3; i++)
60 stdev[i] = in.spacing(i);
64 template <
class HeaderType>
65 Smooth (
const HeaderType& in,
const vector<default_type>& stdev_in):
78 void set_extent (
const vector<uint32_t>& new_extent)
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");
86 if (new_extent.size() == 1)
87 for (
size_t i = 0; i < 3; i++)
94 set_stdev (vector<default_type> (3, stdev_in));
98 void set_zero_boundary (
bool do_zero_boundary) {
105 void set_stdev (
const vector<default_type>& std_dev)
107 for (
size_t i = 0; i < std_dev.size(); ++i)
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];
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];
122 template <
class InputImageType,
class OutputImageType,
typename ValueType =
float>
123 void operator() (InputImageType& input, OutputImageType& output)
127 std::shared_ptr <Image<ValueType> > out;
129 std::unique_ptr<ProgressBar> progress;
131 size_t axes_to_smooth = 0;
132 for (vector<default_type>::const_iterator i =
stdev.begin(); i !=
stdev.end(); ++i)
135 progress.reset (
new ProgressBar (
message, axes_to_smooth + 1));
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));
153 template <
class ImageType>
154 void operator() (ImageType& in_and_output)
156 std::unique_ptr<ProgressBar> progress;
158 size_t axes_to_smooth = 0;
159 for (vector<default_type>::const_iterator i =
stdev.begin(); i !=
stdev.end(); ++i)
162 progress.reset (
new ProgressBar (
message, axes_to_smooth + 1));
165 for (
size_t dim = 0; dim < 3; dim++) {
166 if (
stdev[dim] > 0) {
167 vector<size_t> axes (in_and_output.ndim(), dim);
169 for (
size_t i = 0; i < in_and_output.ndim(); ++i) {
174 DEBUG (
"smoothing dimension " +
str(dim) +
" in place with stride order: " +
str(axes));
176 ThreadedLoop (in_and_output, axes, 1).run (smooth, in_and_output);
185 vector<default_type>
stdev;
189 template <
class ImageType>
190 class SmoothFunctor1D { MEMALIGN (SmoothFunctor1D)
192 SmoothFunctor1D (ImageType& image,
196 bool zero_boundary_in =
false):
199 zero_boundary (zero_boundary_in),
200 spacing (image.spacing(axis_in)),
201 buffer_size (image.size(axis_in)) {
202 buffer.resize(buffer_size);
208 radius = (
extent - 1) / 2;
214 void compute_kernel() {
215 if ((radius < 1) || stdev <= 0.0)
217 kernel.resize (2 * radius + 1);
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];
223 for (ssize_t c = 0; c < kernel.size(); c++) {
224 kernel[c] /= norm_factor;
231 void operator () (ImageType& image) {
235 const ssize_t pos = image.index (axis);
239 for (ssize_t k = 0; k < buffer_size; ++k) {
240 image.index(axis) = k;
241 buffer(k) = image.value();
243 image.index (axis) = pos;
247 if (pos == 0 || pos == image.size(axis) - 1) {
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;
255 ssize_t c = (pos < radius) ? radius - pos : 0;
256 ssize_t kernel_size =
to - from + 1;
258 value_type result = kernel.segment(c, kernel_size).dot(buffer.segment(from, kernel_size));
260 if (!std::isfinite(result)) {
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];
270 result /= av_weights;
271 }
else if (kernel_size != kernel.size()) {
272 result /= kernel.segment(c, kernel_size).sum();
274 image.value() = result;
281 Eigen::VectorXd kernel;
282 const bool zero_boundary;
285 Eigen::VectorXd buffer;
static constexpr uint8_t Float32
vector< default_type > stdev
vector< uint32_t > extent
const vector< size_t > stride_order
static Image scratch(const Header &template_header, const std::string &label="scratch image")
matrix_type stdev(const matrix_type &measurements, const matrix_type &design)
constexpr I ceil(const T x)
template function with cast to different type
MR::default_type value_type
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.
double default_type
the default type used throughout MRtrix
std::string str(const T &value, int precision=0)
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)
void threaded_copy(InputImageType &source, OutputImageType &destination, const vector< size_t > &axes, size_t num_axes_in_thread=1)
T to(const std::string &string)