Download this example

Geometric Properties of Labeled Region

Synopsis

Get geometric properties of labeled regions in an image.

Results

Output In QuickView

Output:

Image (0x7fd274c5f3c0)
RTTI typeinfo:   itk::Image<unsigned int, 2u>
Reference Count: 1
Modified Time: 46
Debug: Off
Object Name:
Observers:
  none
Source: (none)
Source output name: (none)
Release Data: Off
Data Released: False
Global Release Data: Off
PipelineMTime: 0
UpdateMTime: 0
RealTimeStamp: 0 seconds
LargestPossibleRegion:
  Dimension: 2
  Index: [0, 0]
  Size: [20, 20]
BufferedRegion:
  Dimension: 2
  Index: [0, 0]
  Size: [20, 20]
RequestedRegion:
  Dimension: 2
  Index: [0, 0]
  Size: [20, 20]
Spacing: [1, 1]
Origin: [0, 0]
Direction:
1 0
0 1

IndexToPointMatrix:
1 0
0 1

PointToIndexMatrix:
1 0
0 1

Inverse Direction:
1 0
0 1

PixelContainer:
  ImportImageContainer (0x7fd274c5f5b0)
    RTTI typeinfo:   itk::ImportImageContainer<unsigned long, unsigned int>
    Reference Count: 1
    Modified Time: 47
    Debug: Off
    Object Name:
    Observers:
      none
    Pointer: 0x7fd27502d800
    Container manages memory: true
    Size: 400
    Capacity: 400
Number of labels: 4

Label: 0
Volume: 344
Integrated Intensity: 44128
Centroid: [9.06977, 9.54942]
Weighted Centroid: [9.44273, 9.96021]
Axes Length: [23.6173, 23.9599]
MajorAxisLength: 23.9599
MinorAxisLength: 23.6173
Eccentricity: 0.168492
Elongation: 1.0145
Orientation: 2.74768
Bounding box: [0, 19, 0, 19]


Label: 85
Volume: 16
Integrated Intensity: 2544
Centroid: [7.5, 7.5]
Weighted Centroid: [7.44811, 7.47602]
Axes Length: [4.6188, 4.6188]
MajorAxisLength: 4.6188
MinorAxisLength: 4.6188
Eccentricity: 0
Elongation: 1
Orientation: 1.5708
Bounding box: [6, 9, 6, 9]

Label: 127
Volume: 25
Integrated Intensity: 3263
Centroid: [14, 14]
Weighted Centroid: [13.6883, 14.1192]
Axes Length: [5.7735, 5.7735]
MajorAxisLength: 5.7735
MinorAxisLength: 5.7735
Eccentricity: 0
Elongation: 1
Orientation: 1.5708
Bounding box: [12, 16, 12, 16]

Label: 191
Volume: 15
Integrated Intensity: 1840
Centroid: [14, 3]
Weighted Centroid: [13.8647, 3.10978]
Axes Length: [3.4641, 5.7735]
MajorAxisLength: 5.7735
MinorAxisLength: 3.4641
Eccentricity: 0.8
Elongation: 1.66667
Orientation: 0
Bounding box: [12, 16, 2, 4]

Code

C++

#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkImageFileReader.h"
#include "itkLabelGeometryImageFilter.h"
#include "itkLabelToRGBImageFilter.h"

#ifdef ENABLE_QUICKVIEW
#  include "QuickView.h"
#endif

#include <sstream>

using ImageType = itk::Image<unsigned int, 2>;
using RGBPixelType = itk::RGBPixel<unsigned char>;
using RGBImageType = itk::Image<RGBPixelType, 2>;

static void
CreateIntensityImage(ImageType::Pointer image);
static void
CreateLabelImage(ImageType::Pointer image);

