Download this example

K Means Cluster of Pixels in Image

Warning

Fix Problem Contains problems not fixed from original wiki.

Synopsis

Compute kmeans clusters of pixels in an image

Results

Note

Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>

Code

C++

#include "itkImage.h"
#include "itkListSample.h"
#include "itkVector.h"
#include "itkImageKmeansModelEstimator.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkImageToListSampleAdaptor.h"
#include "itkDistanceToCentroidMembershipFunction.h"
#include "itkSampleClassifierFilter.h"
#include "itkMinimumDecisionRule.h"
#include "itkImageFileWriter.h"

using MeasurementVectorType = itk::Vector<unsigned char, 3>;
using ColorImageType = itk::Image<MeasurementVectorType, 2>;
using ScalarImageType = itk::Image<unsigned char, 2>;

static void
CreateImage(ColorImageType::Pointer image);

int
main()
{
  // Create a demo image
  auto image = ColorImageType::New();
  CreateImage(image);

  // Compute pixel clusters using KMeans
  using MembershipFunctionType = itk::Statistics::DistanceToCentroidMembershipFunction<MeasurementVectorType>;
  using MembershipFunctionPointer = MembershipFunctionType::Pointer;
  using MembershipFunctionPointerVector = std::vector<MembershipFunctionPointer>;

  using ImageKmeansModelEstimatorType = itk::ImageKmeansModelEstimator<ColorImageType, MembershipFunctionType>;

  auto kmeansEstimator = ImageKmeansModelEstimatorType::New();
  kmeansEstimator->SetInputImage(image);
  kmeansEstimator->SetNumberOfModels(3);
  kmeansEstimator->SetThreshold(0.01);
  kmeansEstimator->SetOffsetAdd(0.01);
  kmeansEstimator->SetOffsetMultiply(0.01);
  kmeansEstimator->SetMaxSplitAttempts(10);
  kmeansEstimator->Update();

  // Classify each pixel
  using SampleType = itk::Statistics::ListSample<MeasurementVectorType>;
  using ClassifierType = itk::Statistics::SampleClassifierFilter<SampleType>;
  auto classifier = ClassifierType::New();

  using DecisionRuleType = itk::Statistics::MinimumDecisionRule;
  auto decisionRule = DecisionRuleType::New();

  classifier->SetDecisionRule(decisionRule);
  classifier->SetNumberOfClasses(3);

  using ClassLabelVectorObjectType = ClassifierType::ClassLabelVectorObjectType;
  using ClassLabelVectorType = ClassifierType::ClassLabelVectorType;
  using MembershipFunctionVectorObjectType = ClassifierType::MembershipFunctionVectorObjectType;
  using MembershipFunctionVectorType = ClassifierType::MembershipFunctionVectorType;

  // Setup membership functions
  MembershipFunctionPointerVector kmeansMembershipFunctions = kmeansEstimator->GetMembershipFunctions();

  MembershipFunctionVectorObjectType::Pointer membershipFunctionsVectorObject =
    MembershipFunctionVectorObjectType::New();
  classifier->SetMembershipFunctions(membershipFunctionsVectorObject);

  MembershipFunctionVectorType & membershipFunctionsVector = membershipFunctionsVectorObject->Get();

  for (auto & kmeansMembershipFunction : kmeansMembershipFunctions)
  {
    membershipFunctionsVector.push_back(kmeansMembershipFunction.GetPointer());
  }

  // Setup class labels
  auto classLabelsObject = ClassLabelVectorObjectType::New();
  classifier->SetClassLabels(classLabelsObject);

  ClassLabelVectorType & classLabelsVector = classLabelsObject->Get();
  classLabelsVector.push_back(50);
  classLabelsVector.push_back(150);
  classLabelsVector.push_back(250);

  // Perform the classification
  using SampleAdaptorType = itk::Statistics::ImageToListSampleAdaptor<ColorImageType>;
  auto sample = SampleAdaptorType::New();
  sample->SetImage(image);

  classifier->SetInput(sample);
  classifier->Update();

  // Prepare the output image
  auto outputImage = ScalarImageType::New();
  outputImage->SetRegions(image->GetLargestPossibleRegion());
  outputImage->Allocate();
  outputImage->FillBuffer(0);

  // Setup the membership iterator
  const ClassifierType::MembershipSampleType *        membershipSample = classifier->GetOutput();
  ClassifierType::MembershipSampleType::ConstIterator membershipIterator = membershipSample->Begin();

  // Setup the output image iterator - this is automatically synchronized with the membership iterator since the sample
  // is an adaptor
  itk::ImageRegionIteratorWithIndex<ScalarImageType> outputIterator(outputImage,
                                                                    outputImage->GetLargestPossibleRegion());
  outputIterator.GoToBegin();

  while (membershipIterator != membershipSample->End())
  {
    int classLabel = membershipIterator.GetClassLabel();
    // std::cout << "Class label: " << classLabel << std::endl;
    outputIterator.Set(classLabel);
    ++membershipIterator;
    ++outputIterator;
  }

  itk::WriteImage(image, "input.mha");
  itk::WriteImage(outputImage, "output.mha");


  return EXIT_SUCCESS;
}

