ITK/Examples/Broken/Images/MeanSquaresImageToImageMetric

This example extracts a patch from the input image and then attempts to find the best position of the patch in the input image (this process is clearly for demonstration only - we know the "ground truth" so we can compare it to the output).

The current output is (0,0) - this is incorrect. We expect either (50,50) (the bottom left corner of the extracted patch) or (75,75) (the center of the extracted patch).

MeanSquaresImageToImageMetric.cxx

<source lang="cpp">

1. include <itkImage.h>
2. include <itkTranslationTransform.h>
3. include <itkImageFileReader.h>
4. include <itkImageRegistrationMethod.h>
5. include <itkLinearInterpolateImageFunction.h>
6. include <itkRegularStepGradientDescentOptimizer.h>
7. include <itkMeanSquaresImageToImageMetric.h>
8. include <itkRegionOfInterestImageFilter.h>
9. include <itkCastImageFilter.h>
10. include <itkRescaleIntensityImageFilter.h>
11. include <itkImageRegionIterator.h>

1. include <itkImageToVTKImageFilter.h>

1. include "vtkImageViewer.h"
2. include "itkImageFileWriter.h"
3. include "vtkRenderWindowInteractor.h"
4. include "vtkSmartPointer.h"
5. include "vtkImageActor.h"
6. include "vtkInteractorStyleImage.h"
7. include "vtkRenderer.h"
8. include "vtkSphereSource.h"
9. include "vtkPolyDataMapper.h"
10. include "vtkActor.h"
11. include "vtkProperty.h"

1. include <iostream>
2. include <string>

typedef itk::Image<float, 2> FloatImageType;

