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 <iostream>
20 #include "itkIndex.h"
21 #include "itkImage.h"
22 #include "itkRGBPixel.h"
23 #include "itkImageFileReader.h"
24 #include "itkImageFileWriter.h"
25 #include "itkImageRegionIterator.h"
26 #include "itkImageDuplicator.h"
27 #include "itkImageRegionIteratorWithIndex.h"
28 #include "itkLineIterator.h"
29 #include "itkMultiThreaderBase.h"
30 #include "itkRegionOfInterestImageFilter.h"
31 #include "itkMaskFeaturePointSelectionFilter.h"
32 #include "itkBlockMatchingImageFilter.h"
33 #include "itkScalarToRGBColormapImageFilter.h"
34 #include "itkTranslationTransform.h"
35 #include "itkResampleImageFilter.h"
36
37
itkBlockMatchingImageFilterTest(int argc,char * argv[])38 int itkBlockMatchingImageFilterTest( int argc, char * argv[] )
39 {
40 if( argc < 2 )
41 {
42 std::cerr << "Usage: " << std::endl;
43 std::cerr << " itkitkBlockMatchingImageFilterTest inputImageFile outputImageFile [Mask File]" << std::endl;
44 return EXIT_FAILURE;
45 }
46
47 constexpr double selectFraction = 0.01;
48
49 using InputPixelType = unsigned char;
50 using OutputPixelType = itk::RGBPixel<InputPixelType>;
51 static constexpr unsigned int Dimension = 3;
52
53 using InputImageType = itk::Image< InputPixelType, Dimension >;
54 using OutputImageType = itk::Image< OutputPixelType, Dimension >;
55
56 // Parameters used for FS and BM
57 using RadiusType = InputImageType::SizeType;
58 RadiusType blockRadius;
59 blockRadius.Fill( 2 );
60
61 RadiusType searchRadius;
62 searchRadius.Fill( 7 );
63
64 using ReaderType = itk::ImageFileReader< InputImageType >;
65
66 //Set up the reader
67 ReaderType::Pointer reader = ReaderType::New();
68 reader->SetFileName( argv[1] );
69 try
70 {
71 reader->Update();
72 }
73 catch( itk::ExceptionObject & e )
74 {
75 std::cerr << "Error in reading the input image: " << e << std::endl;
76 return EXIT_FAILURE;
77 }
78
79 // Reduce region of interest by SEARCH_RADIUS
80 using RegionOfInterestFilterType = itk::RegionOfInterestImageFilter< InputImageType, InputImageType >;
81
82 RegionOfInterestFilterType::Pointer regionOfInterestFilter = RegionOfInterestFilterType::New();
83
84 regionOfInterestFilter->SetInput( reader->GetOutput() );
85
86 RegionOfInterestFilterType::RegionType regionOfInterest = reader->GetOutput()->GetLargestPossibleRegion();
87
88 RegionOfInterestFilterType::RegionType::IndexType regionOfInterestIndex = regionOfInterest.GetIndex();
89 regionOfInterestIndex += searchRadius;
90 regionOfInterest.SetIndex( regionOfInterestIndex );
91
92 RegionOfInterestFilterType::RegionType::SizeType regionOfInterestSize = regionOfInterest.GetSize();
93 regionOfInterestSize -= searchRadius + searchRadius;
94 regionOfInterest.SetSize( regionOfInterestSize );
95
96 regionOfInterestFilter->SetRegionOfInterest( regionOfInterest );
97 regionOfInterestFilter->Update();
98
99 using FeatureSelectionFilterType = itk::MaskFeaturePointSelectionFilter< InputImageType >;
100 using PointSetType = FeatureSelectionFilterType::FeaturePointsType;
101
102 // Feature Selection
103 FeatureSelectionFilterType::Pointer featureSelectionFilter = FeatureSelectionFilterType::New();
104
105 featureSelectionFilter->SetInput( regionOfInterestFilter->GetOutput() );
106 featureSelectionFilter->SetSelectFraction( selectFraction );
107 featureSelectionFilter->SetBlockRadius( blockRadius );
108 featureSelectionFilter->ComputeStructureTensorsOff();
109
110 // Create transformed image from input to match with
111 using TranslationTransformType = itk::TranslationTransform< double, Dimension >;
112 TranslationTransformType::Pointer transform = TranslationTransformType::New();
113 TranslationTransformType::OutputVectorType translation;
114 // move each pixel in input image 5 pixels along first(0) dimension
115 translation[0] = 20.0;
116 translation[1] = 0.0;
117 translation[2] = 0.0;
118 transform->Translate(translation);
119
120 using ResampleImageFilterType = itk::ResampleImageFilter< InputImageType, InputImageType >;
121 ResampleImageFilterType::Pointer resampleFilter = ResampleImageFilterType::New();
122 resampleFilter->SetTransform( transform );
123 resampleFilter->SetInput( reader->GetOutput() );
124 resampleFilter->SetReferenceImage( reader->GetOutput() );
125 resampleFilter->UseReferenceImageOn();
126
127 using BlockMatchingFilterType = itk::BlockMatchingImageFilter< InputImageType >;
128 BlockMatchingFilterType::Pointer blockMatchingFilter = BlockMatchingFilterType::New();
129
130 // inputs (all required)
131 blockMatchingFilter->SetFixedImage( resampleFilter->GetOutput() );
132 blockMatchingFilter->SetMovingImage( reader->GetOutput() );
133 blockMatchingFilter->SetFeaturePoints( featureSelectionFilter->GetOutput() );
134
135 // parameters (all optional)
136 blockMatchingFilter->SetBlockRadius( blockRadius );
137 blockMatchingFilter->SetSearchRadius( searchRadius );
138
139 std::cout << "Block matching: " << blockMatchingFilter << std::endl;
140 try
141 {
142 blockMatchingFilter->Update();
143 }
144 catch ( itk::ExceptionObject &err )
145 {
146 std::cerr << err << std::endl;
147 return EXIT_FAILURE;
148 }
149
150 // Exercise the following methods
151 BlockMatchingFilterType::DisplacementsType * displacements = blockMatchingFilter->GetDisplacements();
152 if( displacements == nullptr )
153 {
154 std::cerr << "GetDisplacements() failed." << std::endl;
155 return EXIT_FAILURE;
156 }
157 BlockMatchingFilterType::SimilaritiesType * similarities = blockMatchingFilter->GetSimilarities();
158 if( similarities == nullptr )
159 {
160 std::cerr << "GetSimilarities() failed." << std::endl;
161 return EXIT_FAILURE;
162 }
163
164 // create RGB copy of input image
165 using RGBFilterType = itk::ScalarToRGBColormapImageFilter< InputImageType, OutputImageType >;
166 RGBFilterType::Pointer colormapImageFilter = RGBFilterType::New();
167
168 colormapImageFilter->SetColormap( RGBFilterType::Grey );
169 colormapImageFilter->SetInput( reader->GetOutput() );
170 try
171 {
172 colormapImageFilter->Update();
173 }
174 catch ( itk::ExceptionObject &err )
175 {
176 std::cerr << err << std::endl;
177 return EXIT_FAILURE;
178 }
179
180 OutputImageType::Pointer outputImage = colormapImageFilter->GetOutput();
181
182 // Highlight the feature points identified in the output image
183 using PointIteratorType = PointSetType::PointsContainer::ConstIterator;
184 using PointDataIteratorType = BlockMatchingFilterType::DisplacementsType::PointDataContainer::ConstIterator;
185
186 PointIteratorType pointItr = featureSelectionFilter->GetOutput()->GetPoints()->Begin();
187 PointIteratorType pointEnd = featureSelectionFilter->GetOutput()->GetPoints()->End();
188 PointDataIteratorType displItr = displacements->GetPointData()->Begin();
189
190 // define colors
191 OutputPixelType red;
192 red.SetRed( 255 );
193 red.SetGreen( 0 );
194 red.SetBlue( 0 );
195
196 OutputPixelType green;
197 green.SetRed( 0 );
198 green.SetGreen( 255 );
199 green.SetBlue( 0 );
200
201 OutputPixelType blue;
202 blue.SetRed( 0 );
203 blue.SetGreen( 0 );
204 blue.SetBlue( 255 );
205
206 OutputImageType::IndexType index;
207 while ( pointItr != pointEnd )
208 {
209 if ( outputImage->TransformPhysicalPointToIndex(pointItr.Value(), index) )
210 {
211 OutputImageType::IndexType displ;
212 outputImage->TransformPhysicalPointToIndex( pointItr.Value() + displItr.Value(), displ );
213
214 // draw line between old and new location of a point in blue
215 itk::LineIterator< OutputImageType > lineIter( outputImage, index, displ );
216 for ( lineIter.GoToBegin(); !lineIter.IsAtEnd(); ++lineIter )
217 {
218 lineIter.Set( blue );
219 }
220
221 // mark old location of a point in green
222 outputImage->SetPixel(index, green);
223
224 // mark new location of a point in red
225 outputImage->SetPixel(displ, red);
226 }
227 pointItr++;
228 displItr++;
229 }
230
231 //Set up the writer
232 using WriterType = itk::ImageFileWriter< OutputImageType >;
233 WriterType::Pointer writer = WriterType::New();
234
235 writer->SetFileName( argv[2] );
236 writer->SetInput( outputImage );
237 try
238 {
239 writer->Update();
240 }
241 catch( itk::ExceptionObject & e )
242 {
243 std::cerr << "Error in writing the output image:" << e << std::endl;
244 return EXIT_FAILURE;
245 }
246
247 return EXIT_SUCCESS;
248 }
249