Skip to content

Commit

Permalink
ENH: Added GetCircles() function
Browse files Browse the repository at this point in the history
  • Loading branch information
Julien Jomier committed Jul 25, 2003
1 parent 47d668c commit 49f7706
Show file tree
Hide file tree
Showing 2 changed files with 200 additions and 15 deletions.
42 changes: 37 additions & 5 deletions Code/BasicFilters/itkHoughTransform2DCirclesImageFilter.h
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
/*=========================================================================
Program: itkUNC
Program: Insight Segmentation & Registration Toolkit
Module: itkHoughTransform2DCirclesImageFilter.h
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) 2002 CADDLab @ UNC. All rights reserved.
See itkUNCCopyright.txt for details.
Copyright (c) 2002 Insight Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
Expand All @@ -29,7 +29,7 @@ namespace itk
{

/**
* \class HoughTransformLines2DImageFilter
* \class HoughTransformCircles2DImageFilter
* \brief Performs the Hough Transform to find circles in a 2D image.
*
* This filter derives from the base class ImageToImageFilter
Expand All @@ -43,6 +43,7 @@ namespace itk
* When the filter found a "correct" point, it computes the gradient at this
* point and draw a regular narrow-banded circle using the minimum and maximum
* radius given by the user, and fill in the array of radii.
* The SweepAngle value can be adjusted to improve the segmentation.
*
* \ingroup ImageFeatureExtraction
*
Expand Down Expand Up @@ -83,6 +84,11 @@ class ITK_EXPORT HoughTransform2DCirclesImageFilter :
/** Typedef to describe the output image region type. */
typedef typename InputImageType::RegionType OutputImageRegionType;

/** Circle typedef */
typedef EllipseSpatialObject<2> CircleType;
typedef typename CircleType::Pointer CirclePointer;
typedef std::list<CirclePointer> CirclesListType;

/** Run-time type information (and related methods). */
itkTypeMacro( HoughTransform2DCirclesImageFilter, ImageToImageFilter );

Expand Down Expand Up @@ -117,6 +123,26 @@ class ITK_EXPORT HoughTransform2DCirclesImageFilter :
/** Get the scale value */
itkGetMacro(SigmaGradient,double);

/** Get the list of circles. This recomputes the circles */
CirclesListType & GetCircles(unsigned int n=0);

/** Set/Get the number of circles to extract */
itkSetMacro(NumberOfCircles,unsigned int);
itkGetMacro(NumberOfCircles,unsigned int);

/** Set/Get the radius of the disc to remove from the accumulator
* for each circle found */
itkSetMacro(DiscRadiusRatio,float);
itkGetMacro(DiscRadiusRatio,float);

/** Set the variance of the gaussian bluring for the accumulator */
itkSetMacro(Variance,float);
itkGetMacro(Variance,float);

/** Set the sweep angle */
itkSetMacro(SweepAngle,float);
itkGetMacro(SweepAngle,float);

protected:

HoughTransform2DCirclesImageFilter();
Expand All @@ -129,12 +155,18 @@ class ITK_EXPORT HoughTransform2DCirclesImageFilter :

private:

float m_SweepAngle;
double m_MinimumRadius;
double m_MaximumRadius;
double m_Threshold;
double m_SigmaGradient;
OutputImagePointer m_RadiusImage;

CirclesListType m_CirclesList;
unsigned int m_NumberOfCircles;
float m_DiscRadiusRatio;
float m_Variance;
unsigned long m_OldModifiedTime;
unsigned long m_OldNumberOfCircles;

};

