[Insight-developers] How dose the MeanSquaresImagetoImageMetric compute the derivative?

[Luis Ibanez]

Hi Xingxing,


1) Metric = cost function that takes as arguments the array
             of Transform parameters, and returns a scalar value.

    Let's call them

        P = array of paramters from the Transform
        T = transform
        M = metric

    so we can say that one evaluation of the metric is described
    by M(P)


2) The Mean squares metric is computed by mapping points from
    the fixed image coordinate system into the moving image
    coordinate system.


3) When we look for the derivative of the metric we are talking
    about


    D(M(P) = [ d(M(P))/dP0, d(M(P))/dP1, ...d(M(P))/dPn]

    That is, the derivative of the metric is the array whose
    components are the variation of M(P) with respect to each
    one of the transform components.


4) The way a transform component Pi affects the computation
    of the metric is by making the points from the fixed image
    land on slightly different locations in the moving image.

    The Jacobian of the Transform describes how much the
    coordinates of the mapped point will vary as a function of
    the changes in the Transform parameters:

       pm = point in the moving image
       pf = point in the fixed image

       pm = T( pf )

    The Jacobian Matrix of T with respect to P is then given
    by the (i,j) components:

       Jij  =  d(pm(i)) / dPj

     That is:

       Jij = when mapping the point pf, from the fixed image
             to the moving image, how much the coordinate "i"
             of the point in the resulting pm point will change
             if we perturbe the transform parameter "j".
             moving image


5)  The effect of perturbing the landing point of pm in the
     moving image is that now the metric is going to use a
     slighlty different intensity value for the contribution of
     this point to the computation of the mean square differences.

     Let's call

        FI = Fixed Image
        MI = Moving Image

     The amount of this intensity change is given by the gradient
     of the moving image intensities = GMI, which is a covariant
     vector of N components (N = image dimension).


6)  The computation of the Mean Square metric is given by


              M(P) = (1/NP) Sum ( FI(pf) - MI( pm ) )^2

      where NP = total number of pixels

      which can be rewritten as

              M(P) = (1/NP) Sum ( FI(pf) - MI( T(P,pf) ) )^2


7)  You can now put this together using the chain rule:


      d(M(P))/dPi =

        (1/NP) Sum_j [ ( FI(pf) - MI( T(P,pf) ) * Sum_i ( Ji * GMI_i )


      Sorry for the notation,
      but this is as far as I can get without using LaTeX.


8) The code that implements the expression in numeral (7)
    is in the file:

        src/Insight/Code/Algorithms/
              itkMeanSquaresImageToImageMetric.txx


     in lines 196-228, whose essential pieces are:


       const TransformJacobianType & jacobian =
         this->m_Transform->GetJacobian( inputPoint );

      const GradientPixelType gradient =
         this->GetGradientImage()->GetPixel( mappedIndex );

       for(unsigned int par=0; par<ParametersDimension; par++)
         {
         RealType sum = NumericTraits< RealType >::Zero;
         for(unsigned int dim=0; dim<ImageDimension; dim++)
           {
           sum += 2.0 * diff * jacobian( dim, par ) * gradient[dim];
           }
         derivative[par] += sum;
         }



9) For more details, please read the ITK Software Guide

      http://www.itk.org/ItkSoftwareGuide.pdf

    Section 8.8.2, "Transform General Properties",
    in pdf page 427.





        Regards,



           Luis




------------------------
Xingxing Pan wrote:
> I've spent the whole day reading the the source code of the metric 
> class and its ancestors, trying to understand the machanism of 
> computering the derivative. The Jacobian matrix and the derivative image 
> of the moving image are impenetrable.
> Could someone give me a explanation of the machnism or a link to some 
> reference materials?
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> Powered by www.kitware.com
> 
> Visit other Kitware open-source projects at http://www.kitware.com/opensource/opensource.html
> 
> Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ
> 
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-developers