1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18 #ifndef itkVnlHalfHermitianToRealInverseFFTImageFilter_hxx
19 #define itkVnlHalfHermitianToRealInverseFFTImageFilter_hxx
20
21 #include "itkImageRegionIteratorWithIndex.h"
22 #include "itkProgressReporter.h"
23 #include "itkVnlHalfHermitianToRealInverseFFTImageFilter.h"
24
25 namespace itk
26 {
27
28 template< typename TInputImage, typename TOutputImage >
29 void
30 VnlHalfHermitianToRealInverseFFTImageFilter< TInputImage, TOutputImage >
GenerateData()31 ::GenerateData()
32 {
33 // Get pointers to the input and output.
34 typename InputImageType::ConstPointer inputPtr = this->GetInput();
35 typename OutputImageType::Pointer outputPtr = this->GetOutput();
36
37 if ( !inputPtr || !outputPtr )
38 {
39 return;
40 }
41
42 // We don't have a nice progress to report, but at least this simple line
43 // reports the beginning and the end of the process.
44 ProgressReporter progress( this, 0, 1 );
45
46 const InputSizeType inputSize = inputPtr->GetLargestPossibleRegion().GetSize();
47 const InputIndexType inputIndex = inputPtr->GetLargestPossibleRegion().GetIndex();
48 const OutputSizeType outputSize = outputPtr->GetLargestPossibleRegion().GetSize();
49 const OutputIndexType outputIndex = outputPtr->GetLargestPossibleRegion().GetIndex();
50
51 // Allocate output buffer memory
52 outputPtr->SetBufferedRegion( outputPtr->GetRequestedRegion() );
53 outputPtr->Allocate();
54
55 unsigned int vectorSize = 1;
56 for ( unsigned int i = 0; i < ImageDimension; i++ )
57 {
58 if ( !VnlFFTCommon::IsDimensionSizeLegal( outputSize[i] ) )
59 {
60 itkExceptionMacro(<< "Cannot compute FFT of image with size "
61 << outputSize << ". VnlHalfHermitianToRealInverseFFTImageFilter operates "
62 << "only on images whose size in each dimension has"
63 << "only a combination of 2,3, and 5 as prime factors." );
64 }
65 vectorSize *= outputSize[i];
66 }
67
68 // VNL requires the full complex result of the transform, so we
69 // produce it here from the half complex image assumed when the output is real.
70 SignalVectorType signal( vectorSize );
71 ImageRegionIteratorWithIndex< OutputImageType > oIt( outputPtr,
72 outputPtr->GetLargestPossibleRegion() );
73
74 OutputIndexValueType maxXIndex = inputIndex[0] +
75 static_cast< OutputIndexValueType >( inputSize[0] );
76 unsigned int si = 0;
77 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
78 {
79 typename OutputImageType::IndexType index = oIt.GetIndex();
80 if ( index[0] >= maxXIndex )
81 {
82 // Flip the indices in each dimension
83 for (unsigned int i = 0; i < ImageDimension; ++i)
84 {
85 if ( index[i] != outputIndex[i] )
86 {
87 index[i] = outputSize[i] - index[i] + 2 * outputIndex[i];
88 }
89 }
90 signal[si] = std::conj( inputPtr->GetPixel( index ) );
91 }
92 else
93 {
94 signal[si] = inputPtr->GetPixel( index );
95 }
96 si++;
97 }
98
99 OutputPixelType *out = outputPtr->GetBufferPointer();
100
101 // call the proper transform, based on compile type template parameter
102 VnlFFTCommon::VnlFFTTransform< OutputImageType > vnlfft( outputSize );
103 vnlfft.transform( signal.data_block(), 1 );
104
105 // Copy the VNL output back to the ITK image. Extract the real part
106 // of the signal. Ideally, the normalization by the number of
107 // elements should have been accounted for by the VNL inverse
108 // Fourier transform, but it is not. So, we take care of it by
109 // dividing the signal by the vectorSize.
110 for ( unsigned int i = 0; i < vectorSize; i++ )
111 {
112 out[i] = signal[i].real() / vectorSize;
113 }
114 }
115
116 template< typename TInputImage, typename TOutputImage >
117 SizeValueType
118 VnlHalfHermitianToRealInverseFFTImageFilter< TInputImage, TOutputImage >
GetSizeGreatestPrimeFactor() const119 ::GetSizeGreatestPrimeFactor() const
120 {
121 return VnlFFTCommon::GREATEST_PRIME_FACTOR;
122 }
123
124
125 }
126 #endif
127