int
main(int argc, char * argv[])
{
  auto labelImage = ImageType::New();
  auto intensityImage = ImageType::New();
  int  label = 1;

  if (argc < 2)
  {
    // Create a label image that is 0 in the background and where the
    // objects are labeled
    CreateLabelImage(labelImage);
    labelImage->Print(std::cout);
    // Create an intensity image.  Some LabelGeometry features (such as
    // weighted centroid and integrated intensity) depend on an
    // intensity image.
    CreateIntensityImage(intensityImage);
  }
  else if (argc > 3)
  {
    labelImage = itk::ReadImage<ImageType>(argv[1]);

    intensityImage = itk::ReadImage<ImageType>(argv[2]);

    label = std::stoi(argv[3]);
  }

  // NOTE: As of April 8, 2015 the filter does not work with non-zero
  // origins
  double origin[2] = { 0.0, 0.0 };
  labelImage->SetOrigin(origin);
  intensityImage->SetOrigin(origin);

  using LabelGeometryImageFilterType = itk::LabelGeometryImageFilter<ImageType>;
  auto labelGeometryImageFilter = LabelGeometryImageFilterType::New();
  labelGeometryImageFilter->SetInput(labelImage);
  labelGeometryImageFilter->SetIntensityInput(intensityImage);

  // These generate optional outputs.
  labelGeometryImageFilter->CalculatePixelIndicesOn();
  labelGeometryImageFilter->CalculateOrientedBoundingBoxOn();
  labelGeometryImageFilter->CalculateOrientedLabelRegionsOn();
  labelGeometryImageFilter->CalculateOrientedIntensityRegionsOn();

  labelGeometryImageFilter->Update();
  LabelGeometryImageFilterType::LabelsType allLabels = labelGeometryImageFilter->GetLabels();

  using RGBFilterType = itk::LabelToRGBImageFilter<ImageType, RGBImageType>;
  auto rgbLabelImage = RGBFilterType::New();
  rgbLabelImage->SetInput(labelImage);

  using RGBFilterType = itk::LabelToRGBImageFilter<ImageType, RGBImageType>;
  auto rgbOrientedImage = RGBFilterType::New();
  rgbOrientedImage->SetInput(labelGeometryImageFilter->GetOrientedLabelImage(allLabels[label]));

#ifdef ENABLE_QUICKVIEW
  QuickView viewer;
  viewer.ShareCameraOff();
  viewer.InterpolateOff();

  viewer.AddImage(rgbLabelImage->GetOutput(),
                  true,
                  argc > 3 ? itksys::SystemTools::GetFilenameName(argv[1]) : "Generated label image");
  viewer.AddImage(intensityImage.GetPointer(),
                  true,
                  argc > 3 ? itksys::SystemTools::GetFilenameName(argv[2]) : "Generated intensity image");

  std::stringstream desc;
  desc << "Oriented Label: " << allLabels[label];
  viewer.AddImage(rgbOrientedImage->GetOutput(), true, desc.str());

  std::stringstream desc2;
  desc2 << "Oriented Intensity: " << allLabels[label];
  viewer.AddImage(labelGeometryImageFilter->GetOrientedIntensityImage(allLabels[label]), true, desc2.str());
  viewer.Visualize();
#endif

  LabelGeometryImageFilterType::LabelsType::iterator allLabelsIt;
  std::cout << "Number of labels: " << labelGeometryImageFilter->GetNumberOfLabels() << std::endl;
  std::cout << std::endl;

  for (allLabelsIt = allLabels.begin(); allLabelsIt != allLabels.end(); ++allLabelsIt)
  {
    LabelGeometryImageFilterType::LabelPixelType labelValue = *allLabelsIt;
    std::cout << "\tLabel: " << (int)labelValue << std::endl;
    std::cout << "\tVolume: " << labelGeometryImageFilter->GetVolume(labelValue) << std::endl;
    std::cout << "\tIntegrated Intensity: " << labelGeometryImageFilter->GetIntegratedIntensity(labelValue)
              << std::endl;
    std::cout << "\tCentroid: " << labelGeometryImageFilter->GetCentroid(labelValue) << std::endl;
    std::cout << "\tWeighted Centroid: " << labelGeometryImageFilter->GetWeightedCentroid(labelValue) << std::endl;
    std::cout << "\tAxes Length: " << labelGeometryImageFilter->GetAxesLength(labelValue) << std::endl;
    std::cout << "\tMajorAxisLength: " << labelGeometryImageFilter->GetMajorAxisLength(labelValue) << std::endl;
    std::cout << "\tMinorAxisLength: " << labelGeometryImageFilter->GetMinorAxisLength(labelValue) << std::endl;
    std::cout << "\tEccentricity: " << labelGeometryImageFilter->GetEccentricity(labelValue) << std::endl;
    std::cout << "\tElongation: " << labelGeometryImageFilter->GetElongation(labelValue) << std::endl;
    std::cout << "\tOrientation: " << labelGeometryImageFilter->GetOrientation(labelValue) << std::endl;
    std::cout << "\tBounding box: " << labelGeometryImageFilter->GetBoundingBox(labelValue) << std::endl;

    std::cout << std::endl << std::endl;
  }

  return EXIT_SUCCESS;
}