void
CreateImage(ColorImageType::Pointer image)
{
  // Create a black image with a red square and a green square
  ColorImageType::RegionType region;
  ColorImageType::IndexType  start;
  start[0] = 0;
  start[1] = 0;

  ColorImageType::SizeType size;
  size[0] = 200;
  size[1] = 300;

  region.SetSize(size);
  region.SetIndex(start);

  image->SetRegions(region);
  image->Allocate();

  itk::ImageRegionIterator<ColorImageType> imageIterator(image, region);

  itk::Vector<unsigned char, 3> redPixel;
  redPixel[0] = 255;
  redPixel[1] = 0;
  redPixel[2] = 0;

  itk::Vector<unsigned char, 3> greenPixel;
  greenPixel[0] = 0;
  greenPixel[1] = 255;
  greenPixel[2] = 0;

  itk::Vector<unsigned char, 3> blackPixel;
  blackPixel[0] = 0;
  blackPixel[1] = 0;
  blackPixel[2] = 0;

  while (!imageIterator.IsAtEnd())
  {
    if (imageIterator.GetIndex()[0] > 100 && imageIterator.GetIndex()[0] < 150 && imageIterator.GetIndex()[1] > 100 &&
        imageIterator.GetIndex()[1] < 150)
    {
      imageIterator.Set(redPixel);
    }
    else if (imageIterator.GetIndex()[0] > 50 && imageIterator.GetIndex()[0] < 70 && imageIterator.GetIndex()[1] > 50 &&
             imageIterator.GetIndex()[1] < 70)
    {
      imageIterator.Set(greenPixel);
    }
    else
    {
      imageIterator.Set(blackPixel);
    }

    ++imageIterator;
  }
}

Classes demonstrated

template<typename TInputImage, typename TMembershipFunction>
class ImageKmeansModelEstimator : public itk::ImageModelEstimatorBase<TInputImage, TMembershipFunction>

Base class for ImageKmeansModelEstimator object.

itkImageKmeansModelEstimator generates the kmeans model (cluster centers). This object performs clustering of data sets into different clusters, either using user-provided seed points as an initial guess or generating the clusters using a recursive approach when the user provides the number of desired clusters. Each cluster is represented by its cluster center. The two algorithms used are the generalized Lloyd algorithm (GLA) and the Linde-Buzo-Gray algorithm. The cluster centers are also referred to as codewords and a table of cluster centers is referred as a codebook.

As required by the GLA algorithm, the initial seed cluster should contain approximate centers of clusters. The GLA algorithm genrates updated cluster centers that result in a lower distortion than the input seed cluster when the input vectors are mapped/classified/labelled using the given codebooks.

If no codebook is provided, the Linde-Buzo-Gray algorithm is used. This algorithm uses the GLA algorithm at its core to generate the centroids of the input vectors (data). However, since there is no initial codebook, LBG first creates a one word codebook (or centroid of one cluster comprising of all the input training vectors). The LBG uses codeword or centroid splitting to create an increasing number of clusters. Each new set of clusters are optimized using the GLA algorithm. The number of clusters increases as $2^{n}$ n= 0, 1, … The codebook is expected to be in the form of a vnl matrix, where there are N rows, each row representing the cluster mean of a given cluster. The number of columns in the codebook should be equal to the input image vector dimension.

The threshold parameter controls the ‘’optimality’’ of the returned codebook, where optimality is related to the least possible mean-squared error distortion that can be found by the algorithm. For larger thresholds, the result will be less optimal. For smaller thresholds, the result will be more optimal. If a more optimal result is desired, then the algorithm will take longer to complete. A reasonable threshold value is 0.01.

If, during the operation of the algorithm, there are any unused clusters or cells, the m_OffsetAdd and m_OffsetMultiply parameters are used to split the cells with the highest distortion. This function will attempt to fill empty cells up to 10 times (unless the overall distortion is zero). Using 0.01 is a reasonable default values for the m_OffsetAdd and m_OffsetMultiply parameters.

If the GLA is unable to resolve the data into the desired number of clusters or cells, only the codewords which were used will be returned.

In terms of clustering, codewords are cluster centers, and a codebook is a table containing all cluster centers. The GLA produces results that are equivalent to the K-means clustering algorithm.

For more information about the algorithms, see A. Gersho and R. M. Gray, {Vector Quantization and Signal Compression}, Kluwer Academic Publishers, Boston, MA, 1992.

This object supports data handling of multiband images. The object accepts the input image in vector format only, where each pixel is a vector and each element of the vector corresponds to an entry from 1 particular band of a multiband dataset. A single band image is treated as a vector image with a single element for every vector.

This function is templated over the type of input image. In addition, a second parameter for the MembershipFunction needs to be specified. In this case a Membership function that store cluster centroids models needs to be specified.

The Update() function enables the calculation of the various models, creates the membership function objects and populates them.

Note: There is a second implementation of k-means algorithm in ITK under the itk::Statistics namespace. While this algorithm (GLA/LBG based algorithm) is memory efficient, the other algorithm is time efficient.

See

KdTreeBasedKmeansEstimator, WeightedCentroidKdTreeGenerator, KdTree

See

ScalarImageKmeansImageFilter

ITK Sphinx Examples:

• All ITK Sphinx Examples

• K Means Cluster Of Pixels In Image

See itk::ImageKmeansModelEstimator for additional documentation.