Normalized Correlation of Masked Image#

Synopsis#

Normalized correlation of a masked image.

Results#

mask.png

mask.png#

image1.png

image1.png#

image2.png

image2.png#

correlation.mha

correlation.mha#

correlation.png

correlation.png#

Output:

Maximum location fixed: [20, 6]
Maximum value: 1

Code#

C++#

// Maximum value is inf??

#include "itkImage.h"
#include "itkNormalizedCorrelationImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"
#include "itkImageKernelOperator.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkImageFileWriter.h"
#include "itkMinimumMaximumImageCalculator.h"

#include <iostream>
#include <string>

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

void
CreateMask(MaskType * const mask);
void
CreateImage(ImageType * const image);
void
CreateImageOfSquare(ImageType * const image, const itk::Index<2> & cornerOfSquare);

int
main()
{
  // Setup mask
  auto mask = MaskType::New();
  CreateMask(mask);

  itk::WriteImage(mask.GetPointer(), "mask.png");

  // Setup image1
  auto          image1 = ImageType::New();
  itk::Index<2> cornerOfSquare1;
  cornerOfSquare1[0] = 3;
  cornerOfSquare1[1] = 8;
  CreateImageOfSquare(image1, cornerOfSquare1);
  itk::WriteImage(image1.GetPointer(), "image1.png");

  // Setup image2
  itk::Index<2> offset;
  offset[0] = 20;
  offset[1] = 6;

  auto          image2 = ImageType::New();
  itk::Index<2> cornerOfSquare2;
  cornerOfSquare2[0] = cornerOfSquare1[0] + offset[0];
  cornerOfSquare2[1] = cornerOfSquare1[1] + offset[1];
  CreateImageOfSquare(image2, cornerOfSquare2);

  itk::WriteImage(image2.GetPointer(), "image2.png");

  // Create a kernel from an image
  itk::ImageKernelOperator<unsigned char> kernelOperator;
  kernelOperator.SetImageKernel(image1);

  // The radius of the kernel must be the radius of the patch, NOT the size of the patch
  itk::Size<2> radius;
  radius[0] = image1->GetLargestPossibleRegion().GetSize()[0] / 2;
  radius[1] = image1->GetLargestPossibleRegion().GetSize()[1] / 2;

  if (radius[0] % 2 == 0 || radius[1] % 2 == 0)
  {
    std::cerr << "Input must have odd dimensions!" << std::endl;
    return EXIT_FAILURE;
  }

  kernelOperator.CreateToRadius(radius);

  // Perform normalized correlation
  // <input type, mask type, output type>
  using CorrelationFilterType =
    itk::NormalizedCorrelationImageFilter<ImageType, MaskType, FloatImageType, unsigned char>;
  auto correlationFilter = CorrelationFilterType::New();
  correlationFilter->SetInput(image2);
  correlationFilter->SetMaskImage(mask);
  correlationFilter->SetTemplate(kernelOperator);
  correlationFilter->Update();

  itk::WriteImage(correlationFilter->GetOutput(), "correlation.mha");

  using RescaleFilterType = itk::RescaleIntensityImageFilter<FloatImageType, ImageType>;
  auto rescaleFilter = RescaleFilterType::New();
  rescaleFilter->SetInput(correlationFilter->GetOutput());
  rescaleFilter->SetOutputMinimum(0);
  rescaleFilter->SetOutputMaximum(255);
  rescaleFilter->Update();
  itk::WriteImage(rescaleFilter->GetOutput(), "correlation.png");

  using MinimumMaximumImageCalculatorType = itk::MinimumMaximumImageCalculator<FloatImageType>;
  MinimumMaximumImageCalculatorType::Pointer minimumMaximumImageCalculatorFilter =
    MinimumMaximumImageCalculatorType::New();
  minimumMaximumImageCalculatorFilter->SetImage(correlationFilter->GetOutput());
  minimumMaximumImageCalculatorFilter->Compute();

  itk::Index<2> maximumCorrelationPatchCenter = minimumMaximumImageCalculatorFilter->GetIndexOfMaximum();

  // This is the value we expect!
  std::cout << "Maximum location fixed: " << maximumCorrelationPatchCenter - radius << std::endl;

  // If the images can be perfectly aligned, the value is 1
  std::cout << "Maximum value: " << minimumMaximumImageCalculatorFilter->GetMaximum() << std::endl;

  return EXIT_SUCCESS;
}

void
CreateMask(MaskType * const mask)
{
  ImageType::IndexType start;
  start.Fill(0);

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

  ImageType::RegionType region(start, size);

  mask->SetRegions(region);
  mask->Allocate();
  mask->FillBuffer(255); // Make the whole mask "valid"

  // unsigned int squareSize = 8;
  unsigned int squareSize = 3;

  itk::Index<2> cornerOfSquare = { { 3, 8 } };

  // Remove pixels from the mask in a small square. The correlationw will not be computed at these pixels.
  itk::ImageRegionIterator<MaskType> maskIterator(mask, region);

  while (!maskIterator.IsAtEnd())
  {
    if (maskIterator.GetIndex()[0] > cornerOfSquare[0] && maskIterator.GetIndex()[0] < cornerOfSquare[0] + squareSize &&
        maskIterator.GetIndex()[1] > cornerOfSquare[1] && maskIterator.GetIndex()[1] < cornerOfSquare[1] + squareSize)
    {
      maskIterator.Set(0);
    }

    ++maskIterator;
  }
}

void
CreateImage(ImageType * const image)
{
  ImageType::IndexType start;
  start.Fill(0);

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

  ImageType::RegionType region(start, size);

  image->SetRegions(region);
  image->Allocate();
  image->FillBuffer(0);
}

void
CreateImageOfSquare(ImageType * const image, const itk::Index<2> & cornerOfSquare)
{
  ImageType::IndexType start;
  start.Fill(0);

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

  ImageType::RegionType region(start, size);

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

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

  unsigned int squareSize = 8;

  while (!imageIterator.IsAtEnd())
  {
    if (imageIterator.GetIndex()[0] > cornerOfSquare[0] &&
        imageIterator.GetIndex()[0] < cornerOfSquare[0] + squareSize &&
        imageIterator.GetIndex()[1] > cornerOfSquare[1] && imageIterator.GetIndex()[1] < cornerOfSquare[1] + squareSize)
    {
      imageIterator.Set(255);
    }

    ++imageIterator;
  }
}

Classes demonstrated#

template<typename TInputImage, typename TMaskImage, typename TOutputImage, typename TOperatorValueType = typename TOutputImage::PixelType>
class NormalizedCorrelationImageFilter : public itk::NeighborhoodOperatorImageFilter<TInputImage, TOutputImage, TOperatorValueType>

Computes the normalized correlation of an image and a template.

This filter calculates the normalized correlation between an image and the template. Normalized correlation is frequently use in feature detection because it is invariant to local changes in contrast.

The filter can be given a mask. When presented with an input image and a mask, the normalized correlation is only calculated at those pixels under the mask.

See

Image

See

Neighborhood

See

NeighborhoodOperator

See

NeighborhoodIterator

ITK Sphinx Examples:

See itk::NormalizedCorrelationImageFilter for additional documentation.