[Insight-users] mean square metric
Luis Ibanez
luis.ibanez at kitware.com
Sun Apr 6 17:06:43 EDT 2008
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
>>
>>
>
>
More information about the Insight-users
mailing list