int main(int argc, char *argv[]) {

 if(argc < 2)
   {
   std::cerr << "Required: filename" << std::endl;
   return EXIT_FAILURE;
   }
 
 std::string filename = argv[1];

 typedef itk::Image<unsigned char, 2> UnsignedCharImageType;
 
 typedef itk::ImageFileReader<FloatImageType> ReaderType;
 typedef itk::ImageToVTKImageFilter<UnsignedCharImageType> ConnectorType;

 // Read the image
 ReaderType::Pointer reader = ReaderType::New();
 reader->SetFileName(filename.c_str());
 reader->Update();

 // Setup the original image actor
 typedef itk::CastImageFilter< FloatImageType,
                               UnsignedCharImageType > CastFilterType;
 CastFilterType::Pointer originalCastFilter = CastFilterType::New();
 originalCastFilter->SetInput(reader->GetOutput());
 
 ConnectorType::Pointer originalConnector = ConnectorType::New();
 originalConnector->SetInput(originalCastFilter->GetOutput());

 vtkSmartPointer<vtkImageActor> originalActor =
   vtkSmartPointer<vtkImageActor>::New();
 originalActor->SetInput(originalConnector->GetOutput());

 // Extract a small region
 typedef itk::RegionOfInterestImageFilter< FloatImageType,
                                           FloatImageType > ExtractFilterType;

 ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();

 FloatImageType::IndexType start;
 start[0] = 50;
 start[1] = 50;

 FloatImageType::SizeType size;
 size[0] = 100;
 size[1] = 100;

 FloatImageType::RegionType desiredRegion;
 desiredRegion.SetSize(size);
 desiredRegion.SetIndex(start);

 extractFilter->SetRegionOfInterest(desiredRegion);
 extractFilter->SetInput(reader->GetOutput());
 extractFilter->Update();

 // Display extracted region (kernel)
 CastFilterType::Pointer kernelCastFilter = CastFilterType::New();
 kernelCastFilter->SetInput(extractFilter->GetOutput());
 
 ConnectorType::Pointer extractedConnector = ConnectorType::New();
 extractedConnector->SetInput(kernelCastFilter->GetOutput());

 vtkSmartPointer<vtkImageActor> extractedActor =
   vtkSmartPointer<vtkImageActor>::New();
 extractedActor->SetInput(extractedConnector->GetOutput());

 // Perform registration
 typedef itk::MeanSquaresImageToImageMetric<FloatImageType, FloatImageType> MetricType; // <TFixedImage, TMovingImage>
 typedef itk::TranslationTransform<double, 2> TransformType; // Compiler error if you change double to float??!
 typedef itk::RegularStepGradientDescentOptimizer OptimizerType;
 typedef itk::LinearInterpolateImageFunction<FloatImageType, double> InterpolatorType; // Compiler error if you change double to float??!
 typedef itk::ImageRegistrationMethod<FloatImageType, FloatImageType> RegistrationType;

 MetricType::Pointer metric = MetricType::New();
 TransformType::Pointer transform = TransformType::New();
 OptimizerType::Pointer optimizer = OptimizerType::New();
 InterpolatorType::Pointer interpolator = InterpolatorType::New();
 
 RegistrationType::Pointer registration = RegistrationType::New();
 registration->SetMetric(metric);
 registration->SetOptimizer(optimizer);
 registration->SetInterpolator(interpolator);
 registration->SetTransform(transform);
 registration->SetFixedImage(reader->GetOutput());
 registration->SetMovingImage(extractFilter->GetOutput());
 registration->SetFixedImageRegion(reader->GetOutput()->GetLargestPossibleRegion());

 RegistrationType::ParametersType initialParameters(transform->GetNumberOfParameters());
 initialParameters[0] = 0;
 initialParameters[1] = 0;

 registration->SetInitialTransformParameters(initialParameters);

 registration->Update();

 RegistrationType::ParametersType finalParameters = registration->GetLastTransformParameters();
 std::cout << "Final parameters: " << finalParameters << std::endl;

 vtkSmartPointer<vtkSphereSource> sphereSource =
   vtkSmartPointer<vtkSphereSource>::New();
 sphereSource->SetCenter(finalParameters[0], finalParameters[1], 0);
 sphereSource->SetRadius(10);
 sphereSource->Update();

 vtkSmartPointer<vtkPolyDataMapper> mapper =
   vtkSmartPointer<vtkPolyDataMapper>::New();
 mapper->SetInputConnection(sphereSource->GetOutputPort());

 vtkSmartPointer<vtkActor> sphereActor =
   vtkSmartPointer<vtkActor>::New();
 sphereActor->SetMapper(mapper);
 sphereActor->GetProperty()->SetColor(1,0,0);

 // Visualize
 vtkSmartPointer<vtkRenderWindow> renderWindow =
   vtkSmartPointer<vtkRenderWindow>::New();
 renderWindow->SetSize(900, 300);

 vtkSmartPointer<vtkRenderWindowInteractor> interactor =
   vtkSmartPointer<vtkRenderWindowInteractor>::New();
 interactor->SetRenderWindow(renderWindow);

 double leftViewport[4] = {0.0, 0.0, 0.5, 1.0};
 double rightViewport[4] = {0.5, 0.0, 1.0, 1.0};

 vtkSmartPointer<vtkRenderer> leftRenderer =
   vtkSmartPointer<vtkRenderer>::New();
 renderWindow->AddRenderer(leftRenderer);
 leftRenderer->SetViewport(leftViewport);
 leftRenderer->SetBackground(.6, .5, .4);

 vtkSmartPointer<vtkRenderer> rightRenderer =
   vtkSmartPointer<vtkRenderer>::New();
 renderWindow->AddRenderer(rightRenderer);
 rightRenderer->SetViewport(rightViewport);
 rightRenderer->SetBackground(.4, .5, .6);

 leftRenderer->AddActor(originalActor);
 leftRenderer->AddActor(sphereActor); // Display the location of the best position of the moving image in the fixed image
 
 rightRenderer->AddActor(extractedActor);
 
 leftRenderer->ResetCamera();
 rightRenderer->ResetCamera();

 renderWindow->Render();

 vtkSmartPointer<vtkInteractorStyleImage> style =
   vtkSmartPointer<vtkInteractorStyleImage>::New();
 interactor->SetInteractorStyle(style);

 interactor->Start();

 return EXIT_SUCCESS;

}

</source>

CMakeLists.txt

<source lang="cmake"> cmake_minimum_required(VERSION 2.6)

PROJECT(MeanSquaresImageToImageMetric)

include_directories(/home/doriad/src/ITK/Wrapping/WrapITK/ExternalProjects/ItkVtkGlue/src/)

1. include_directories(${ITK_SOURCE_DIR}/Wrapping/WrapITK/ExternalProjects/ItkVtkGlue/src/)

FIND_PACKAGE(VTK REQUIRED) INCLUDE(${VTK_USE_FILE})

FIND_PACKAGE(ITK REQUIRED) INCLUDE(${ITK_USE_FILE})

ADD_EXECUTABLE(MeanSquaresImageToImageMetric MeanSquaresImageToImageMetric.cxx) TARGET_LINK_LIBRARIES(MeanSquaresImageToImageMetric vtkHybrid ITKBasicFilters ITKCommon ITKIO ITKNumerics)

</source>