ITK hough检测圆

#include "itkHoughTransform2DCirclesImageFilter.h"
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkImageRegionIterator.h"
#include "itkThresholdImageFilter.h"
#include "itkMinimumMaximumImageCalculator.h"
#include <itkGradientMagnitudeImageFilter.h>
#include <itkDiscreteGaussianImageFilter.h>
#include <list>
#include "itkCastImageFilter.h"
#include "vnl/vnl_math.h"

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageToVTKImageFilter.h"
#include "vtkImageViewer.h"

#include "vtkImageCast.h"

int main( int argc, char *argv[] )
{
    // Software Guide : BeginCodeSnippet
    typedef   unsigned char   PixelType;
    typedef   float           AccumulatorPixelType;
    const     unsigned int    Dimension = 2;
    typedef itk::Image< PixelType, Dimension >  ImageType;
    ImageType::IndexType localIndex;
    typedef itk::Image< AccumulatorPixelType, Dimension > AccumulatorImageType;
    
    typedef  itk::ImageFileReader< ImageType > ReaderType;
    ReaderType::Pointer reader = ReaderType::New();
    reader->SetFileName("/Users/mac/Desktop/CT数据/tsing/image_116" );
    reader->Update();
    
    ImageType::Pointer localImage = reader->GetOutput();
    
    // Software Guide : BeginCodeSnippet
    std::cout << "Computing Hough Map" << std::endl;
    
    typedef itk::HoughTransform2DCirclesImageFilter<PixelType,
    AccumulatorPixelType> HoughTransformFilterType;
    HoughTransformFilterType::Pointer houghFilter = HoughTransformFilterType::New();
    // Software Guide : BeginCodeSnippet
    houghFilter->SetInput( reader->GetOutput() );
    
    houghFilter->SetNumberOfCircles(1);
    houghFilter->SetMinimumRadius(5);
    houghFilter->SetMaximumRadius(45);
    
    houghFilter->Update();
    AccumulatorImageType::Pointer localAccumulator = houghFilter->GetOutput();
    // Software Guide : BeginCodeSnippet
    HoughTransformFilterType::CirclesListType circles;
    circles = houghFilter->GetCircles();
    std::cout << "Found " << circles.size() << " circle(s)." << std::endl;
    
    // Software Guide : BeginCodeSnippet
    typedef  unsigned char                            OutputPixelType;
    typedef  itk::Image< OutputPixelType, Dimension > OutputImageType;
    
    OutputImageType::Pointer  localOutputImage = OutputImageType::New();
    
    OutputImageType::RegionType region;
    region.SetSize(localImage->GetLargestPossibleRegion().GetSize());
    region.SetIndex(localImage->GetLargestPossibleRegion().GetIndex());
    localOutputImage->SetRegions( region );
    localOutputImage->SetOrigin(localImage->GetOrigin());
    localOutputImage->SetSpacing(localImage->GetSpacing());
    localOutputImage->Allocate();
    localOutputImage->FillBuffer(0);
    // Software Guide : BeginCodeSnippet
    typedef HoughTransformFilterType::CirclesListType CirclesListType;
    CirclesListType::const_iterator itCircles = circles.begin();
    
    while( itCircles != circles.end() )
    {
        std::cout << "Center: ";
        std::cout << (*itCircles)->GetObjectToParentTransform()->GetOffset()
        << std::endl;
        std::cout << "Radius: " << (*itCircles)->GetRadius()[0] << std::endl;
        
        // Software Guide : BeginCodeSnippet
        for(double angle = 0;angle <= 2*vnl_math::pi; angle += vnl_math::pi/60.0 )
        {
            localIndex[0] =
            (long int)((*itCircles)->GetObjectToParentTransform()->GetOffset()[0]
                       + (*itCircles)->GetRadius()[0]*vcl_cos(angle));
            localIndex[1] =
            (long int)((*itCircles)->GetObjectToParentTransform()->GetOffset()[1]
                       + (*itCircles)->GetRadius()[0]*vcl_sin(angle));
            OutputImageType::RegionType outputRegion =
            localOutputImage->GetLargestPossibleRegion();
            
            if( outputRegion.IsInside( localIndex ) )
            {
                localOutputImage->SetPixel( localIndex, 255 );
            }
        }
        itCircles++;
    }
    
    
    typedef itk::ImageToVTKImageFilter<ImageType> ConnectorType;
    ConnectorType::Pointer connector= ConnectorType::New();
    connector->SetInput( localOutputImage );
    connector->Update();
    vtkImageViewer2* viewer= vtkImageViewer2::New();
  
  
    viewer->SetInputData( connector->GetOutput() );
    viewer->SetSize(1200,1200);

    ConnectorType::Pointer connector2= ConnectorType::New();
    connector2->SetInput( reader->GetOutput() );
    connector2->Update();
    vtkImageViewer2* viewer2= vtkImageViewer2::New();
    vtkRenderWindowInteractor* renderWindowInteractor2=
    vtkRenderWindowInteractor::New();
    viewer2->SetupInteractor( renderWindowInteractor2);
    viewer->SetupInteractor( renderWindowInteractor2);
    viewer2->SetInputData( connector2->GetOutput() );
    viewer2->SetSize(1200,1200);
    viewer2->Render();
    viewer->Render();
    renderWindowInteractor2->Start();
    
    return 0;
}