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 
19 #include "itkComplexToRealImageFilter.h"
20 #include "itkComplexToRealImageAdaptor.h"
21 #include "itkMath.h"
22 #include "itkSubtractImageFilter.h"
23 #include "itkTestingMacros.h"
24 
itkComplexToRealFilterAndAdaptorTest(int,char * [])25 int itkComplexToRealFilterAndAdaptorTest( int, char* [] )
26 {
27 
28   // Define the dimension of the images
29   constexpr unsigned int ImageDimension = 3;
30 
31   // Declare the types of the images
32   using InputPixelType = std::complex< float >;
33   using OutputPixelType = float;
34 
35   using InputImageType = itk::Image< InputPixelType, ImageDimension >;
36   using OutputImageType = itk::Image< OutputPixelType, ImageDimension >;
37 
38   // Declare appropriate Iterator types for each image
39   using InputIteratorType = itk::ImageRegionIteratorWithIndex<InputImageType>;
40   using OutputIteratorType = itk::ImageRegionIteratorWithIndex<OutputImageType>;
41 
42   // Declare the type of the index to access images
43   using IndexType = itk::Index<ImageDimension>;
44 
45   // Declare the type of the size
46   using SizeType = itk::Size<ImageDimension>;
47 
48   // Declare the type of the Region
49   using RegionType = itk::ImageRegion<ImageDimension>;
50 
51   // Create two images
52   InputImageType::Pointer inputImage  = InputImageType::New();
53 
54   // Define their size, and start index
55   SizeType size;
56   size[0] = 2;
57   size[1] = 2;
58   size[2] = 2;
59 
60   IndexType start;
61   start[0] = 0;
62   start[1] = 0;
63   start[2] = 0;
64 
65   RegionType region;
66   region.SetIndex( start );
67   region.SetSize( size );
68 
69   // Initialize Image A
70   inputImage->SetLargestPossibleRegion( region );
71   inputImage->SetBufferedRegion( region );
72   inputImage->SetRequestedRegion( region );
73   inputImage->Allocate();
74   // Create one iterator for the Input Image (this is a light object)
75   InputIteratorType it( inputImage, inputImage->GetBufferedRegion() );
76 
77   // Initialize the content of Image A
78   InputPixelType value( 13, 25 );
79   it.GoToBegin();
80   while( !it.IsAtEnd() )
81   {
82     it.Set( value );
83     ++it;
84   }
85 
86   // Declare the type for the ComplexToReal filter
87   using FilterType = itk::ComplexToRealImageFilter< InputImageType,
88                                OutputImageType >;
89 
90   // Create the Filter
91   FilterType::Pointer filter = FilterType::New();
92 
93   EXERCISE_BASIC_OBJECT_METHODS( filter, ComplexToRealImageFilter,
94     UnaryGeneratorImageFilter );
95 
96   // Set the input image
97   filter->SetInput( inputImage );
98 
99   // Execute the filter
100   filter->Update();
101 
102   // Get the filter output
103   OutputImageType::Pointer outputImage = filter->GetOutput();
104 
105   // Create an iterator for going through the image output
106   OutputIteratorType ot( outputImage, outputImage->GetRequestedRegion() );
107 
108   // Check the content of the result image
109   const OutputImageType::PixelType epsilon = 1e-6;
110   ot.GoToBegin();
111   it.GoToBegin();
112   while( !ot.IsAtEnd() )
113     {
114     const InputImageType::PixelType  input  = it.Get();
115     const OutputImageType::PixelType output = ot.Get();
116     const OutputImageType::PixelType real  = input.real();
117     if( !itk::Math::FloatAlmostEqual( real, output, 10, epsilon ) )
118       {
119       std::cerr.precision( static_cast< int >( itk::Math::abs( std::log10( epsilon ) ) ) );
120       std::cerr << "Error " << std::endl;
121       std::cerr << " real( " << input << ") = " << real << std::endl;
122       std::cerr << " differs from " << output;
123       std::cerr << " by more than " << epsilon << std::endl;
124       return EXIT_FAILURE;
125       }
126     ++ot;
127     ++it;
128     }
129 
130   //
131   // Test the itk::ComplexToRealImageAdaptor
132   //
133 
134   using AdaptorType = itk::ComplexToRealImageAdaptor< InputImageType,
135                           OutputImageType::PixelType >;
136 
137   AdaptorType::Pointer realAdaptor = AdaptorType::New();
138 
139   EXERCISE_BASIC_OBJECT_METHODS( realAdaptor, ComplexToRealImageAdaptor,
140     ImageAdaptor );
141 
142   realAdaptor->SetImage( inputImage );
143 
144   using DiffFilterType = itk::SubtractImageFilter<
145                         OutputImageType,
146                         AdaptorType,
147                         OutputImageType >;
148 
149   DiffFilterType::Pointer diffFilter = DiffFilterType::New();
150 
151   diffFilter->SetInput1( outputImage );
152   diffFilter->SetInput2( realAdaptor );
153 
154   diffFilter->Update();
155 
156   // Get the filter output
157   OutputImageType::Pointer diffImage = diffFilter->GetOutput();
158 
159   // Check the content of the diff image
160   //
161 
162   // Create an iterator for going through the image output
163   OutputIteratorType dt( diffImage, diffImage->GetRequestedRegion() );
164 
165   dt.GoToBegin();
166   while( !dt.IsAtEnd() )
167     {
168     const OutputImageType::PixelType diff = dt.Get();
169     if( std::fabs( diff ) > epsilon )
170       {
171       std::cerr.precision( static_cast< int >( itk::Math::abs( std::log10( epsilon ) ) ) );
172       std::cerr << "Error comparing results with Adaptors" << std::endl;
173       std::cerr << " difference = " << diff << std::endl;
174       std::cerr << " differs from 0 ";
175       std::cerr << " by more than " << epsilon << std::endl;
176       return EXIT_FAILURE;
177       }
178     ++dt;
179     }
180 
181   return EXIT_SUCCESS;
182 }
183