Expand Down
173 changes: 163 additions & 10 deletions Code/BasicFilters/itkHoughTransform2DCirclesImageFilter.txx
View file
Original file line number Diff line number Diff line change
@@ -1,28 +1,29 @@
/*=========================================================================
Program: itkUNC
Program: Insight Segmentation & Registration Toolkit
Module: itkHoughTransform2DCirclesImageFilter.txx
Language: C++
Date: $Date$
Version: $Revision$
Copyright (c) 2002 CADDLab @ UNC. All rights reserved.
See itkUNCCopyright.txt for details.
Copyright (c) 2002 Insight Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/

#ifndef __itkHoughTransform2DCirclesImageFilter_txx
#define __itkHoughTransform2DCirclesImageFilter_txx

#include "itkHoughTransform2DCirclesImageFilter.h"
#include <itkImageRegionIteratorWithIndex.h>
#include <itkDiscreteGaussianImageFilter.h>
#include "itkGaussianDerivativeImageFunction.h"



#ifndef PI
#define PI 3.1415926535897932384626433832795
#endif
Expand All @@ -37,8 +38,13 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>
{
m_Threshold = 0;
m_MinimumRadius = 0; // by default
m_MaximumRadius = 5; // by default
m_MaximumRadius = 10; // by default
m_SigmaGradient = 1; // Scale of the DoG filter
m_DiscRadiusRatio = 1;
m_Variance = 10;
m_OldModifiedTime = 0;
m_OldNumberOfCircles = 0;
m_SweepAngle = 0.0;
}

template<typename TInputPixelType, typename TOutputPixelType>
Expand All @@ -59,6 +65,15 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>
// Get the input and output pointers
InputImageConstPointer m_InputImage = this->GetInput(0);
OutputImagePointer m_OutputImage = this->GetOutput(0);
/*
typedef fltk::Image2DViewer<float> ViewerType;
ViewerType::Pointer viewer = ViewerType::New();
viewer->SetLabel("Intermediate Image");
viewer->SetImage(this->GetInput(0));
viewer->Show();
Fl::run();
return;*/

// Allocate the output
typename InputImageType::RegionType region;
Expand All @@ -79,7 +94,7 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>

typedef GaussianDerivativeImageFunction<InputImageType> DoGFunctionType;
typename DoGFunctionType::Pointer DoGFunction = DoGFunctionType::New();
DoGFunction->SetInputImage(m_InputImage.GetPointer());
DoGFunction->SetInputImage(m_InputImage);
DoGFunction->SetSigma(m_SigmaGradient);

m_RadiusImage = OutputImageType::New();
Expand Down Expand Up @@ -117,13 +132,17 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>
Vx /= norm;
Vy /= norm;

for(double angle = -m_SweepAngle;angle<=m_SweepAngle;angle+=0.05)
{
double i = m_MinimumRadius;
double distance;

do{

index[0] = (long int)(point[0]+i*Vx);
index[1] = (long int)(point[1]+i*Vy);


index[0] = (long int)(point[0]-i*(Vx*cos(angle)+Vy*sin(angle)));
index[1] = (long int)(point[1]-i*(Vx*sin(angle)+Vy*cos(angle)));

distance = sqrt( (index[1]-point[1])*(index[1]-point[1])
+(index[0]-point[0])*(index[0]-point[0])
Expand All @@ -145,7 +164,9 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>
(index[1]>0) &&
(index[1]<static_cast<typename Index<2>::IndexValueType>(size[1]))
&& (distance < m_MaximumRadius)
);
);
}

}

}
Expand All @@ -154,6 +175,127 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>

}


