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 itkHoughTransform2DCirclesImageFilter_hxx
19 #define itkHoughTransform2DCirclesImageFilter_hxx
20
21 #include "itkHoughTransform2DCirclesImageFilter.h"
22 #include "itkImageRegionIteratorWithIndex.h"
23 #include "itkDiscreteGaussianImageFilter.h"
24 #include "itkGaussianDerivativeImageFunction.h"
25 #include "itkMinimumMaximumImageCalculator.h"
26 #include "itkMath.h"
27
28 namespace itk
29 {
30 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
31 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
HoughTransform2DCirclesImageFilter()32 ::HoughTransform2DCirclesImageFilter()
33 {
34 this->SetNumberOfRequiredInputs( 1 );
35 }
36
37 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
38 void
39 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
SetRadius(double radius)40 ::SetRadius(double radius)
41 {
42 this->SetMinimumRadius(radius);
43 this->SetMaximumRadius(radius);
44 }
45
46 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
47 void
48 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
EnlargeOutputRequestedRegion(DataObject * output)49 ::EnlargeOutputRequestedRegion(DataObject *output)
50 {
51 Superclass::EnlargeOutputRequestedRegion(output);
52 output->SetRequestedRegionToLargestPossibleRegion();
53 }
54
55 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
56 void
57 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
GenerateInputRequestedRegion()58 ::GenerateInputRequestedRegion()
59 {
60 Superclass::GenerateInputRequestedRegion();
61 if ( this->GetInput() )
62 {
63 InputImagePointer image =
64 const_cast< InputImageType * >( this->GetInput() );
65 image->SetRequestedRegionToLargestPossibleRegion();
66 }
67 }
68
69
70 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
71 void
72 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
VerifyPreconditions()73 ::VerifyPreconditions() ITKv5_CONST
74 {
75 Superclass::VerifyPreconditions();
76
77 if ( ! (m_GradientNormThreshold >= 0.0) )
78 {
79 itkExceptionMacro("Failed precondition: GradientNormThreshold >= 0.");
80 }
81 }
82
83
84 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
85 void
86 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
GenerateData()87 ::GenerateData()
88 {
89 // Get the input and output pointers
90 const InputImageConstPointer inputImage = this->GetInput(0);
91 const OutputImagePointer outputImage = this->GetOutput(0);
92
93 // Allocate the output
94 this->AllocateOutputs();
95 outputImage->FillBuffer(0);
96
97 using DoGFunctionType = GaussianDerivativeImageFunction< InputImageType >;
98 const auto DoGFunction = DoGFunctionType::New();
99 DoGFunction->SetSigma(m_SigmaGradient);
100 DoGFunction->SetUseImageSpacing(m_UseImageSpacing);
101 // Set input image _after_ setting the other GaussianDerivative properties,
102 // to avoid multiple kernel recomputation within GaussianDerivativeImageFunction.
103 DoGFunction->SetInputImage(inputImage);
104
105 m_RadiusImage = RadiusImageType::New();
106
107 m_RadiusImage->SetRegions( outputImage->GetLargestPossibleRegion() );
108 m_RadiusImage->SetOrigin( inputImage->GetOrigin() );
109 m_RadiusImage->SetSpacing( inputImage->GetSpacing() );
110 m_RadiusImage->SetDirection( inputImage->GetDirection() );
111 m_RadiusImage->Allocate( true ); // initialize buffer to zero
112
113 ImageRegionConstIteratorWithIndex< InputImageType > image_it( inputImage,
114 inputImage->GetRequestedRegion() );
115 image_it.GoToBegin();
116
117 const ImageRegion< 2 > & region = outputImage->GetRequestedRegion();
118
119 while ( !image_it.IsAtEnd() )
120 {
121 if ( image_it.Get() > m_Threshold )
122 {
123 const Index< 2 > inputIndex = image_it.GetIndex();
124 const typename DoGFunctionType::VectorType grad =
125 DoGFunction->DoGFunctionType::EvaluateAtIndex( inputIndex );
126
127 double Vx = grad[0];
128 double Vy = grad[1];
129
130 const double norm = std::sqrt(Vx * Vx + Vy * Vy);
131
132 // if the gradient is not flat (using GradientNormThreshold to estimate flatness)
133 if ( norm > m_GradientNormThreshold )
134 {
135 Vx /= norm;
136 Vy /= norm;
137
138 for ( double angle = -m_SweepAngle; angle <= m_SweepAngle; angle += 0.05 )
139 {
140 double i = m_MinimumRadius;
141 double distance;
142
143 do
144 {
145 const Index< 2 > outputIndex =
146 {{
147 Math::Round<IndexValueType>( inputIndex[0] - i * ( Vx * std::cos(angle) + Vy * std::sin(angle) ) ),
148 Math::Round<IndexValueType>( inputIndex[1] - i * ( Vx * std::sin(angle) + Vy * std::cos(angle) ) )
149 }};
150
151 if ( region.IsInside(outputIndex) )
152 {
153 distance = std::sqrt(static_cast<double>((outputIndex[0] - inputIndex[0]) * (outputIndex[0] - inputIndex[0])
154 + (outputIndex[1] - inputIndex[1]) * (outputIndex[1] - inputIndex[1])));
155
156 ++outputImage->GetPixel(outputIndex);
157 m_RadiusImage->GetPixel(outputIndex) += distance;
158 }
159 else
160 {
161 break;
162 }
163 ++i;
164 }
165 while ( distance < m_MaximumRadius );
166 }
167 }
168 }
169 ++image_it;
170 }
171
172 // Compute the average radius
173 ImageRegionConstIterator< OutputImageType > output_it( outputImage,
174 outputImage->GetLargestPossibleRegion() );
175 ImageRegionIterator< RadiusImageType > radius_it( m_RadiusImage,
176 m_RadiusImage->GetLargestPossibleRegion() );
177 output_it.GoToBegin();
178 radius_it.GoToBegin();
179 while ( !output_it.IsAtEnd() )
180 {
181 if ( output_it.Get() > 1 )
182 {
183 radius_it.Value() /= output_it.Get();
184 }
185 ++output_it;
186 ++radius_it;
187 }
188 }
189
190 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
191 typename HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >::CirclesListType &
192 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
GetCircles()193 ::GetCircles()
194 {
195 // Make sure that all the required inputs exist and have a non-null value
196 this->VerifyPreconditions();
197
198 if ( this->GetMTime() == m_OldModifiedTime )
199 {
200 // If the filter has not been updated
201 return m_CirclesList;
202 }
203
204 if( m_RadiusImage.IsNull() )
205 {
206 itkExceptionMacro(<<"Update() must be called before GetCircles().");
207 }
208
209 m_CirclesList.clear();
210
211 if ( m_NumberOfCircles > 0 )
212 {
213 // Blur the accumulator in order to find the maximum
214 using InternalImageType = Image< float, 2 >;
215
216 // The variable "outputImage" is only used as input to gaussianFilter.
217 // It should not be modified, because GetOutput(0) should not be changed.
218 const auto outputImage = OutputImageType::New();
219 outputImage->Graft(this->GetOutput(0));
220
221 const auto gaussianFilter = DiscreteGaussianImageFilter< OutputImageType, InternalImageType >::New();
222
223 gaussianFilter->SetInput(outputImage); // The output is the accumulator image
224 gaussianFilter->SetVariance(m_Variance);
225 gaussianFilter->SetUseImageSpacing(m_UseImageSpacing);
226
227 gaussianFilter->Update();
228 const InternalImageType::Pointer postProcessImage = gaussianFilter->GetOutput();
229
230 const auto minMaxCalculator = MinimumMaximumImageCalculator< InternalImageType >::New();
231 ImageRegionIterator< InternalImageType > it_input( postProcessImage,
232 postProcessImage->GetLargestPossibleRegion() );
233
234 CirclesListSizeType circles = 0;
235
236 // Find maxima
237 // Break out of "forever loop" as soon as the requested number of circles is found.
238 for(;;)
239 {
240 minMaxCalculator->SetImage(postProcessImage);
241 minMaxCalculator->ComputeMaximum();
242
243 if ( minMaxCalculator->GetMaximum() <= 0 )
244 {
245 // When all pixel values in 'postProcessImage' are zero or less, no more circles
246 // should be found. Note that a zero in 'postProcessImage' might correspond to a
247 // zero in the accumulator image, or it might be that the pixel is within a
248 // removed disc around a previously found circle center.
249 break;
250 }
251
252 const InternalImageType::IndexType indexOfMaximum = minMaxCalculator->GetIndexOfMaximum();
253
254 // Create a Circle Spatial Object
255 const auto Circle = CircleType::New();
256 Circle->SetId(static_cast<int>( circles ));
257 Circle->SetRadiusInObjectSpace( m_RadiusImage->GetPixel( indexOfMaximum ) );
258
259 CircleType::PointType center;
260 center[0] = indexOfMaximum[0];
261 center[1] = indexOfMaximum[1];
262 Circle->SetCenterInObjectSpace(center);
263 Circle->Update();
264
265 m_CirclesList.push_back(Circle);
266
267 circles++;
268 if ( circles >= m_NumberOfCircles )
269 {
270 break;
271 }
272
273 // Remove a black disc from the Hough space domain
274 for ( double angle = 0; angle <= 2 * itk::Math::pi; angle += itk::Math::pi / 1000 )
275 {
276 for ( double length = 0; length < m_DiscRadiusRatio *
277 Circle->GetRadiusInObjectSpace()[0]; length += 1 )
278 {
279 const Index< 2 > index =
280 {{
281 Math::Round<IndexValueType>( indexOfMaximum[0] + length * std::cos(angle) ),
282 Math::Round<IndexValueType>( indexOfMaximum[1] + length * std::sin(angle) )
283 }};
284
285 if ( postProcessImage->GetLargestPossibleRegion().IsInside(index) )
286 {
287 postProcessImage->SetPixel(index, 0);
288 }
289 }
290 }
291 }
292 }
293
294 m_OldModifiedTime = this->GetMTime();
295 return m_CirclesList;
296 }
297
298 template< typename TInputPixelType, typename TOutputPixelType, typename TRadiusPixelType >
299 void
300 HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType, TRadiusPixelType >
PrintSelf(std::ostream & os,Indent indent) const301 ::PrintSelf(std::ostream & os, Indent indent) const
302 {
303 Superclass::PrintSelf(os, indent);
304
305 os << indent << "Threshold: " << m_Threshold << std::endl;
306 os << indent << "Gradient Norm Threshold: " << m_GradientNormThreshold << std::endl;
307 os << indent << "Minimum Radius: " << m_MinimumRadius << std::endl;
308 os << indent << "Maximum Radius: " << m_MaximumRadius << std::endl;
309 os << indent << "Derivative Scale : " << m_SigmaGradient << std::endl;
310 os << indent << "Number Of Circles: " << m_NumberOfCircles << std::endl;
311 os << indent << "Disc Radius Ratio: " << m_DiscRadiusRatio << std::endl;
312 os << indent << "Accumulator blur variance: " << m_Variance << std::endl;
313 os << indent << "Sweep angle : " << m_SweepAngle << std::endl;
314 os << indent << "UseImageSpacing: " << m_UseImageSpacing << std::endl;
315
316 itkPrintSelfObjectMacro( RadiusImage );
317
318 os << indent << "CirclesList: " << std::endl;
319 unsigned int i = 0;
320 auto it = m_CirclesList.begin();
321 while( it != m_CirclesList.end() )
322 {
323 os << indent << "[" << i << "]: " << *it << std::endl;
324 ++it;
325 ++i;
326 }
327
328 os << indent << "OldModifiedTime: "
329 << NumericTraits< ModifiedTimeType >::PrintType( m_OldModifiedTime )
330 << std::endl;
331 }
332 } // end namespace
333
334 #endif
335