[Insight-developers] RE: Demons voodoo

[Luis Ibanez]

Hi Jim,

Here are some comments regarding the implementation
of the Demons algorithm.


 >
> According to the demon's paper (page 15 and equation 4 on page 17), the goal
> of the demons algorithm is to "make an image g appear more like an image f".
> Equation 4 is essentially
>  
>     v = (g -f) grad(f) / (gradmag(f)^2 + (g - f))
> 


1) The equation 4 in the INRIA report paper is actually:

       v = (g -f) grad(f) / (gradmag(f)^2 + (g - f)^2 )

   Note (g-f) are squared in the denominator.

   This equation (in the paper) was unbalanced in units in the
   denominator. This was fixed in the ITK implementation by
   introducing a factor K that is equal to the squared diagonal
   of a pixel in physical units.


About the interpretation of the deformation field, you certainly
can make both cases. (A) say that V is the vector that will move
a pixel in "f"(fixed image) onto the corresponding pixel in "g"
(moving image). (B) say that -V is the vector that will move the
corresponding pixel of "g" into this pixel of "f".

As Lydia pointed out, the second option is the one that allows
to use the deformation field for resampling "g" and produce an
image in the coordinate system of "f" that doesn't have holes.

That's why the implementation is using -V, as you showed in your
second point

 >
>  dvi =  (f - g) grad(f) / (gradmag(f)^2 + (f - g))
>  

   (where (f-g) in the denominator should be (f-g)^2 )

So the only difference between the first equation and
the second is the sign of the vector.

>
> Mapping the code to the equations above, f is the fixed image and g is the
> moving image. Note that this is calculating a vector at a point pi that is
> associated with image f.  This seems to be calculating a vector field that
> is keeping track of "where a pixel in f goes" as f is warped to g.
>  

This is the case for all the transforms in  the registration framework.
The transform is always telling for every pixel of the fixed image,
where to find its corresponding pixel in the moving image. This is how
at resampling time we can walk over the space of the fixed image and
find corresponding pixels to resample from the moving image, therefore
generating a mapped moving image without holes nor duplicated pixels.

The fact that the gradient of "f" is used, indicates that the fixed
image space is the one that must be used as reference frame.
The current implementation allows to take the displacement field
directly from the Demons filter and use it in the WarpImageFilter.

--

As Thirion mention in the paper, "Demons" is a generic name for a
family of algorithms that apply local forces in order to compute
deformations, there are many possible implementatios of Demons
algorithms.  In Thirion's paper there are at least three different
possible implementations.  Corinne Matmann recently contributed
a variant that is using equation (5), this was named itkSymmetric
ForcesDemonsRegistrationFuntion/Filter.  Using equation (9) will
lead to another variant, for which we simply will have to figure
out an appropriate name.

Your suggested implementation can be easily constructured from
Corinne's variation by switching the interpretation of "f" and "g"
and removing the use of one of the gradients. This should probably
be introduced then as a third type of Demons-like algorithm.

We could call it something like

     itkForwardDemonsRegistrationFunction
     itkForwardDemonsRegistrationFilter

  or

     itkBackwardDemonsRegistrationFunction
      (depending on what your preference is for interpreting
       the direction of the displacement field).


This interpretation will not be suitable for resampling,
however users in the list have asked for such functionlity
in order to do things that do not require resampling.




     Luis