ITK  5.2.0
Insight Toolkit
SphinxExamples/src/Filtering/Convolution/NormalizedCorrelationOfMaskedImage/Code.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// Maximum value is inf??
#include "itkImage.h"
#include <iostream>
#include <string>
using ImageType = 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
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 =
CorrelationFilterType::Pointer correlationFilter = CorrelationFilterType::New();
correlationFilter->SetInput(image2);
correlationFilter->SetMaskImage(mask);
correlationFilter->SetTemplate(kernelOperator);
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");
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)
{
start.Fill(0);
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)
{
start.Fill(0);
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)
{
start.Fill(0);
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();
}
itk::ImageKernelOperator
A NeighborhoodOperator whose coefficients are from an image.
Definition: itkImageKernelOperator.h:50
itk::Index
Represent a n-dimensional index in a n-dimensional image.
Definition: itkIndex.h:66
itk::Size< 2 >
itk::MinimumMaximumImageCalculator
Computes the minimum and the maximum intensity values of an image.
Definition: itkMinimumMaximumImageCalculator.h:44
itkImageFileReader.h
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itkImage.h
itkMinimumMaximumImageCalculator.h
itk::Index::Fill
void Fill(IndexValueType value)
Definition: itkIndex.h:270
itk::Size::Fill
void Fill(SizeValueType value)
Definition: itkSize.h:211
itk::ImageRegionIterator
A multi-dimensional iterator templated over image type that walks a region of pixels.
Definition: itkImageRegionIterator.h:78
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itk::ImageFileWriter
Writes image data to a single file.
Definition: itkImageFileWriter.h:87
itk::Index::GetIndex
const IndexValueType * GetIndex() const
Definition: itkIndex.h:228
itk::GTest::TypedefsAndConstructors::Dimension2::RegionType
ImageBaseType::RegionType RegionType
Definition: itkGTestTypedefsAndConstructors.h:54
itkRescaleIntensityImageFilter.h
itkImageFileWriter.h
itk::ImageKernelOperator::SetImageKernel
void SetImageKernel(const ImageType *kernel)
itkRegionOfInterestImageFilter.h
itk::NormalizedCorrelationImageFilter
Computes the normalized correlation of an image and a template.
Definition: itkNormalizedCorrelationImageFilter.h:56
itk::RescaleIntensityImageFilter
Applies a linear transformation to the intensity levels of the input Image.
Definition: itkRescaleIntensityImageFilter.h:154
itkImageKernelOperator.h
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:86
itkNormalizedCorrelationImageFilter.h
itk::Size::GetSize
const SizeValueType * GetSize() const
Definition: itkSize.h:169