[Insight-users] problems computing similarity between images

[Luis Ibanez]

This is a multi-part message in MIME format.
--------------020309040903080602040700
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit


Hi Petra,

Thanks for posting the detailed code of your test.
After running it on debugging here is what I can see:

The ExtractImageFilter is not appropriate for the
purpose of this test, given that the image you get
as output has its regions set *exactly* as the
"extraction region".

That is, if you requested to extract a region with
start index (15,27) and size (300,400), form a full
image of region (0,0)-O(1000,1000), the resulting
image at the output will have regions (buffered,
requested and largest) equal to

              (15,27)-(1000,1000),

and the origin set to the same origin of your input
image.  When this image maps indices to coordinates
in physical space, its mappings are placed exactly
in the same location as the original image.

In other words, the extracted region is placed in
the same spatial location where it was when it used
to belong to the larger image.

You could use the RegionOfInterest image filter,
but beware that this one will put the region to
(0,0)-(300,400) and the origin to 15*spacingX,
27*spacingY. So, it will also be mapped to the
same physical location.    :-/


-------


What I would suggest you is to use the RegionOfInterest
ImageFilter AND use your initial transform for specifying
that you want to move that extracted region to be on top
of the other region (the fixed image region).

Let's say that you want to register two regions of size
(300,400) and with starting indices:

Region1 start index = ( 10, 20)
Region2 start index = (120,230)

You will use Region1 to set the FixedImageRegion() on the
registration method, and Region2 for the RegionOfInterest
ImageFilter.

Then, prepare the initial transform to be the one
that moves Region2 on top of Region1. That will be

   TranslationParam[0] = ( 120 - 10 ) * spacing[0];
   TranslationParam[1] = ( 230 - 30 ) * spacing[1];

and set this parameters to the metric...

Then you will actually be comparing Region2 to Region1.



Please find attached the modified code that do this
preparatoin on the initial transform, as well as its
respective CMakeLists.txt file.

The program is expecting now the starting index of
region2 as command line argument.

The final word about the test would be to use the output
transform for resampling region2 on the coordinate system
of the fixed image.

You will find several examples on how  to do this on the
directory:

        Insight/Examples/Registration

You will find details about these examples in the
software guide

     http://www . itk . org/ItkSoftwareGuide . pdf

- ExtractImageFilter in Section 7.5
- RegionOfInterestImageFilter in Section 7.4
- ResampleImageFilter in Section 6.7.1
- RegistrationMethods in Chapter 8



  Please let us know if you have further questions.



    Thanks




----------------------
Petra Schneider wrote:
> Hello,
> 
> my intention is very simple. I need to compute the similarity between different
> regions of two images using NormalizedCorrelationImageToImageMetric (see my
> code below).
> But there seems to be a mistake. In some cases, I get the optimal metric value,
> if I compare different regions of the same image. What's wrong? Need help.....
> 
> 
> Thanks,
> 	Petra
> 
> 

--------------020309040903080602040700
Content-Type: text/plain;
 name="NormalizedCorrelation.cxx"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="NormalizedCorrelation.cxx"

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkNormalizedCorrelationImageToImageMetric.h"
#include "itkImageRegistrationMethod.h"
#include "itkTranslationTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkRegularStepGradientDescentOptimizer.h"
#include "itkExtractImageFilter.h"
#include "itkRegionOfInterestImageFilter.h"


int main( int argc, char ** args )
{

  if( argc < 5 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << "NormalizedCorrelation   imagefileFixed imagefileMoving  startX startY" << std::endl;
    return -1;
    }

	typedef unsigned char  pixel_value;
	static const unsigned int  image_dimension ( 2 );
	typedef itk :: Image < pixel_value, image_dimension > 	 image_type;

	typedef itk :: ImageFileReader < image_type >	 reader_type;
	reader_type :: Pointer	reader_fixed ( reader_type :: New() );
	reader_type :: Pointer	reader_moving ( reader_type :: New() );


	reader_fixed -> SetFileName ( args[ 1 ] );
	image_type :: Pointer image_fixed = reader_fixed -> GetOutput();
	image_fixed -> Update();

  reader_moving -> SetFileName ( args[ 2 ] );
	image_type :: Pointer image_moving = reader_moving -> GetOutput();
	image_moving -> Update();


	typedef itk :: TranslationTransform < double, image_dimension > transform_type;
	typedef itk :: LinearInterpolateImageFunction < image_type, double > interpolator_type;
	typedef itk :: NormalizedCorrelationImageToImageMetric < image_type, image_type >  metric_type;
	typedef itk :: RegionOfInterestImageFilter < image_type, image_type > extractor_type;


	transform_type :: Pointer transform ( transform_type :: New() );
	interpolator_type :: Pointer interpolator ( interpolator_type :: New());
	metric_type :: Pointer metric ( metric_type :: New() );
	extractor_type :: Pointer extractor ( extractor_type :: New() );




	image_type :: RegionType region_moving;
	image_type :: RegionType region_fixed;
	image_type :: IndexType index_moving;
	image_type :: IndexType index_fixed;
	image_type :: SizeType size_moving;
	image_type :: SizeType size_fixed;


	index_moving[ 0 ] = atoi( args[3] );		
  index_moving[ 1 ] = atoi( args[4] );
	
  index_fixed[ 0 ]  = 0;		index_fixed[  1 ] = 0;
	size_moving[ 0 ] = 50;		size_moving[  1 ] = 55;
	size_fixed[ 0 ]  = 50;		size_fixed[   1 ] = 55;


	region_moving.SetSize( size_moving );
	region_moving.SetIndex( index_moving );
	region_fixed.SetSize ( size_fixed );
	region_fixed.SetIndex ( index_fixed );


	extractor -> SetInput ( image_moving );
	extractor -> SetRegionOfInterest( region_moving );
	extractor -> Update();
	const image_type :: Pointer moving ( extractor -> GetOutput() );

  std::cout << "Moving image is" << std::endl;
  moving->Print( std::cout );

  std::cout << "Fixed image is" << std::endl;
  image_fixed->Print( std::cout );

  const double * spacing = image_fixed->GetSpacing();

	typedef transform_type :: ParametersType  parameters_type;
	parameters_type parameters ( transform -> GetNumberOfParameters() );

  parameters[0] = ( index_moving[0] - index_fixed[0] ) * spacing[0];
  parameters[1] = ( index_moving[1] - index_fixed[1] ) * spacing[1];

	metric -> SetFixedImage ( image_fixed );
	metric -> SetFixedImageRegion( image_fixed -> GetRequestedRegion() );
	metric -> SetMovingImage( moving );
	metric -> SetTransform ( transform );
	metric -> SetTransformParameters ( parameters );
	metric -> SetInterpolator ( interpolator );
	metric -> Initialize();


  std::cout << metric -> GetValue ( parameters ) << std::endl;

}


--------------020309040903080602040700
Content-Type: text/plain;
 name="CMakeLists.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
 filename="CMakeLists.txt"

PROJECT(myProject)

FIND_PACKAGE(ITK)
IF(ITK_FOUND)
	INCLUDE(${ITK_USE_FILE})
ELSE(ITK_FOUND)
	MESSAGE(FATAL_ERROR "ITK not found. Please set ITK_DIR.")
ENDIF(ITK_FOUND)

ADD_EXECUTABLE(myProject NormalizedCorrelation.cxx )

TARGET_LINK_LIBRARIES(myProject ITKIO ITKCommon ITKNumerics)


--------------020309040903080602040700--