[Insight-users] mean square metric
Fab_83
Fabian_Eisa at gmx.de
Wed Apr 9 04:33:15 EDT 2008
Hi Louis,
thanks for your answer. I see why the pixels outside the domain shouldn't be
considered. I changed my code but still I don't get the same values than
itk. Enclose I posted my new code. The function SumOfAbsolutDistance (I
know actually the name should be SumOfSquareDifferenz but after the
validation I want to calculate the Sum of absolut distance with this
function) should do exactly what itk is doing in the metric
itk::MeanSquaresImageToImageMetric. Itk is calculating the metric for every
quatered voxel translation and my function for whole voxel translation
steps. That means that there should be 10 values that are the same (every
4th value in itk should be the same than the function value).
I spend plenty of time to find the problem but I am puzzled.
Thanks a lot,
Fab
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkTranslationTransform.h"
#include "itkNearestNeighborInterpolateImageFunction.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkResampleImageFilter.h"
#include "itkImageRegionIteratorWithIndex.h"
#include "itkMeanSquaresImageToImageMetric.h"
#include "itkNormalizedCorrelationImageToImageMetric.h"
#include "itkMeanReciprocalSquareDifferenceImageToImageMetric.h"
#include "itkMutualInformationImageToImageMetric.h"
#include "itkMattesMutualInformationImageToImageMetric.h"
#include "itkNormalizedMutualInformationHistogramImageToImageMetric.h"
#include "itkMeanSquaresHistogramImageToImageMetric.h"
#include "itkCorrelationCoefficientHistogramImageToImageMetric.h"
#include "itkMatchCardinalityImageToImageMetric.h"
#include "itkKappaStatisticImageToImageMetric.h"
#include "itkGradientDifferenceImageToImageMetric.h"
#include <iostream>
using namespace std;
float faktor=1;
const int rangex = 10*faktor;
const int rangey = 0;
const int rangez = 0;
const unsigned int Dimension = 3;
typedef short PixelType;
typedef itk::Image< PixelType, Dimension > ImageType;
typedef itk::TranslationTransform< double, Dimension > TransformType;
TransformType::Pointer transform = TransformType::New();
typedef itk::LinearInterpolateImageFunction< ImageType, double >
InterpolatorType;
InterpolatorType::Pointer interpolator = InterpolatorType::New();
typedef itk::ImageRegionIteratorWithIndex<ImageType> IterType;
typedef itk::ResampleImageFilter<ImageType, ImageType > FilterType;
FilterType::Pointer filter = FilterType::New();
ImageType::Pointer fixedImage;
ImageType::Pointer movingImage;
template <class MetricType>
void generatefile(typename MetricType::Pointer metric, string filename)
{
double result=0;
//double trans=1;
double hilf=0;
double previousvalue=0, relVal, Steigung=0;
double metricvalue;
double worstcase, bestcase, relworstcase;
metric->SetTransform( ::transform ); //:: nur lokal
metric->SetInterpolator( interpolator );
metric->SetFixedImage( fixedImage );
metric->SetMovingImage( movingImage );
metric->SetFixedImageRegion( fixedImage->GetBufferedRegion() );
try
{
metric->Initialize();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception catched: here is something wrong with the memory
!" << std::endl;
std::cerr << excep << std::endl;
exit(1);
}
MetricType::TransformParametersType displacement( Dimension );
ofstream output(filename.c_str());
displacement[0]=0;
displacement[1]=0;
displacement[2]=0;
bestcase=metric->GetValue(displacement);
displacement[0]=rangex/faktor;
worstcase=metric->GetValue(displacement);
relworstcase=worstcase-bestcase;
output<<"original Normiert"<<std::endl;
float frac;
for( float dx = 0; dx <= rangex; dx++ )
{
for( int dy = -rangey; dy <= rangey; dy++ )
{
for (int dz = -rangez; dz<=rangez; dz++)
{
frac=dx/faktor;
displacement[0]= frac;
displacement[1]=dy;
displacement[2]=dz;
metricvalue = metric->GetValue( displacement );
result=(metricvalue-previousvalue)/metricvalue;
previousvalue=metricvalue;
relVal=(metricvalue-bestcase)/relworstcase;
hilf=relVal;
output <<metricvalue<<" "<<relVal<<" "<<endl;
}
}
}
}
void SumOfAbsolutDistance ()
{
typedef itk::NumericTraits<PixelType>::RealType PixelRealType;
typedef itk::ImageFileWriter< ImageType > WriterType;
WriterType::Pointer writer = WriterType::New();
writer->SetFileName( "Test.mhd" );
TransformType::Pointer transform = TransformType::New();
float fraction=1;
double sum=0;
double N=0;
ImageType::Pointer Image;
ImageType::Pointer ResultImage = ImageType::New();
//////////////////
int layers=100;
int size=130;
ImageType::IndexType start;
start.Fill(0);
ImageType::SizeType dim; //dimension
dim[0]=size; //size along X
dim[1]=size; //size along Y
dim[2]=layers; //size along Z
ImageType::SpacingType spacing;
spacing.Fill(1);
ResultImage->SetSpacing( spacing );
ImageType::PointType origin;
origin.Fill(0);
ResultImage->SetOrigin( origin );
ImageType::RegionType region; // creates the region
region.SetSize(dim);
region.SetIndex(start);
ResultImage->SetRegions(region); // the region is passed to the object
"image"
ResultImage->Allocate();
/////////////////////////
typedef itk::ImageRegionIteratorWithIndex<ImageType> IterType;
ImageType::RegionType reg = ResultImage->GetRequestedRegion();
ImageType::PixelType FixedPixelValue;
ImageType::PixelType MovingPixelValue;
ImageType::PixelType ResultPixel;
filter->SetInterpolator(interpolator);
filter->SetDefaultPixelValue( 0 );
filter->SetOutputParametersFromImage(movingImage);
filter->SetInput( movingImage );
ImageType::SpacingType spc=movingImage->GetSpacing();
TransformType::ParametersType translation = transform->GetParameters();
ofstream outputFile;
outputFile.open("SSD_SelfCalculated.txt");
ImageType::SizeType sz = reg.GetSize();
for (int i=0; i<10; i++)
{
translation[0] = fraction*spc[0]*i;
translation[1] = spc[1]*0;
translation[2] = spc[2]*0;
transform->SetParameters( translation );
filter->SetTransform( transform );
Image=filter->GetOutput();
filter->Update();
--sz[0];
reg.SetSize(sz);
IterType it(ResultImage, reg);
sum = 0;
N = 0;
for (it.GoToBegin(); !it.IsAtEnd(); ++it)
{
IterType::IndexType ind = it.GetIndex();
FixedPixelValue=fixedImage->GetPixel(ind);
MovingPixelValue=Image->GetPixel(ind);
double _diff=FixedPixelValue-MovingPixelValue;
ResultPixel = static_cast<PixelType>(_diff+0.5);
it.Set(ResultPixel);
////////////////////////////// SSD ////////////////////////
sum += _diff * _diff;
N++;
}
outputFile << sum/N << std::endl;
}
outputFile.close();
}
int main( int argc, char * argv[] )
{
if( argc < 3 )
{
std::cerr << "Usage: " << std::endl;
std::cerr << argv[0] << " fixedImage movingImage " << std::endl;
return EXIT_FAILURE;
}
typedef itk::ImageFileReader< ImageType > ReaderType;
ReaderType::Pointer fixedReader = ReaderType::New();
ReaderType::Pointer movingReader = ReaderType::New();
fixedReader->SetFileName( argv[ 1 ] );
movingReader->SetFileName( argv[ 2 ] );
try
{
fixedReader->Update();
movingReader->Update();
}
catch( itk::ExceptionObject & excep )
{
std::cerr << "Exception catched !" << std::endl;
std::cerr << excep << std::endl;
}
fixedImage = fixedReader->GetOutput();
movingImage = movingReader->GetOutput();
::transform->SetIdentity();
//transform->SetIdentity();
SumOfAbsolutDistance();
//////////////////// Mean Squares //////////////////////
{
typedef itk::MeanSquaresImageToImageMetric< ImageType, ImageType >
SSDMetricType;
SSDMetricType::Pointer ssdmetric = SSDMetricType::New();
generatefile<SSDMetricType>(ssdmetric, "SSD.txt");
}
return EXIT_SUCCESS;
}
Luis Ibanez wrote:
>
>
> Hi Fab,
>
>
> Please read the code implementing the GetValue() method of the
> itkMeanSquaresImageToImageMetric.txx file. In particular, look
> at lines 77-108.
>
>
> ITK computes the metric by using only the pixels that satisfy
> *all* the following conditions:
>
>
> 1) inside the FixedImageRegion
> 2) inside the Fixed image mask (if set)
> 3) after mapping the point in to the moving image space
> it must be inside the Moving image mask (if set).
> 4) it must be inside the buffered of the moving image
>
>
> In your computation of the mean squares, you shouldn't count
> pixels that fall outside the domain of one of the image,
> simply because you don't know what should be the pixel value
> outside of the image. You could speculate it to be zero, but,
> that will be incorrect if you are dealing with a CT scan, where
> those values ouside should rather be -1024.
>
>
> Regards,
>
>
> Luis
>
>
>
> -----------------
> Fab_83 wrote:
>> Hi Louis,
>> thanks for your reply.
>>
>> 1) All my types are double.
>> 2) Thanks for the suggestion. I am using ImageIterators now. That is
>> really
>> the much easier way.
>>
>> I figured out that after the translation in the different metrics, itk
>> does
>> not consider the pixels at the edge. That means for example if I am
>> translating the moving image 5 pixels in one direction these 5pixels at
>> the
>> edge are not considered in the computing by itk. That actually means the
>> total number of calculated values is decreasing with each translation
>> step.
>> If that is right, that is probably the reason for the different values I
>> get
>> when using my own program for calculating the metrics.
>> If that is wrong would you please tell me what itk is doing with the
>> pixels
>> at the edge after the translation?
>>
>> Thanks,
>> Fab
>>
>>
>>
>> Luis Ibanez wrote:
>>
>>>Hi Fabian,
>>>
>>>Thanks for posting your code.
>>>
>>>1) What is the type of ImageType::PixelType ?
>>>
>>> In your code you are multiplying in this
>>> type before accumulating in the "double sum"
>>>
>>> sum += pixelImValue*pixelImValue;
>>>
>>> If PixelType is something that "char", then
>>> the product is overflowing the type and
>>> wrapping its value around.
>>>
>>> You should probably convert the value to
>>> RealType before performing the product.
>>>
>>> typedef itk::NumericTraits<PixelType>::RealType PixelRealType;
>>>
>>> PixelRealType realPixel = pixelImValue;
>>> sum += realPixel * realPixel;
>>>
>>>
>>>2) Just for coding style you could use ImageIterators
>>> instead of three nested loops. It will be *a lot*
>>> faster than using the GetPixel() method.
>>>
>>>
>>>Please let us know if the changes listed in (1) help,
>>>
>>>
>>>
>>> Thanks,
>>>
>>>
>>> Luis
>>>
>>>
>>>--------------
>>>Fab_83 wrote:
>>>
>>>>Hi Louis,
>>>>
>>>>enclosed I send you my code for calculating the MS. I took a look in itk
>>>>and
>>>>I think it is exactly what itk is doing. But the results are not
>>>>compareable. In this code I am reading in the differenceimages, which I
>>>>create by using LinearInterpolator and TranslationTransform (as itk).
>>>>Maybe
>>>>the problem is at the edges of the translated images. I am setting these
>>>>values zero, but actually itk is doing the same.
>>>>Nevertheless I think it is a problem with the edge because my MS values
>>>>are
>>>>increasing with every translation. And considering that the images are
>>>>just
>>>>Gausian distributed noise the MS should stay almost the same.
>>>>Maybe it is an overflow problem. My images are 512x512x150 and the
>>>>squared
>>>>values are around 1000. "Double" is only just able to handle such big
>>>>values. But itk is also using double, so the results should be the same.
>>>>
>>>>The region for calculating the image is the same. Moreover I am using
the
>>>>same Interpolator and Transform.
>>>>
>>>>//////////////////////////////////////////////
>>>>
>>>> double sum=0;
>>>> double N=0;
>>>>
>>>>
>>>>
>>>> ImageType::PixelType pixelImValue;
>>>>
>>>>
>>>> std::ofstream outputFile;
>>>> outputFile.open("SumofSquareDistance.txt");
>>>>
>>>> for (int i=1; i<NumberofInputImages; i++)
>>>>
>>>> {
>>>>
>>>> inputfilename = argv[i];
>>>>
>>>> reader->SetFileName( inputfilename );
>>>>
>>>> image = reader->GetOutput();
>>>>
>>>> reader->Update();
>>>>
>>>>
>>>> for(z=0; z<layers; z++)
>>>> {
>>>> for (x=0; x<size; x++)
>>>> {
>>>> for(y=0; y<size; y++)
>>>> {
>>>> pixelIndex[0]=x;
>>>> pixelIndex[1]=y;
>>>> pixelIndex[2]=z;
>>>>
>>>> pixelImValue = image->GetPixel(pixelIndex);
>>>>
>>>>
>>>> sum +=
>>>> pixelImValue*pixelImValue;
>>>>
>>>> N++;
>>>>
>>>>
>>>> }
>>>> }
>>>> }
>>>>
>>>> outputFile << sum/N << std::endl;
>>>>
>>>> }
>>>>
>>>> outputFile.close();
>>>>
>>>>
>>>>////////////////////////////////////////////////
>>>>
>>>>kind regards, Fab
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>
>>>>Luis Ibanez wrote:
>>>>
>>>>
>>>>>Hi Fab,
>>>>>
>>>>>Could you please post your source code to the list ?
>>>>>
>>>>>The main suspects are:
>>>>>
>>>>>A) You may not be using the same region
>>>>> of support for computing the metric
>>>>>
>>>>>
>>>>>B) You may not be using the same interpolator
>>>>> Note that intepolators perform an implicit
>>>>> smoothing of the moving image and therefore
>>>>> will alter the final results of any metric
>>>>> with respect to an instance in which you
>>>>> use a different interpolator.
>>>>>
>>>>>
>>>>>C) You may not be using the same Transform.
>>>>>
>>>>>
>>>>>
>>>>>Please let us know what Transform and Interpolator
>>>>>are you using.
>>>>>
>>>>>
>>>>> Thanks
>>>>>
>>>>>
>>>>> Luis
>>>>>
>>>>>
>>>>>-------------
>>>>>Fab_83 wrote:
>>>>>
>>>>>
>>>>>>hi everybody,
>>>>>>
>>>>>>I have used itk::MeanSquaresImageToImageMetric as a metric for my
>>>>>>registration. Then I wanted to evaluate the calculated values and
wrote
>>>>>>my
>>>>>>own routine to calculate the MS. I simply created two images with
noise,
>>>>>>made an translation and calculated the differenceimages. From these
>>>>>>differenceimages I calculated the MS with the same equation than itk
but
>>>>>>I
>>>>>>get different values than itk. The strangest thing is, that I get
values
>>>>>>which get bigger while itk gets values which stay almost the same.
>>>>>>
>>>>>>Perhabs sombody has an idea and can help..
>>>>>>
>>>>>>Thanks a lot,
>>>>>>Fab
>>>>>>
>>>>>
>>>>>_______________________________________________
>>>>>Insight-users mailing list
>>>>>Insight-users at itk.org
>>>>>http://www.itk.org/mailman/listinfo/insight-users
>>>>>
>>>>>
>>>>
>>>>
>>>_______________________________________________
>>>Insight-users mailing list
>>>Insight-users at itk.org
>>>http://www.itk.org/mailman/listinfo/insight-users
>>>
>>>
>>
>>
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
>
--
View this message in context: http://www.nabble.com/mean-square-metric-tp16418466p16582068.html
Sent from the ITK - Users mailing list archive at Nabble.com.
More information about the Insight-users
mailing list