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

[Luis Ibanez]

Hi Christoph,


       You are right,  This is indeed a bug.

I missed to notice that the while loop was controlled by
the end() of the samples container instead of the IsAtEnd()
of the region iterator.


Karthik already has a fix for it and it will be committed
before ITK 3.0.



     Regards,



         Luis


----------------------------
Christoph Niedermayr wrote:
> 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
>>>
>>>
>>
>>
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
> 
>