1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
/* Copyright (c) 2008-2022 the MRtrix3 contributors.
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* Covered Software is provided under this License on an "as is"
* basis, without warranty of any kind, either expressed, implied, or
* statutory, including, without limitation, warranties that the
* Covered Software is free of defects, merchantable, fit for a
* particular purpose or non-infringing.
* See the Mozilla Public License v. 2.0 for more details.
*
* For more details, see http://www.mrtrix.org/.
*/
#ifndef __image_adapter_gaussian1D_h__
#define __image_adapter_gaussian1D_h__
#include "adapter/base.h"
namespace MR
{
namespace Adapter
{
template <class ImageType>
class Gaussian1D :
public Base<Gaussian1D<ImageType>,ImageType>
{ MEMALIGN (Gaussian1D<ImageType>)
public:
using base_type = Base<Gaussian1D<ImageType>, ImageType>;
using value_type = typename ImageType::value_type;
using base_type::name;
using base_type::size;
using base_type::spacing;
using base_type::index;
Gaussian1D (const ImageType& parent,
default_type stdev_in = 1.0,
size_t axis_in = 0,
size_t extent = 0,
bool zero_boundary = false) :
base_type (parent),
stdev (stdev_in),
axis (axis_in),
zero_boundary (zero_boundary) {
if (!extent)
radius = ceil(2 * stdev / spacing(axis));
else if (extent == 1)
radius = 0;
else
radius = (extent - 1) / 2;
compute_kernel();
}
value_type value ()
{
if (!kernel.size())
return base_type::value();
const ssize_t pos = index (axis);
if (zero_boundary)
if (pos == 0 || pos == size(axis) - 1)
return 0.0;
const ssize_t from = (pos < radius) ? 0 : pos - radius;
const ssize_t to = (pos + radius) >= size(axis) ? size(axis) - 1 : pos + radius;
value_type result = 0.0;
value_type av_weights = 0.0;
ssize_t c = (pos < radius) ? radius - pos : 0;
for (ssize_t k = from; k <= to; ++k, ++c) {
index (axis) = k;
value_type neighbour_value = base_type::value();
if (std::isfinite (neighbour_value)) {
av_weights += kernel[c];
result += value_type (base_type::value()) * kernel[c];
}
}
result /= av_weights;
index (axis) = pos;
return result;
}
protected:
void compute_kernel()
{
if ((radius < 1) || stdev <= 0.0)
return;
kernel.resize (2 * radius + 1);
default_type norm_factor = 0.0;
for (size_t c = 0; c < kernel.size(); ++c) {
kernel[c] = exp(-((c-radius) * (c-radius) * spacing(axis) * spacing(axis)) / (2 * stdev * stdev));
norm_factor += kernel[c];
}
for (size_t c = 0; c < kernel.size(); c++) {
kernel[c] /= norm_factor;
}
}
default_type stdev;
ssize_t radius;
size_t axis;
vector<default_type> kernel;
const bool zero_boundary;
};
}
}
#endif
|
