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