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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
|
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
#ifndef itkFFTConvolutionImageFilter_h
#define itkFFTConvolutionImageFilter_h
#include "itkConvolutionImageFilterBase.h"
#include "itkProgressAccumulator.h"
#include "itkHalfHermitianToRealInverseFFTImageFilter.h"
#include "itkRealToHalfHermitianForwardFFTImageFilter.h"
#include "itkZeroFluxNeumannBoundaryCondition.h"
namespace itk
{
/**
* \class FFTConvolutionImageFilter
* \brief Convolve a given image with an arbitrary image kernel using
* multiplication in the Fourier domain.
*
* This filter produces output equivalent to the output of the
* ConvolutionImageFilter. However, it takes advantage of the
* convolution theorem to accelerate the convolution computation when
* the kernel is large.
*
* \warning This filter ignores the spacing, origin, and orientation
* of the kernel image and treats them as identical to those in the
* input image.
*
* This code was adapted from the Insight Journal contribution:
*
* "FFT Based Convolution"
* by Gaetan Lehmann
* https://www.insight-journal.org/browse/publication/717
*
* \ingroup ITKConvolution
* \sa ConvolutionImageFilter
* \sa InverseDeconvolutionImageFilter
* \sa IterativeDeconvolutionImageFilter
*
*/
template <typename TInputImage,
typename TKernelImage = TInputImage,
typename TOutputImage = TInputImage,
typename TInternalPrecision = double>
class ITK_TEMPLATE_EXPORT FFTConvolutionImageFilter
: public ConvolutionImageFilterBase<TInputImage, TKernelImage, TOutputImage>
{
public:
ITK_DISALLOW_COPY_AND_MOVE(FFTConvolutionImageFilter);
using Self = FFTConvolutionImageFilter;
using Superclass = ConvolutionImageFilterBase<TInputImage, TKernelImage, TOutputImage>;
using Pointer = SmartPointer<Self>;
using ConstPointer = SmartPointer<const Self>;
/** Method for creation through the object factory. */
itkNewMacro(Self);
/** \see LightObject::GetNameOfClass() */
itkOverrideGetNameOfClassMacro(FFTConvolutionImageFilter);
/** Dimensionality of input and output data is assumed to be the same. */
static constexpr unsigned int ImageDimension = TInputImage::ImageDimension;
using InputImageType = TInputImage;
using OutputImageType = TOutputImage;
using KernelImageType = TKernelImage;
using InputPixelType = typename InputImageType::PixelType;
using OutputPixelType = typename OutputImageType::PixelType;
using KernelPixelType = typename KernelImageType::PixelType;
using InputIndexType = typename InputImageType::IndexType;
using OutputIndexType = typename OutputImageType::IndexType;
using KernelIndexType = typename KernelImageType::IndexType;
using InputSizeType = typename InputImageType::SizeType;
using OutputSizeType = typename OutputImageType::SizeType;
using KernelSizeType = typename KernelImageType::SizeType;
using SizeValueType = typename InputSizeType::SizeValueType;
using InputRegionType = typename InputImageType::RegionType;
using OutputRegionType = typename OutputImageType::RegionType;
using KernelRegionType = typename KernelImageType::RegionType;
/** Internal types used by the FFT filters. */
using InternalImageType = Image<TInternalPrecision, TInputImage::ImageDimension>;
using InternalRegionType = typename InternalImageType::RegionType;
using InternalSizeType = typename InternalImageType::SizeType;
using InternalIndexType = typename InternalImageType::IndexType;
using InternalImagePointerType = typename InternalImageType::Pointer;
using InternalComplexType = std::complex<TInternalPrecision>;
using InternalComplexImageType = Image<InternalComplexType, TInputImage::ImageDimension>;
using InternalComplexImagePointerType = typename InternalComplexImageType::Pointer;
/** Typedef to describe the boundary condition. */
using typename Superclass::BoundaryConditionType;
itkSetMacro(SizeGreatestPrimeFactor, SizeValueType);
itkGetMacro(SizeGreatestPrimeFactor, SizeValueType);
protected:
FFTConvolutionImageFilter();
~FFTConvolutionImageFilter() override = default;
/** Because the inputs are real, we can use the specialized filters
* for real-to-complex Fourier transforms. */
using FFTFilterType = RealToHalfHermitianForwardFFTImageFilter<InternalImageType, InternalComplexImageType>;
using IFFTFilterType = HalfHermitianToRealInverseFFTImageFilter<InternalComplexImageType, InternalImageType>;
/** Convolution uses a spatial region equivalent to the
* output region padded by the kernel radius on all sides.
* The input requested region is expanded by the kernel radius
* within the bounds of the input largest possible region.
*
* \sa ProcessObject::GenerateInputRequestedRegion() */
void
GenerateInputRequestedRegion() override;
/** This filter uses a minipipeline to compute the output. */
void
GenerateData() override;
/** Prepare the input images for operations in the Fourier
* domain. This includes resizing the input and kernel images,
* normalizing the kernel if requested, shifting the kernel, and
* taking the Fourier transform of the padded inputs. */
void
PrepareInputs(const InputImageType * input,
const KernelImageType * kernel,
InternalComplexImagePointerType & preparedInput,
InternalComplexImagePointerType & preparedKernel,
ProgressAccumulator * progress,
float progressWeight);
/** Prepare the input image. This includes padding the image and
* taking the Fourier transform of the padded image. */
void
PrepareInput(const InputImageType * input,
InternalComplexImagePointerType & preparedInput,
ProgressAccumulator * progress,
float progressWeight);
/** Pad the input image. */
void
PadInput(const InputImageType * input,
InternalImagePointerType & paddedInput,
ProgressAccumulator * progress,
float progressWeight);
/** Take the Fourier transform of the padded input. */
void
TransformPaddedInput(const InternalImageType * paddedInput,
InternalComplexImagePointerType & transformedInput,
ProgressAccumulator * progress,
float progressWeight);
/** Prepare the kernel. This includes resizing the input and kernel
* images, normalizing the kernel if requested, shifting the kernel,
* and taking the Fourier transform of the padded kernel. */
void
PrepareKernel(const KernelImageType * kernel,
InternalComplexImagePointerType & preparedKernel,
ProgressAccumulator * progress,
float progressWeight);
/** Produce output from the final Fourier domain image. */
void
ProduceOutput(InternalComplexImageType * paddedOutput, ProgressAccumulator * progress, float progressWeight);
/** Crop the padded version of the output. */
void
CropOutput(InternalImageType * paddedOutput, ProgressAccumulator * progress, float progressWeight);
/** Get the radius of the kernel image. Used to pad the input image
* for convolution. */
KernelSizeType
GetKernelRadius() const;
/** Get padding around the region of interest that results from FFT
* factoring requirements. FFT typically requires that image side lengths
* are factorable only by a fixed set of prime numbers (often 2, 3, and 5).
* After the input image is padded for the kernel width and cropped to the
* region of interest the result is then padded for FFT execution. This value
* is reused for kernel padding and output cropping. */
InternalSizeType
GetFFTPadSize() const;
/** Get whether the X dimension has an odd size. */
bool
GetXDimensionIsOdd() const;
void
PrintSelf(std::ostream & os, Indent indent) const override;
private:
SizeValueType m_SizeGreatestPrimeFactor{};
InternalSizeType m_FFTPadSize{ { 0 } };
InternalRegionType m_PaddedInputRegion{};
};
} // namespace itk
#ifndef ITK_MANUAL_INSTANTIATION
# include "itkFFTConvolutionImageFilter.hxx"
#endif
#endif
|
