[Insight-users] Gradient computation bug in SingleValuedVnlCostFunctionAdaptor when scales are used (patch included)

[Luis Ibanez]

Hi Tom,

Thanks for pointing this out.

You are right, the current code doesn't perform the correct
mathematical operation that corresponds to scaling the argument
of a cost function.


However, that was not the intent of introducing the scaling
in these classes.


The motivation for this scaling operation was to extend the
functionality of VNL optimizers to be usable with domains outside
of the cannonical [-1:+1] domain.


Since these functions are used for the purpose of optimization,
we intended simply to shrink (or expand) the domain of the original
problem in order to map to the [-1:+1] that the algorithmic
implementation of vnl optimizers expect.

To make sense of it, you may find interesting to look at both the
ITK wrappers for VNL optimizer, and the ITK wrapper for the cost
function.


Take for example,
the itkAmoebaOptimizer.cxx:

   it its call to GetValue() in line 84, it uses the scales
   to multiply the parameters *before* they are passes to
   the cost function.

      84:      parameters[i] *= scales[i];

   In its call to StartOptimization(), before invoking the
   VNL optimizer, we scale the parameters of the transform

     218:       parameters[i] *= scales[i];

   then we invoke the VNL optimizer

     227:      m_VnlOptimizer->minimize( parameters );

    and then we restore the scale of the transform parameters

    242:       parameters[i] /= scales[i];


The goal, is to create a situation where the VNL optimizer
"thinks" that is optimizing a function whose domain is close
to [-1:+1], as VNL algorithms expect or assume.


However, the VNL optimizer is calling the cost function
in order to get values from it. That cost function will
happen to be an ITK ImageMetric (most of the time), and
therefore the VNL cost function adaptor must "restore"
the scale of the parametric space to what the metric
expects as parameters. That's why the cost function adaptor

     itkSingleValuedVnlCostFunctionAdaptor.cxx

in its call to "f()" rescales the parameters *back* by using
in line 61:

    61:   parameters[i] = inparameters[i]/m_Scales[i];

and then calling the actual ITK cost function in line

    69:   m_CostFunction->GetValue( parameters ));


In this way, ITK metrics will still work in a space where
the domain of translations may be [-500mm, +500mm] and
[-0.1degrees,+0.1degrees], while the VNL optimizer works
in a space with a domain close to [-1:+1].


The point is well taken in that we could have applied the
scaling to the resulting values of the gradients. Something
that could conveniently be done in the method:

        ConvertExternalToInternalGradient()

in line 176 of itkSingleValuedVnlCostFunctionAdaptor.cxx


This additional (and mathematically correct scaling) will
indicate the following:


    A) if we have an image of 100x100 pixels with spacing
       2.0mm x 2.0mm and we are optimizing the Mean square
       metric against another image of similar dimensions,
       we see this problem happening in a domain of
       [200x200]mm.  If the Metric changes by 35 intensity
       values in average when you translate the images by
       5 millimeters, that corresponds to a derivative of

                        7 iv / mm


    B) The scaling conversion will be equivalent to replace
       the units of mm with units of "image width".  That is
       somehow normalizing the space to the size of the image.

       In this context, the conversion

                F(x) = f(x/s)

                grad(F(x)) = grad(f(x/s))/s

       will give us the value of the gradient on the units
       of "image width", as:

                     (7/200)  iv / iw



      At this point, we could make this correction,
      but this may require to recalibrate the parameters
      of most of image registration examples and tests
      that use VNL optimizers...


Have you run an experimental build with this change
to evaluate the impact of the modification ?


----


The operation is not done in the same way either across optimizers.

For example,
if you look at the GradientDescent and RegularStepGradientDescent
optimizers, the scales are applied only to the *resulting* gradient
vector returned from the cost function, and then, the resulting vector
is used during the computation of the step to be taken. In these two
optimizers, the scale factor is not used for feeding the values
inside the cost function.



     Regards,


         Luis


-----------------------
Tom Vercauteren wrote:
> Hi all,
> 
> I think I found a bug in
> SingleValuedVnlCostFunctionAdaptor::ConvertExternalToInternalGradient.
> Here is a brief description:
> 
> SingleValuedVnlCostFunctionAdaptor can be used to wrap an itk cost
> function into a vnl cost function.
> 
> Additionally the user can specify a set of scales such that the
> function to be optimized is F(x) = f(x./s) instead of f(x).
> 
> The gradient of F, grad(F)(x) is thus given by grad(F)(x) = grad(f)(x./s)./s
> 
> Instead of that, the cost function adpator returns grad(f)(x./s)
> 
> 
> I have filed a bug report in the bug tracker:
> http://www.itk.org/Bug/view.php?id=6887
> 
> Let me know if I am missing something.
> 
> Best regards,
> Tom Vercauteren
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>