/** Get the list of circles. This recomputes the circles */
template<typename TInputPixelType, typename TOutputPixelType>
typename HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>::CirclesListType &
HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>
::GetCircles(unsigned int n)
{
if((this->GetMTime() == m_OldModifiedTime) && (n == m_OldNumberOfCircles)) // if the filter has not been updated
{
return m_CirclesList;
}

m_CirclesList.clear();

/** Blur the accumulator in order to find the maximum */
typedef Image<float,2> InternalImageType;

OutputImagePointer outputImage = OutputImageType::New();

typename OutputImageType::RegionType region;
region.SetSize(this->GetOutput(0)->GetLargestPossibleRegion().GetSize());
region.SetIndex(this->GetOutput(0)->GetLargestPossibleRegion().GetIndex());
outputImage->SetRegions( region );
outputImage->SetOrigin(this->GetOutput(0)->GetOrigin());
outputImage->SetSpacing(this->GetOutput(0)->GetSpacing());
outputImage->Allocate();
outputImage->FillBuffer(0);

ImageRegionConstIteratorWithIndex< OutputImageType > image_it( this->GetOutput(0), this->GetOutput(0)->GetRequestedRegion() );
image_it.GoToBegin();

ImageRegionIterator< InternalImageType > it( outputImage, outputImage->GetRequestedRegion() );

while( !image_it.IsAtEnd() )
{
it.Set(image_it.Get());
++image_it;
++it;
}


typedef itk::DiscreteGaussianImageFilter<OutputImageType,InternalImageType> GaussianFilterType;
typename GaussianFilterType::Pointer gaussianFilter = GaussianFilterType::New();

gaussianFilter->SetInput(outputImage); // the output is the accumulator image
double variance[2];
variance[0] = m_Variance;
variance[1] = m_Variance;
gaussianFilter->SetVariance(variance);
gaussianFilter->Update();
InternalImageType::Pointer m_PostProcessImage = gaussianFilter->GetOutput();

typedef itk::MinimumMaximumImageCalculator<InternalImageType> MinMaxCalculatorType;
typename MinMaxCalculatorType::Pointer minMaxCalculator = MinMaxCalculatorType::New();
itk::ImageRegionIterator<InternalImageType> it_input(m_PostProcessImage,m_PostProcessImage->GetLargestPossibleRegion());

itk::Index<2> m_index;

unsigned int circles=0;
bool found;

// Find maxima
do
{
minMaxCalculator->SetImage(m_PostProcessImage);
minMaxCalculator->ComputeMaximum();
InternalImageType::PixelType max = minMaxCalculator->GetMaximum();

found = false;
for(it_input.GoToBegin();!it_input.IsAtEnd();++it_input)
{
if(it_input.Get() == max)
{
// Create a Line Spatial Object
CirclePointer Circle = CircleType::New();
Circle->SetId(circles);
Circle->SetRadius(m_RadiusImage->GetPixel(it_input.GetIndex()));

CircleType::VectorType center;
center[0] = it_input.GetIndex()[0];
center[1] = it_input.GetIndex()[1];
Circle->GetObjectToParentTransform()->SetOffset(center);
Circle->ComputeBoundingBox();

m_CirclesList.push_back(Circle);

// Remove a black disc from the hough space domain
for(double angle = 0; angle <= 2*PI ; angle += PI/1000)
{
for(double lenght = 0; lenght < m_DiscRadiusRatio*Circle->GetRadius()[0];lenght += 1)
{
m_index[0] = (long int)(it_input.GetIndex()[0] + lenght * cos(angle));
m_index[1] = (long int)(it_input.GetIndex()[1] + lenght * sin(angle));
if( ((m_index[0]<=(long)m_PostProcessImage->GetLargestPossibleRegion().GetSize()[0])
&& (m_index[0]>=0)
&& (m_index[1]<=(long)m_PostProcessImage->GetLargestPossibleRegion().GetSize()[1])
&& (m_index[1]>=0)
)
)
{
m_PostProcessImage->SetPixel(m_index,0);
}
}
}
minMaxCalculator->SetImage(m_PostProcessImage);
minMaxCalculator->ComputeMaximum();
max = minMaxCalculator->GetMaximum();

circles++;
found = true;
if(circles == m_NumberOfCircles) break;
}
}
} while((circles<m_NumberOfCircles) && (found));

m_OldModifiedTime = this->GetMTime();
m_OldNumberOfCircles = m_CirclesList.size();
return m_CirclesList;
}


/** Print Self information */
template<typename TInputPixelType, typename TOutputPixelType>
void
Expand All @@ -176,9 +318,20 @@ HoughTransform2DCirclesImageFilter< TInputPixelType, TOutputPixelType>

std::cout << "Radius Image Information : "
<< m_RadiusImage << std::endl;
}

std::cout << "Number Of Circles: "
<< m_NumberOfCircles << std::endl;

std::cout << "Disc Radius: "
<< m_DiscRadiusRatio << std::endl;

std::cout << "Accumulator blur variance: "
<< m_Variance << std::endl;

std::cout << "Sweep angle : "
<< m_SweepAngle << std::endl;

}
} // end namespace

#endif

0 comments on commit 49f7706

Please sign in to comment.