File: FFTImageFilterFourierDomainFiltering.cxx

package info (click to toggle)
insighttoolkit4 4.13.3withdata-dfsg1-4
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 489,260 kB
  • sloc: cpp: 557,342; ansic: 146,850; fortran: 34,788; python: 16,572; sh: 2,187; lisp: 2,070; tcl: 993; java: 362; perl: 200; makefile: 129; csh: 81; pascal: 69; xml: 19; ruby: 10
file content (258 lines) | stat: -rw-r--r-- 8,652 bytes parent folder | download | duplicates (5)
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
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  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
 *
 *         http://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.
 *
 *=========================================================================*/

//  Software Guide : BeginLatex
//
//  One of the most common image processing operations performed in the Fourier
//  Domain is the masking of the spectrum in order to eliminate a range of
//  spatial frequencies from the input image. This operation is typically
//  performed by taking the input image, computing its Fourier transform using
//  a FFT filter, masking the resulting image in the Fourier domain with a
//  mask, and finally taking the result of the masking and computing its
//  inverse Fourier transform.
//
//  This typical process is illustrated in the example below.
//
//  \index{itk::Forward\-FFT\-Image\-Filter}
//  \index{itk::Vnl\-Forward\-FFT\-Image\-Filter}
//  \index{itk::FFTW\-Forward\-FFT\-Image\-Filter}
//  \index{itk::Mask\-Image\-Filter}
//
//  Software Guide : EndLatex


#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"

// Software Guide : BeginLatex
//
// We start by including the headers of the FFT filters and the Mask image
// filter. Note that we use two different types of FFT filters here. The first
// one expects as input an image of real pixel type (real in the sense of
// complex numbers) and produces as output a complex image. The second FFT
// filter expects as in put a complex image and produces a real image as
// output.
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet
#include "itkVnlForwardFFTImageFilter.h"
#include "itkVnlInverseFFTImageFilter.h"
#include "itkMaskImageFilter.h"
// Software Guide : EndCodeSnippet


int main( int argc, char * argv [] )
{
  if( argc < 4 )
    {
    std::cerr << "Usage: " << argv[0] << " inputScalarImage  inputMaskImage";
    std::cerr << " outputFilteredImage" << std::endl;
    }

  // Software Guide : BeginLatex
  //
  // The first decision to make is related to the pixel type and dimension of the
  // images on which we want to compute the Fourier transform.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef float  InputPixelType;
  const unsigned int Dimension = 2;

  typedef itk::Image< InputPixelType, Dimension > InputImageType;
  // Software Guide : EndCodeSnippet


  // Software Guide : BeginLatex
  //
  // Then we select the pixel type to use for the mask image and instantiate the
  // image type of the mask.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef unsigned char  MaskPixelType;

  typedef itk::Image< MaskPixelType, Dimension > MaskImageType;
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // Both the input image and the mask image can be read from files or could be
  // obtained as the output of a preprocessing pipeline. We omit here the details
  // of reading the image since the process is quite standard.
  //
  // Software Guide : EndLatex

  typedef itk::ImageFileReader< InputImageType >    InputReaderType;
  typedef itk::ImageFileReader< MaskImageType  >    MaskReaderType;

  InputReaderType::Pointer inputReader = InputReaderType::New();
  MaskReaderType::Pointer  maskReader  = MaskReaderType::New();

  inputReader->SetFileName( argv[1] );
  maskReader->SetFileName( argv[2] );

  // Software Guide : BeginLatex
  //
  // Now the \doxygen{VnlForwardFFTImageFilter} can be instantiated.
  // Like most ITK filters, the FFT filter is instantiated using the full image type.
  // By not setting the output image type, we decide to use the default one provided
  // by the filter. Using this type we construct one instance of the filter.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef itk::VnlForwardFFTImageFilter< InputImageType >  FFTFilterType;

  FFTFilterType::Pointer fftFilter = FFTFilterType::New();

  fftFilter->SetInput( inputReader->GetOutput() );
  // Software Guide : EndCodeSnippet


  // Software Guide : BeginLatex
  //
  // Since our purpose is to perform filtering in the frequency domain by
  // altering the weights of the image spectrum, we need a filter that will
  // mask the Fourier transform of the input image with a binary image. Note that the
  // type of the spectral image is taken here from the traits of the FFT filter.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef FFTFilterType::OutputImageType    SpectralImageType;

  typedef itk::MaskImageFilter< SpectralImageType,
                                    MaskImageType,
                                    SpectralImageType >  MaskFilterType;

  MaskFilterType::Pointer maskFilter = MaskFilterType::New();
  // Software Guide : EndCodeSnippet


  // Software Guide : BeginLatex
  //
  // We connect the inputs to the mask filter by taking the outputs from the
  // first FFT filter and from the reader of the Mask image.
  //
  // Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet
  maskFilter->SetInput1( fftFilter->GetOutput() );
  maskFilter->SetInput2( maskReader->GetOutput() );
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // For the purpose of verifying the aspect of the spectrum after being filtered
  // with the mask, we can write out the output of the Mask filter to a file.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef itk::ImageFileWriter< SpectralImageType > SpectralWriterType;
  SpectralWriterType::Pointer spectralWriter = SpectralWriterType::New();
  spectralWriter->SetFileName("filteredSpectrum.mhd");
  spectralWriter->SetInput( maskFilter->GetOutput() );
  spectralWriter->Update();
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // The output of the mask filter will contain the \emph{filtered} spectrum
  // of the input image. We must then apply an inverse Fourier transform on it in
  // order to obtain the filtered version of the input image. For that purpose we
  // create another instance of the FFT filter.
  //
  // Software Guide : EndLatex


  // Software Guide : BeginCodeSnippet
  typedef itk::VnlInverseFFTImageFilter<
    SpectralImageType >  IFFTFilterType;

  IFFTFilterType::Pointer fftInverseFilter = IFFTFilterType::New();

  fftInverseFilter->SetInput( maskFilter->GetOutput() );
  // Software Guide : EndCodeSnippet

  // Software Guide : BeginLatex
  //
  // The execution of the pipeline can be triggered by invoking the
  // \code{Update()} method in this last filter.  Since this invocation can
  // eventually throw an exception, the call must be placed inside a try/catch
  // block.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  try
    {
    fftInverseFilter->Update();
    }
  catch( itk::ExceptionObject & excp )
    {
    std::cerr << "Error: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
    }
  // Software Guide : EndCodeSnippet


  // Software Guide : BeginLatex
  //
  // The result of the filtering can now be saved into an image file, or be
  // passed to a subsequent processing pipeline. Here we simply write it out to
  // an image file.
  //
  // Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  typedef itk::ImageFileWriter< InputImageType > WriterType;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName( argv[3] );
  writer->SetInput( fftInverseFilter->GetOutput() );
  // Software Guide : EndCodeSnippet


  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excp )
    {
    std::cerr << "Error writing the real image: " << std::endl;
    std::cerr << excp << std::endl;
    return EXIT_FAILURE;
    }

  // Software Guide : BeginLatex
  //
  // Note that this example is just a minimal illustration of the multiple types
  // of processing that are possible in the Fourier domain.
  //
  //  Software Guide : EndLatex

  return EXIT_SUCCESS;
}