Examples/Segmentation/HoughTransform2DLinesImageFilter.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{HoughTransform2DLinesImageFilter} to find straight lines in a
// 2-dimensional image.
//
// First, we include the header files of the filter.
//
// Software Guide : EndLatex
 
 
// Software Guide : BeginCodeSnippet
#include "itkHoughTransform2DLinesImageFilter.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 "itkCastImageFilter.h"
 
int
main(int argc, char * argv[])
{
  if (argc < 4)
  {
    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 << " numberOfLines " << 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 = float;
  constexpr unsigned int Dimension = 2;
 
  using ImageType = itk::Image<PixelType, Dimension>;
  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;
  }
  const ImageType::Pointer localImage = reader->GetOutput();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  Once the image is loaded, we apply a
  //  \doxygen{GradientMagnitudeImageFilter} to segment edges.  This casts
  //  the input image using a \doxygen{CastImageFilter}.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using CastingFilterType =
    itk::CastImageFilter<ImageType, AccumulatorImageType>;
  auto caster = CastingFilterType::New();
 
  std::cout << "Applying gradient magnitude filter" << std::endl;
 
  using GradientFilterType =
    itk::GradientMagnitudeImageFilter<AccumulatorImageType,
                                      AccumulatorImageType>;
  auto gradFilter = GradientFilterType::New();
 
  caster->SetInput(localImage);
  gradFilter->SetInput(caster->GetOutput());
  gradFilter->Update();
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  The next step is to apply a threshold filter on the gradient magnitude
  //  image to keep only bright values. Only pixels with a high value will be
  //  used by the Hough transform filter.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  std::cout << "Thresholding" << std::endl;
  using ThresholdFilterType = itk::ThresholdImageFilter<AccumulatorImageType>;
  auto threshFilter = ThresholdFilterType::New();
 
  threshFilter->SetInput(gradFilter->GetOutput());
  threshFilter->SetOutsideValue(0);
  constexpr unsigned char threshBelow = 0;
  constexpr unsigned char threshAbove = 255;
  threshFilter->ThresholdOutside(threshBelow, threshAbove);
  threshFilter->Update();
  // Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  We create the HoughTransform2DLinesImageFilter 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::HoughTransform2DLinesImageFilter<AccumulatorPixelType,
                                          AccumulatorPixelType>;
 
  auto houghFilter = HoughTransformFilterType::New();
  // Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  We set the input to the filter to be the output of the
  //  ThresholdImageFilter. We set also the number of lines 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. SetDiscRadius() defines the radius of this disc.
  //
  //  The output of the filter is the accumulator.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  houghFilter->SetInput(threshFilter->GetOutput());
  houghFilter->SetNumberOfLines(std::stoi(argv[3]));
 
  if (argc > 4)
  {
    houghFilter->SetVariance(std::stod(argv[4]));
  }
 
  if (argc > 5)
  {
    houghFilter->SetDiscRadius(std::stod(argv[5]));
  }
  houghFilter->Update();
  const AccumulatorImageType::Pointer localAccumulator =
    houghFilter->GetOutput();
  // Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  We can also get the lines as \doxygen{LineSpatialObject}. The
  //  \code{GetLines()} function return a list of those.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  HoughTransformFilterType::LinesListType lines;
  lines = houghFilter->GetLines();
  std::cout << "Found " << lines.size() << " line(s)." << std::endl;
  // Software Guide : EndCodeSnippet
 
  //  Software Guide : BeginLatex
  //
  //  We can then allocate an image to draw the resulting lines as binary
  //  objects.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using OutputPixelType = unsigned char;
  using OutputImageType = itk::Image<OutputPixelType, Dimension>;
 
  auto localOutputImage = OutputImageType::New();
 
  const OutputImageType::RegionType region(
    localImage->GetLargestPossibleRegion());
  localOutputImage->SetRegions(region);
  localOutputImage->CopyInformation(localImage);
  localOutputImage->Allocate(true); // initialize buffer to zero
  // Software Guide : EndCodeSnippet
 
 
  //  Software Guide : BeginLatex
  //
  //  We iterate through the list of lines and we draw them.
  //
  //  Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  using LineIterator =
    HoughTransformFilterType::LinesListType::const_iterator;
  LineIterator itLines = lines.begin();
  while (itLines != lines.end())
  {
    // Software Guide : EndCodeSnippet
 
    //  Software Guide : BeginLatex
    //
    //  We get the list of points which consists of two points to represent a
    //  straight line.  Then, from these two points, we compute a fixed point
    //  $u$ and a unit vector $\vec{v}$ to parameterize the line.
    //
    //  Software Guide : EndLatex
 
    // Software Guide : BeginCodeSnippet
    using LinePointListType =
      HoughTransformFilterType::LineType::LinePointListType;
 
    LinePointListType                 pointsList = (*itLines)->GetPoints();
    LinePointListType::const_iterator itPoints = pointsList.begin();
 
    double u[2];
    u[0] = (*itPoints).GetPositionInObjectSpace()[0];
    u[1] = (*itPoints).GetPositionInObjectSpace()[1];
    itPoints++;
    double v[2];
    v[0] = u[0] - (*itPoints).GetPositionInObjectSpace()[0];
    v[1] = u[1] - (*itPoints).GetPositionInObjectSpace()[1];
 
    const double norm = std::sqrt(v[0] * v[0] + v[1] * v[1]);
    v[0] /= norm;
    v[1] /= norm;
    // Software Guide : EndCodeSnippet
 
    //  Software Guide : BeginLatex
    //
    //  We draw a white pixels in the output image to represent the line.
    //
    //  Software Guide : EndLatex
 
    // Software Guide : BeginCodeSnippet
    ImageType::IndexType localIndex;
    itk::Size<2>         size =
      localOutputImage->GetLargestPossibleRegion().GetSize();
    const float diag =
      std::sqrt(static_cast<float>(size[0] * size[0] + size[1] * size[1]));
 
    for (auto i = static_cast<int>(-diag); i < static_cast<int>(diag); ++i)
    {
      localIndex[0] = static_cast<long>(u[0] + i * v[0]);
      localIndex[1] = static_cast<long>(u[1] + i * v[1]);
 
      const OutputImageType::RegionType outputRegion =
        localOutputImage->GetLargestPossibleRegion();
 
      if (outputRegion.IsInside(localIndex))
      {
        localOutputImage->SetPixel(localIndex, 255);
      }
    }
    itLines++;
  }
  // 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<OutputImageType>;
  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;
}