[Insight-users] MattesMutualInformationImageToImageMetric Joint PDF update

[Luis Ibanez]

Hi TienGarry,


       Thanks for your detailed question,


It is always great when users take the time of reviewing
the source code, and it is even better when they report
back to the list about sections that look suspicious.


         Please keep doing this.   :-)


     That is what make Open Source Great !


---


In this particular case, there seems to be a section
missing from the source code that you are looking at.

In the current CVS version of this file:

  Insight/Code/Algorithms/
    itkMattesMutualInformationImageToImageMetric.txx


Lines 764- contain the following:

// Since a zero-order BSpline (box car) kernel is used for
// the fixed image marginal pdf, we need only increment the
// fixedImageParzenWindowIndex by value of 1.0.
m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] +=
  static_cast<PDFValueType>( 1 );

/**
 * The region of support of the parzen window determines which bins
 * of the joint PDF are effected by the pair of image values.
 * Since we are using a cubic spline for the moving image parzen
 * window, four bins are affected.  The fixed image parzen window is
 * a zero-order spline (box car) and thus effects only one bin.
 *
 *  The PDF is arranged so that moving image bins corresponds to the
 * zero-th (column) dimension and the fixed image bins corresponds
 * to the first (row) dimension.
 *
 */

// Pointer to affected bin to be updated
    JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
                                   ( (*fiter).FixedImageParzenWindowIndex
    * m_JointPDF->GetOffsetTable()[1] );

// Move the pointer to the first affected bin
int pdfMovingIndex =
 static_cast<int>( movingImageParzenWindowIndex ) - 1;
pdfPtr += pdfMovingIndex;

for (; pdfMovingIndex <=
       static_cast<int>( movingImageParzenWindowIndex) + 2;
       pdfMovingIndex++, pdfPtr++ )
        {


You seem to be missing the comment in lines 770-781.


This comments explains that the Parzen window used to average the
contributions of samples to the JointPDF along the moving image
direction, is a Cubic-BSpline. A Cubic-BSpline has a support of
four points. Therefore the value of the contribution to the JointPDF
is therefore distributed among four (not three) contiguous bins in
the joint histogram. Please not that the for loop variable is going
from:

       =  movingImageParzenWindowIndex - 1;

to:

      <=  movingImageParzenWindowIndex + 2;

Therefore it is visiting four locations:

    1) movingImageParzenWindowIndex - 1
    2) movingImageParzenWindowIndex
    3) movingImageParzenWindowIndex + 1
    4) movingImageParzenWindowIndex + 2


The specific values of the contributions are computed by using
the BSpline kernel:

 m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );

evaluated at the position movingImageParzenWindowArg, which is
computed as the distance between the current bin and the center
position of where the contribution to the JointPDF should have
gone (lines 798-800):

        double movingImageParzenWindowArg =
          static_cast<double>( pdfMovingIndex ) -
          static_cast<double>( movingImageParzenWindowTerm );



Please let us know if this now makes sense,
or if you have any further questions.


BTW: You may also want to update the version of ITK
     that you are looking at.



   Regards,



       Luis



----------------
TienGarry wrote:
> Dear All, 
> 
> I read the code of itkMattesMutualInformationImageToImageMetric.txx, I notice that
> 
> when the joint pdf is updated the affected bins are pdfMovingIndex-1, pdfMovingIndex , 
> 
> and pdfMovingIndex+1, while the fixedimage pdf is updated the affected bin is the FixedImageParzenWindowIndex.
> 
> The reason is the 3-order and zero-order BSpline kenel is used respectively.
> 
> Anyone can explain why the joint pdf is affected by 3 bins?
> 
> Thank you so much.
> 
>  
> 
> 00764// Since a zero-order BSpline (box car) kernel is used for
> 00765       // the fixed image marginal pdf, we need only increment the
> 00766       // fixedImageParzenWindowIndex by value of 1.0.
> 00767       m_FixedImageMarginalPDF[(*fiter).FixedImageParzenWindowIndex] += 
> 00768         static_cast<PDFValueType>( 1 );
> 00769         
> 00783       // Pointer to affected bin to be updated
> 00784       JointPDFValueType *pdfPtr = m_JointPDF->GetBufferPointer() +
> 00785                                    ( (*fiter).FixedImageParzenWindowIndex 
> 00786                                      * m_JointPDF->GetOffsetTable()[1] );
> 00787  
> 00788       // Move the pointer to the first affected bin
> 00789       int pdfMovingIndex = static_cast<int>( movingImageParzenWindowIndex ) - 1;
> 00790       pdfPtr += pdfMovingIndex;
> 00791 
> 00792       for (; pdfMovingIndex <= static_cast<int>( movingImageParzenWindowIndex )
> 00793                                 + 2;
> 00794             pdfMovingIndex++, pdfPtr++ )
> 00795         {
> 00796 
> 00797         // Update PDF for the current intensity pair
> 00798         double movingImageParzenWindowArg = 
> 00799           static_cast<double>( pdfMovingIndex ) - 
> 00800           static_cast<double>( movingImageParzenWindowTerm );
> 00801 
> 00802         *(pdfPtr) += static_cast<PDFValueType>( 
> 00803           m_CubicBSplineKernel->Evaluate( movingImageParzenWindowArg ) );
> 00804 
> 00805         }  //end parzen windowing for loop
> 00806 
> 00807       } //end if-block check sampleOk
> 00808     } // end iterating over fixed image spatial sample container for loop
> 
> 
> ------------------------------------------------------------------------
> 奥运圣火到哪里了？  点击查看！ <http://msn.live.cn/ditu>
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users