ITK  4.8.0
Insight Segmentation and Registration Toolkit
WikiExamples/Images/FFTNormalizedCorrelationImageFilter.cxx
#include "itkImage.h"
#include <iostream>
#include <string>
typedef itk::Image<unsigned char, 2> ImageType;
typedef itk::Image<float, 2> FloatImageType;
//typedef itk::Image<unsigned char, 2> UnsignedCharImageType;
static void CreateImage(ImageType::Pointer image, const itk::Index<2>& cornerOfSquare);
template <typename TImage>
static void WriteImage(TImage* const image, const std::string& filename);
int main(int argc, char *argv[])
{
itk::Index<2> offset;
offset[0] = 5;
offset[1] = 6;
ImageType::Pointer fixedImage = ImageType::New();
itk::Index<2> cornerOfFixedSquare;
cornerOfFixedSquare[0] = 3;
cornerOfFixedSquare[1] = 8;
CreateImage(fixedImage, cornerOfFixedSquare);
WriteImage(fixedImage.GetPointer(), "fixedImage.png");
ImageType::Pointer movingImage = ImageType::New();
itk::Index<2> cornerOfMovingSquare;
cornerOfMovingSquare[0] = cornerOfFixedSquare[0] + offset[0];
cornerOfMovingSquare[1] = cornerOfFixedSquare[1] + offset[1];
CreateImage(movingImage, cornerOfMovingSquare);
WriteImage(movingImage.GetPointer(), "movingImage.png");
// Perform normalized correlation
CorrelationFilterType::Pointer correlationFilter = CorrelationFilterType::New();
correlationFilter->SetFixedImage(fixedImage);
correlationFilter->SetMovingImage(movingImage);
correlationFilter->Update();
WriteImage(correlationFilter->GetOutput(), "correlation.mha");
RescaleFilterType::Pointer rescaleFilter = RescaleFilterType::New();
rescaleFilter->SetInput(correlationFilter->GetOutput());
rescaleFilter->SetOutputMinimum(0);
rescaleFilter->SetOutputMaximum(255);
rescaleFilter->Update();
WriteImage(rescaleFilter->GetOutput(), "correlation.png");
typedef itk::MinimumMaximumImageCalculator<FloatImageType> MinimumMaximumImageCalculatorType;
MinimumMaximumImageCalculatorType::Pointer minimumMaximumImageCalculatorFilter = MinimumMaximumImageCalculatorType::New ();
minimumMaximumImageCalculatorFilter->SetImage(correlationFilter->GetOutput());
minimumMaximumImageCalculatorFilter->Compute();
itk::Index<2> maximumCorrelationPatchCenter = minimumMaximumImageCalculatorFilter->GetIndexOfMaximum();
itk::Size<2> outputSize = correlationFilter->GetOutput()->GetLargestPossibleRegion().GetSize();
itk::Index<2> maximumCorrelationPatchCenterFixed;
maximumCorrelationPatchCenterFixed[0] = outputSize[0]/2 - maximumCorrelationPatchCenter[0];
maximumCorrelationPatchCenterFixed[1] = outputSize[1]/2 - maximumCorrelationPatchCenter[1];
std::cout << "Maximum location: " << maximumCorrelationPatchCenter << std::endl;
std::cout << "Maximum location fixed: " << maximumCorrelationPatchCenterFixed << std::endl; // This is the value we expect!
std::cout << "Maximum value: " << minimumMaximumImageCalculatorFilter->GetMaximum() << std::endl; // If the images can be perfectly aligned, the value is 1
return EXIT_SUCCESS;
}
void CreateImage(ImageType::Pointer 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);
ImageType::IndexValueType 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(TImage* const image, const std::string& filename)
{
typedef itk::ImageFileWriter<TImage> WriterType;
typename WriterType::Pointer writer = WriterType::New();
writer->SetFileName(filename);
writer->SetInput(image);
writer->Update();
}