Extract Contours From Image#

Synopsis#

Extract contours from an image.

Results#

Output:

There are 2 contours

Contour 0:
[10.5, 19]
[10.4603, 18]
[10.4603, 17]
[10.4603, 16]
[10.4603, 15]
[10.4603, 14]
[10.4603, 13]
[10.4603, 12]
[10.4603, 11]
[10.5, 10]
[10, 9.58579]
[9.5, 10]
[9.5397, 11]
[9.5397, 12]
[9.5397, 13]
[9.5397, 14]
[9.5397, 15]
[9.5397, 16]
[9.5397, 17]
[9.5397, 18]
[9.5, 19]
[10, 19.4142]
[10.5, 19]

Contour 1:
[59.4142, 59]
[59, 58.5]
[58.5, 58]
[58, 57.5]
[57.5, 57]
[57, 56.5]
[56.5, 56]
[56, 55.5]
[55.5, 55]
[55, 54.5]
[54.5, 54]
[54, 53.5]
[53.5, 53]
[53, 52.5]
[52.5, 52]
[52, 51.5]
[51.5, 51]
[51, 50.5]
[50.5, 50]
[50, 49.5]
[49.5, 49]
[49, 48.5]
[48.5, 48]
[48, 47.5]
[47.5, 47]
[47, 46.5]
[46.5, 46]
[46, 45.5]
[45.5, 45]
[45, 44.5]
[44.5, 44]
[44, 43.5]
[43.5, 43]
[43, 42.5]
[42.5, 42]
[42, 41.5]
[41.5, 41]
[41, 40.5]
[40.5, 40]
[40, 39.5858]
[39.5858, 40]
[40, 40.5]
[40.5, 41]
[41, 41.5]
[41.5, 42]
[42, 42.5]
[42.5, 43]
[43, 43.5]
[43.5, 44]
[44, 44.5]
[44.5, 45]
[45, 45.5]
[45.5, 46]
[46, 46.5]
[46.5, 47]
[47, 47.5]
[47.5, 48]
[48, 48.5]
[48.5, 49]
[49, 49.5]
[49.5, 50]
[50, 50.5]
[50.5, 51]
[51, 51.5]
[51.5, 52]
[52, 52.5]
[52.5, 53]
[53, 53.5]
[53.5, 54]
[54, 54.5]
[54.5, 55]
[55, 55.5]
[55.5, 56]
[56, 56.5]
[56.5, 57]
[57, 57.5]
[57.5, 58]
[58, 58.5]
[58.5, 59]
[59, 59.4142]
[59.4142, 59]

Code#

C++#

#include "itkImage.h"
#include "itkImageFileWriter.h"
#include "itkContourExtractor2DImageFilter.h"
#include "itkApproximateSignedDistanceMapImageFilter.h"

using UnsignedCharImageType = itk::Image<unsigned char, 2>;
using FloatImageType = itk::Image<float, 2>;

static void
CreateImage(UnsignedCharImageType::Pointer image);

int
main()
{
  auto image = UnsignedCharImageType::New();
  CreateImage(image);

  using ApproximateSignedDistanceMapImageFilterType =
    itk::ApproximateSignedDistanceMapImageFilter<UnsignedCharImageType, FloatImageType>;
  ApproximateSignedDistanceMapImageFilterType::Pointer approximateSignedDistanceMapImageFilter =
    ApproximateSignedDistanceMapImageFilterType::New();
  approximateSignedDistanceMapImageFilter->SetInput(image);
  approximateSignedDistanceMapImageFilter->SetInsideValue(255);
  approximateSignedDistanceMapImageFilter->SetOutsideValue(0);
  approximateSignedDistanceMapImageFilter->Update();

  using ContourExtractor2DImageFilterType = itk::ContourExtractor2DImageFilter<FloatImageType>;
  auto contourExtractor2DImageFilter = ContourExtractor2DImageFilterType::New();
  contourExtractor2DImageFilter->SetInput(approximateSignedDistanceMapImageFilter->GetOutput());
  contourExtractor2DImageFilter->SetContourValue(0);
  contourExtractor2DImageFilter->Update();

  std::cout << "There are " << contourExtractor2DImageFilter->GetNumberOfOutputs() << " contours" << std::endl;


  for (unsigned int i = 0; i < contourExtractor2DImageFilter->GetNumberOfOutputs(); ++i)
  {
    std::cout << "Contour " << i << ": " << std::endl;
    ContourExtractor2DImageFilterType::VertexListType::ConstIterator vertexIterator =
      contourExtractor2DImageFilter->GetOutput(i)->GetVertexList()->Begin();
    while (vertexIterator != contourExtractor2DImageFilter->GetOutput(i)->GetVertexList()->End())
    {
      std::cout << vertexIterator->Value() << std::endl;
      ++vertexIterator;
    }
    std::cout << std::endl;
  }

  return EXIT_SUCCESS;
}

