HoughTransform2DCirclesImageFilter.cxx

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

// Software Guide : BeginLatex
//
// This example illustrates the use of the
// \doxygen{HoughTransform2DCirclesImageFilter} to find circles in a
// 2-dimensional image.
//
// First, we include the header files of the filter.
//
// Software Guide : EndLatex


// Software Guide : BeginCodeSnippet
#include "itkHoughTransform2DCirclesImageFilter.h"
// Software Guide : EndCodeSnippet

#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 "itkMath.h"

int
main(int argc, char * argv[])
{
  if (argc < 6)
  {
    std::cerr << "Missing Parameters " << std::endl;
    std::cerr << "Usage: " << argv[0] << std::endl;
    std::cerr << " inputImage " << std::endl;
    std::cerr << " outputImage" << std::endl;
    std::cerr << " numberOfCircles " << std::endl;
    std::cerr << " radius Min " << std::endl;
    std::cerr << " radius Max " << std::endl;
    std::cerr << " sweep Angle (default = 0)" << std::endl;
    std::cerr << " SigmaGradient (default = 1) " << std::endl;
    std::cerr << " variance of the accumulator blurring (default = 5) "
              << std::endl;
    std::cerr
      << " radius of the disk to remove from the accumulator (default = 10) "
      << std::endl;
    return EXIT_FAILURE;
  }

  //  Software Guide : BeginLatex
  //
  //  Next, we declare the pixel type and image dimension and specify the
  //  image type to be used as input. We also specify the image type of the
  //  accumulator used in the Hough transform filter.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using PixelType = unsigned char;
  using AccumulatorPixelType = unsigned int;
  using RadiusPixelType = float;
  constexpr unsigned int Dimension = 2;
  using ImageType = itk::Image<PixelType, Dimension>;
  ImageType::IndexType localIndex;
  using AccumulatorImageType = itk::Image<AccumulatorPixelType, Dimension>;
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //
  //  We setup a reader to load the input image.
  //
  //  Software Guide : EndLatex
  // Software Guide : BeginCodeSnippet
  using ReaderType = itk::ImageFileReader<ImageType>;
  auto reader = ReaderType::New();
  reader->SetFileName(argv[1]);
  try
  {
    reader->Update();
  }
  catch (const itk::ExceptionObject & excep)
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    return EXIT_FAILURE;
  }
  ImageType::Pointer localImage = reader->GetOutput();
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  We create the HoughTransform2DCirclesImageFilter based on the pixel
  //  type of the input image (the resulting image from the
  //  ThresholdImageFilter).
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  std::cout << "Computing Hough Map" << std::endl;

  using HoughTransformFilterType =
    itk::HoughTransform2DCirclesImageFilter<PixelType,
                                            AccumulatorPixelType,
                                            RadiusPixelType>;
  auto houghFilter = HoughTransformFilterType::New();
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //
  //  We set the input of the filter to be the output of the
  //  ImageFileReader. We set also the number of circles we are looking for.
  //  Basically, the filter computes the Hough map, blurs it using a certain
  //  variance and finds maxima in the Hough map. After a maximum is found,
  //  the local neighborhood, a circle, is removed from the Hough map.
  //  SetDiscRadiusRatio() defines the radius of this disc proportional to
  //  the radius of the disc found.  The Hough map is computed by looking at
  //  the points above a certain threshold in the input image. Then, for each
  //  point, a Gaussian derivative function is computed to find the direction
  //  of the normal at that point. The standard deviation of the derivative
  //  function can be adjusted by SetSigmaGradient(). The accumulator is
  //  filled by drawing a line along the normal and the length of this line
  //  is defined by the minimum radius (SetMinimumRadius()) and the maximum
  //  radius (SetMaximumRadius()).  Moreover, a sweep angle can be defined by
  //  SetSweepAngle() (default 0.0) to increase the accuracy of detection.
  //
  //  The output of the filter is the accumulator.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  houghFilter->SetInput(reader->GetOutput());

  houghFilter->SetNumberOfCircles(std::stoi(argv[3]));
  houghFilter->SetMinimumRadius(std::stod(argv[4]));
  houghFilter->SetMaximumRadius(std::stod(argv[5]));

  if (argc > 6)
  {
    houghFilter->SetSweepAngle(std::stod(argv[6]));
  }
  if (argc > 7)
  {
    houghFilter->SetSigmaGradient(std::stoi(argv[7]));
  }
  if (argc > 8)
  {
    houghFilter->SetVariance(std::stod(argv[8]));
  }
  if (argc > 9)
  {
    houghFilter->SetDiscRadiusRatio(std::stod(argv[9]));
  }

  houghFilter->Update();
  AccumulatorImageType::Pointer localAccumulator = houghFilter->GetOutput();
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //
  //  We can also get the circles as \doxygen{EllipseSpatialObject}. The
  //  \code{GetCircles()} function return a list of those.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  HoughTransformFilterType::CirclesListType circles;
  circles = houghFilter->GetCircles();
  std::cout << "Found " << circles.size() << " circle(s)." << std::endl;
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //
  //  We can then allocate an image to draw the resulting circles as binary
  //  objects.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;

  auto 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(true); // initializes buffer to zero
  // Software Guide : EndCodeSnippet


  //  Software Guide : BeginLatex
  //
  //  We iterate through the list of circles and we draw them.
  //
  //  Software Guide : EndLatex
  // Software Guide : BeginCodeSnippet
  using CirclesListType = HoughTransformFilterType::CirclesListType;
  CirclesListType::const_iterator itCircles = circles.begin();

  while (itCircles != circles.end())
  {
    std::cout << "Center: ";
    std::cout << (*itCircles)->GetCenterInObjectSpace() << std::endl;
    std::cout << "Radius: " << (*itCircles)->GetRadiusInObjectSpace()[0]
              << std::endl;
    // Software Guide : EndCodeSnippet

    //  Software Guide : BeginLatex
    //
    //  We draw white pixels in the output image to represent each circle.
    //
    //  Software Guide : EndLatex

    // Software Guide : BeginCodeSnippet
    for (double angle = 0; angle <= itk::Math::twopi;
         angle += itk::Math::pi / 60.0)
    {
      const HoughTransformFilterType::CircleType::PointType centerPoint =
        (*itCircles)->GetCenterInObjectSpace();
      using IndexValueType = ImageType::IndexType::IndexValueType;
      localIndex[0] = itk::Math::Round<IndexValueType>(
        centerPoint[0] +
        (*itCircles)->GetRadiusInObjectSpace()[0] * std::cos(angle));
      localIndex[1] = itk::Math::Round<IndexValueType>(
        centerPoint[1] +
        (*itCircles)->GetRadiusInObjectSpace()[0] * std::sin(angle));
      OutputImageType::RegionType outputRegion =
        localOutputImage->GetLargestPossibleRegion();

      if (outputRegion.IsInside(localIndex))
      {
        localOutputImage->SetPixel(localIndex, 255);
      }
    }
    itCircles++;
  }
  // Software Guide : EndCodeSnippet

  //  Software Guide : BeginLatex
  //
  //  We setup a writer to write out the binary image created.
  //
  //  Software Guide : EndLatex

  // Software Guide : BeginCodeSnippet
  using WriterType = itk::ImageFileWriter<ImageType>;
  auto writer = WriterType::New();

  writer->SetFileName(argv[2]);
  writer->SetInput(localOutputImage);

  try
  {
    writer->Update();
  }
  catch (const itk::ExceptionObject & excep)
  {
    std::cerr << "Exception caught !" << std::endl;
    std::cerr << excep << std::endl;
    return EXIT_FAILURE;
  }
  // Software Guide : EndCodeSnippet

  return EXIT_SUCCESS;
}