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 itkVnlComplexToComplexFFTImageFilter_hxx
19 #define itkVnlComplexToComplexFFTImageFilter_hxx
20
21 #include "itkVnlComplexToComplexFFTImageFilter.h"
22 #include "itkProgressReporter.h"
23 #include "itkVnlFFTCommon.h"
24 #include "itkImageRegionIteratorWithIndex.h"
25 #include "itkImageAlgorithm.h"
26
27
28 namespace itk
29 {
30
31 template< typename TImage >
32 VnlComplexToComplexFFTImageFilter< TImage >
VnlComplexToComplexFFTImageFilter()33 ::VnlComplexToComplexFFTImageFilter()
34 {
35 this->DynamicMultiThreadingOn();
36 }
37
38
39 template <typename TImage>
40 void
41 VnlComplexToComplexFFTImageFilter< TImage >
BeforeThreadedGenerateData()42 ::BeforeThreadedGenerateData()
43 {
44 const ImageType * input = this->GetInput();
45 ImageType * output = this->GetOutput();
46
47 const typename ImageType::RegionType bufferedRegion = input->GetBufferedRegion();
48 const typename ImageType::SizeType & imageSize = bufferedRegion.GetSize();
49
50 for( unsigned int ii = 0; ii < ImageDimension; ++ii )
51 {
52 if ( !VnlFFTCommon::IsDimensionSizeLegal( imageSize[ii] ) )
53 {
54 itkExceptionMacro(<< "Cannot compute FFT of image with size "
55 << imageSize << ". VnlComplexToComplexFFTImageFilter operates "
56 << "only on images whose size in each dimension has"
57 << "only a combination of 2,3, and 5 as prime factors." );
58 }
59 }
60
61 // Copy the input to the output, and we will work in place on the output.
62 ImageAlgorithm::Copy< ImageType, ImageType >( input, output, bufferedRegion, bufferedRegion );
63
64 using VclPixelType = std::complex< typename PixelType::value_type >;
65 auto * outputBuffer = static_cast< VclPixelType * >( output->GetBufferPointer() );
66
67 // call the proper transform, based on compile type template parameter
68 VnlFFTCommon::VnlFFTTransform< Image< typename PixelType::value_type , ImageDimension > > vnlfft( imageSize );
69 if ( this->GetTransformDirection() == Superclass::INVERSE )
70 {
71 vnlfft.transform( outputBuffer, 1 );
72 }
73 else
74 {
75 vnlfft.transform( outputBuffer, -1 );
76 }
77 }
78
79
80 template <typename TImage>
81 void
82 VnlComplexToComplexFFTImageFilter< TImage >
DynamicThreadedGenerateData(const OutputImageRegionType & outputRegionForThread)83 ::DynamicThreadedGenerateData(const OutputImageRegionType& outputRegionForThread)
84 {
85 // Normalize the output if backward transform
86 if ( this->GetTransformDirection() == Superclass::INVERSE )
87 {
88 using IteratorType = ImageRegionIterator< OutputImageType >;
89 SizeValueType totalOutputSize = this->GetOutput()->GetRequestedRegion().GetNumberOfPixels();
90 IteratorType it(this->GetOutput(), outputRegionForThread);
91 while( !it.IsAtEnd() )
92 {
93 PixelType val = it.Value();
94 val /= totalOutputSize;
95 it.Set(val);
96 ++it;
97 }
98 }
99 }
100
101 } // end namespace itk
102
103 #endif // _itkVnlComplexToComplexFFTImageFilter_hxx
104