void
CreateImage(UnsignedCharImageType::Pointer image)
{
  // Create an image
  itk::Index<2> start;
  start.Fill(0);

  itk::Size<2> size;
  size.Fill(100);

  itk::ImageRegion<2> region(start, size);
  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);

  // Create a line of white pixels
  for (unsigned int i = 40; i < 60; ++i)
  {
    itk::Index<2> pixel;
    pixel.Fill(i);
    image->SetPixel(pixel, 255);
  }

  // Create another line of pixels
  for (unsigned int i = 10; i < 20; ++i)
  {
    itk::Index<2> pixel;
    pixel[0] = 10;
    pixel[1] = i;
    image->SetPixel(pixel, 255);
  }

  itk::WriteImage(image, "image.png");
}

Classes demonstrated#

template<typename TInputImage>
class ContourExtractor2DImageFilter : public itk::ImageToPathFilter<TInputImage, PolyLineParametricPath<2>>

Computes a list of PolyLineParametricPath objects from the contours in a 2D image.

Uses the “marching squares” method to compute a the iso-valued contours of the input 2D image for a given intensity value. Multiple outputs may be produced because an image can have multiple contours at a given level, so it is advised to call GetNumberOfIndexedOutputs() and GetOutput(n) to retrieve all of the contours. The contour value to be extracted can be set with SetContourValue(). Image intensities will be linearly interpolated to provide sub-pixel resolution for the output contours.

The marching squares algorithm is a special case of the marching cubes algorithm (Lorensen, William and Harvey E. Cline. Marching Cubes: A High Resolution 3D Surface Construction Algorithm. Computer Graphics (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170). A simple explanation is available here: http://users.polytech.unice.fr/~lingrand/MarchingCubes/algo.html

There is an ambiguous case in the marching squares algorithm: if a given 2x2-pixel square has two high-valued and two low-valued pixels, each pair diagonally adjacent. (Note that high-valued and low-valued are with respect to the contour value sought when LabelContours is false. When LabelContours is true, high-valued means the label being traced and low-valued means any other label.) In this case, the default behavior is that the low-valued pixels are connected into the same contour via an isthmus that separates the high-valued pixels into separate contours. To reverse this, call VertexConnectHighPixelsOn(). Note that when LabelContours is true, the default behavior will leave all four pixels in separate contours. In this case, VertexConnectHighPixels equal to true can instead create contours that are crossing barbells.

Outputs are not guaranteed to be closed paths: contours which intersect the image edge will be left open. All other paths will be closed. (The closedness of a path can be tested by checking whether the beginning point is the same as the end point.)

Produced paths are oriented. Following the path from beginning to end, image intensity values lower than the contour value are to the left of the path and intensity values greater than the contour value are to the right. In other words, the image gradient at a path segment is (approximately) in the direct of that segment rotated right by 90 degrees, because the image intensity values increase from left-to-right across the segment. This means that the generated contours will circle clockwise around “hills” of above-contour-value intensity, and counter-clockwise around “depressions” of below-contour-value intensity. This convention can be reversed by calling ReverseContourOrientationOn().

By default values are interpreted as intensities relative to a contour value. First calling LabelContoursOn() changes this behavior to instead interpret each value as a label. Boundaries are computed for each label separately and all are returned. The value of LabelContours affects the interpretation of VertexConnectHighPixels; see above.

By default the input image’s largest possible region will be processed; call SetRequestedRegion() to process a different region, or ClearRequestedRegion() to revert to the default value. Note that the requested regions are usually set on the output; however since paths have no notion of a “region”, this must be set at the filter level.

This class was contributed to the Insight Journal by Zachary Pincus. https://www.insight-journal.org/browse/publication/72

See

Image

See

Path

See

PolyLineParametricPath

ITK Sphinx Examples:

See itk::ContourExtractor2DImageFilter for additional documentation.