[Insight-users] BUG in MattesMI (FixedImageMask, UseAllPixelsOn)

[Luis Ibanez]

Hi Chris

                  This is not a bug.


The RandomImageIterator that is used in the implementation
of MattesMutualInformation does not move in an orderly manner
across the Image Buffer, but it never gets out of the region
that it was defined to walk.

Every invocation of the itr++ operator is simply equivalent
to rolling a dice for each one of the index components in the
iterator position.

Calling the operator++ many times will never result in the
index to go out of the image domain. The reason is that the
dice rolling is done by taking the image buffer as a continuous
linear array, and then every call to ++ generates a number between
zero and the number of pixels in the image. The result is used
for computing the index of the pixel that correspond to that
offset position from the origin of the image buffer.

You may find interesting to look at the code of the

    Insight/Code/Common/
      itkImageRandomConstIteratorWithIndex.h
      itkImageRandomConstIteratorWithIndex.txx

where the operator++() is implemented.

You will find that it calls the RandomJump method, that in turn
uses the MersenneTwisterRandomVariateGenerator for generating a
random position *inside* the range of the image region.


Please take a look at this code and let us know if you have
further questions,



    Thanks


       Luis


--------------------------
Christoph Niedermayr wrote:
> Hi all!
> 
> I wanted to use the MattesMutualInformationImageToImageMetric with a
> fixed (and moving) image masks.
> The metric should use all pixels in the mask so i set UseAllPixelsOn().
> After some hours of debugging :( i found the following bug:
> 
> in [1], the metric prepares for taking as much samples as the size of
> the fixed image region
> in [2], the whole fixed image is sampled, but "bad" samples outside the
> mask are skipped (thats ok). the sample iterator only steps forward when
> a "good" sample has been found. the loop only stops when the sample
> container is full (thats bad).
> 
> this means that for each sample outside the mask, the iterator will step
> one index *outside* the image buffer!
> 
> i hope i'm not mistaken!
> chris
> 
> 
> [1]
> itkMattesMutualInformationImageToImageMetric.txx, line 237:
>  if( m_UseAllPixels )
>     {
>     m_NumberOfSpatialSamples = 
>               this->GetFixedImageRegion().GetNumberOfPixels();
>     }
>   
>   /**
>    * Allocate memory for the fixed image sample container.
>    */
>   m_FixedImageSamples.resize( m_NumberOfSpatialSamples );
> 
> 
> 
> [2]
> itkMattesMutualInformationImageToImageMetric.txx, line 501:
>   if( this->m_FixedImageMask )
>     {
> 
>     typename Superclass::InputPointType inputPoint;
> 
>     iter=samples.begin();
> 
>     while( iter <= end )
>       {
>       // Get sampled index
>       FixedImageIndexType index = regionIter.GetIndex();
>       // Check if the Index is inside the mask, translate index to point
>       this->m_FixedImage->TransformIndexToPhysicalPoint( index,
> inputPoint );
> 
>       // If not inside the mask, ignore the point
>       if( !this->m_FixedImageMask->IsInside( inputPoint ) )
>         {
>         ++regionIter; // jump to next pixel
>         continue;
>         }
> 
>       // Get sampled fixed image value
>       (*iter).FixedImageValue = regionIter.Value();
>       // Translate index to point
>       (*iter).FixedImagePointValue = inputPoint;
> 
>       // Jump to random position
>       ++regionIter;
>       ++iter;
>       }
>     }
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
> 
>