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

[Christoph Niedermayr]

Am Montag, den 04.09.2006, 14:54 -0400 schrieb Luis Ibanez:

>                   This is not a bug.

i still believe it is :)

> 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.

this is correct, but only applies to SampleFixedDomain().
As i said, i use UseAllPixelsOn(), which causes SampleFullFixedDomain()
to be called instead.
if you take a look at [1], you will see that SampleFullFixedDomain()
uses a ImageRegionConstIteratorWithIndex.

thanks for the detailed explanation though.

best regards,
chris


[1] itkMattesMutualInformationImageToImageMetric.txx, line 487
template < class TFixedImage, class TMovingImage >
void
MattesMutualInformationImageToImageMetric<TFixedImage,TMovingImage>
::SampleFullFixedImageDomain( FixedImageSpatialSampleContainer&
samples )
{

  // Set up a region interator within the user specified fixed image
region.
  typedef ImageRegionConstIteratorWithIndex<FixedImageType>
RegionIterator;
  RegionIterator regionIter( this->m_FixedImage,
this->GetFixedImageRegion() );


> 
> 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
> > 
> > 
> 
>