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 "itkImageFileReader.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);

template <typename TImage>
void
WriteImage(const TImage * const image, const std::string & filename);

int
main(int, char *[])
{
  // Setup mask
  MaskType::Pointer mask = MaskType::New();
  CreateMask(mask);

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

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

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

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

  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>;
  CorrelationFilterType::Pointer correlationFilter = CorrelationFilterType::New();
  correlationFilter->SetInput(image2);
  correlationFilter->SetMaskImage(mask);
  correlationFilter->SetTemplate(kernelOperator);
  correlationFilter->Update();

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

  using RescaleFilterType = itk::RescaleIntensityImageFilter<FloatImageType, ImageType>;
  RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
  rescaleFilter->SetInput(correlationFilter->GetOutput());
  rescaleFilter->SetOutputMinimum(0);
  rescaleFilter->SetOutputMaximum(255);
  rescaleFilter->Update();
  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;
  }
}

template <typename TImage>
void
WriteImage(const TImage * const image, const std::string & filename)
{
  using WriterType = itk::ImageFileWriter<TImage>;
  typename WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(filename);
  writer->SetInput(image);
  writer->Update();
}

Classes demonstrated