Geometric Properties of Labeled Region#

Synopsis#

Get geometric properties of labeled regions in an image.

Results#

Output:

Number of labels: 3

      Label: 85
      Physical size: 16
      Number of pixels: 16
      Centroid: [7.5, 7.5]
      Equivalent ellipsoid diameter: [4.51352, 4.51352]
      Elongation: 1
      Flatness: 1
      Principal axes: 1 0
0 1

      Bounding box: ImageRegion
  Dimension: 2
  Index: [6, 6]
  Size: [4, 4]

      Integrated intensity: 2487
      Mean intensity: 155.438


      Label: 127
      Physical size: 25
      Number of pixels: 25
      Centroid: [14, 14]
      Equivalent ellipsoid diameter: [5.6419, 5.6419]
      Elongation: 1
      Flatness: 1
      Principal axes: 1 0
0 1

      Bounding box: ImageRegion
  Dimension: 2
  Index: [12, 12]
  Size: [5, 5]

      Integrated intensity: 4310
      Mean intensity: 172.4


      Label: 191
      Physical size: 15
      Number of pixels: 15
      Centroid: [14, 3]
      Equivalent ellipsoid diameter: [3.32063, 5.7515]
      Elongation: 1.73205
      Flatness: 1.73205
      Principal axes: 0 1
-1 -0

      Bounding box: ImageRegion
  Dimension: 2
  Index: [12, 2]
  Size: [5, 3]

      Integrated intensity: 2155
      Mean intensity: 143.667

Code#

C++#

#include "itkImage.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkImageFileReader.h"
#include "itkLabelImageToShapeLabelMapFilter.h"
#include "itkLabelStatisticsImageFilter.h"
#include "itkLabelToRGBImageFilter.h"

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

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();

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

  // Shape attributes (size, centroid, ellipsoid axes, elongation, ...).
  using LabelImageToShapeLabelMapFilterType = itk::LabelImageToShapeLabelMapFilter<ImageType>;
  auto shapeFilter = LabelImageToShapeLabelMapFilterType::New();
  shapeFilter->SetInput(labelImage);
  shapeFilter->Update();

  // Intensity statistics computed over each labeled region.
  using LabelStatisticsImageFilterType = itk::LabelStatisticsImageFilter<ImageType, ImageType>;
  auto statisticsFilter = LabelStatisticsImageFilterType::New();
  statisticsFilter->SetInput(intensityImage);
  statisticsFilter->SetLabelInput(labelImage);
  statisticsFilter->Update();

  auto labelMap = shapeFilter->GetOutput();
  std::cout << "Number of labels: " << labelMap->GetNumberOfLabelObjects() << std::endl;
  std::cout << std::endl;

  for (unsigned int n = 0; n < labelMap->GetNumberOfLabelObjects(); ++n)
  {
    const auto labelObject = labelMap->GetNthLabelObject(n);
    const auto labelValue = labelObject->GetLabel();

    std::cout << "\tLabel: " << static_cast<int>(labelValue) << std::endl;
    std::cout << "\tPhysical size: " << labelObject->GetPhysicalSize() << std::endl;
    std::cout << "\tNumber of pixels: " << labelObject->GetNumberOfPixels() << std::endl;
    std::cout << "\tCentroid: " << labelObject->GetCentroid() << std::endl;
    std::cout << "\tEquivalent ellipsoid diameter: " << labelObject->GetEquivalentEllipsoidDiameter() << std::endl;
    std::cout << "\tElongation: " << labelObject->GetElongation() << std::endl;
    std::cout << "\tFlatness: " << labelObject->GetFlatness() << std::endl;
    std::cout << "\tPrincipal axes: " << labelObject->GetPrincipalAxes() << std::endl;
    std::cout << "\tBounding box: " << labelObject->GetBoundingBox() << std::endl;

    if (statisticsFilter->HasLabel(labelValue))
    {
      std::cout << "\tIntegrated intensity: " << statisticsFilter->GetSum(labelValue) << std::endl;
      std::cout << "\tMean intensity: " << statisticsFilter->GetMean(labelValue) << std::endl;
    }

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

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

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

  return EXIT_SUCCESS;
}

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

  auto size = ImageType::SizeType::Filled(20);

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

  itk::ImageRegionIteratorWithIndex<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{};

  auto size = ImageType::SizeType::Filled(20);

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

  itk::ImageRegionIteratorWithIndex<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 TInputImage, typename TOutputImage = LabelMap<ShapeLabelObject<typename TInputImage::PixelType, TInputImage::ImageDimension>>> class LabelImageToShapeLabelMapFilter : public ImageToImageFilter<TInputImage, LabelMap<ShapeLabelObject<typename TInputImage::PixelType, TInputImage::ImageDimension>>>

Converts a label image to a label map and valuates the shape attributes.

A convenient class that converts a label image to a label map and valuates the shape attribute at once.

This implementation was taken from the Insight Journal paper: https://doi.org/10.54294/q6auw4