[Insight-users] mean square metric
Luis Ibanez
luis.ibanez at kitware.com
Thu Apr 3 11:08:36 EDT 2008
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
>>
>>
>
>
More information about the Insight-users
mailing list