void
CreateIntensityImage(ImageType::Pointer image)
{
  // Create a random image.
  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size.Fill(20);

  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
  image->SetRegions(region);
  image->Allocate();

  itk::ImageRegionIterator<ImageType> imageIterator(image, image->GetLargestPossibleRegion());
  // Make a random image
  // Create an unchanging seed.
  srand(1);
  while (!imageIterator.IsAtEnd())
  {
    int randomNumber = rand() % 255;
    imageIterator.Set(randomNumber);
    ++imageIterator;
  }
}

void
CreateLabelImage(ImageType::Pointer image)
{
  // Create a black image with labeled squares.
  ImageType::IndexType start;
  start.Fill(0);

  ImageType::SizeType size;
  size.Fill(20);

  ImageType::RegionType region;
  region.SetSize(size);
  region.SetIndex(start);
  image->SetRegions(region);
  image->Allocate();

  itk::ImageRegionIterator<ImageType> imageIterator(image, image->GetLargestPossibleRegion());

  // Make a square
  while (!imageIterator.IsAtEnd())
  {
    if ((imageIterator.GetIndex()[0] > 5 && imageIterator.GetIndex()[0] < 10) &&
        (imageIterator.GetIndex()[1] > 5 && imageIterator.GetIndex()[1] < 10))
    {
      imageIterator.Set(85);
    }
    else if ((imageIterator.GetIndex()[0] > 11 && imageIterator.GetIndex()[0] < 17) &&
             (imageIterator.GetIndex()[1] > 11 && imageIterator.GetIndex()[1] < 17))
    {
      imageIterator.Set(127);
    }
    else if ((imageIterator.GetIndex()[0] > 11 && imageIterator.GetIndex()[0] < 17) &&
             (imageIterator.GetIndex()[1] > 1 && imageIterator.GetIndex()[1] < 5))
    {
      imageIterator.Set(191);
    }
    else
    {
      imageIterator.Set(0);
    }

    ++imageIterator;
  }
}

Classes demonstrated

template<typename TLabelImage, typename TIntensityImage = TLabelImage>
class LabelGeometryImageFilter : public itk::ImageToImageFilter<TLabelImage, TIntensityImage>

Given a label map and an optional intensity image, compute geometric features.

This filter enables the measurement of geometric features of all objects in a labeled ND volume. This labeled volume can represent, for instance, a medical image segmented into different anatomical structures or a microscope image segmented into individual cells. This filter is closely related to the itkLabelStatisticsImageFilter, which measures statistics of image regions defined by a labeled mask such as min, max, and mean intensity, intensity standard deviation, and bounding boxes. This filter, however, measures the geometry of the labeled regions themselves.

It calculates features similar to the regionprops command of Matlab. The set of measurements that it enables include: volume, centroid, eigenvalues, eigenvectors, axes lengths, eccentricity, elongation, orientation, bounding box, oriented bounding box, and rotation matrix. These features are based solely on the labeled mask itself. It also calculates integrated intensity and weighted centroid, which are measured on an intensity image under the labeled mask. These features represent the set of currently calculated features, but the framework of the filter is designed so that it can be easily expanded to measure a wide variety of other features.

This work is part of the National Alliance for Medical

Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149. Information on the National Centers for Biomedical Computing can be obtained from http://commonfund.nih.gov/bioinformatics.

Authors

Dirk Padfield and James Miller.

This filter was contributed in the Insight Journal paper: “A Label Geometry Image Filter for Multiple Object Measurement” by Padfield D., Miller J http://www.insight-journal.org/browse/publication/301

ITK Sphinx Examples:

• All ITK Sphinx Examples

• Geometric Properties Of Labeled Region

See itk::LabelGeometryImageFilter for additional documentation.