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