[Suspected Spam]Re: [Insight-developers] OrientedImage and gradient calculations

[Miller, James V (GE, Research)]

Luis, 

You bring up a point that I did not include in the original message, namely
by having the following methods, 

   A) TransformIndexGradientToPhysicalGradient()
   B) TransformLocalGradientToPhysicalGradient()

we now have named the following coordinate systems:

   1) Index coordinate frame.  Integral I, J, K frame.
   2) ContinuousIndex coordinate frame.  Same as the Index coordinate frame
       but with floating point precision.
   3) Physical coordinate frame.  A world coordinate frame that
       accounts for pixel spacing and orientation.
   4) Local coordinate frame.  A pseudo-world coordinate frame that
       accounts for pixel spacing but not orientation.

The term "Local" in this context would be new to ITK. We should make sure
that this is the name we want to associate with this coordinate
frame.

We should probably revisit our gradient calculation filters and either:

   i) Document that they only calculate the gradient wrt the local coordinate
        frame
  ii) Modify the filters so that they can produce gradients in any coordinate
        frame (Index, Local, or Physical).


There are other places in the toolkit where gradients are calculated in place
(without using a filter). For instance, in some of the level sets and optical
flow based registration.  We'll need to determine whether these gradients
are being calculated in the right space.  Right now, "most" places are 
calculating the gradient in the Local coordinate frame.  Before we added
orientation, these filters "thought" they were calculating the gradient
in the Physical coordinate frame.  For the level sets and anisotropic 
diffusion, the gradients should probably in the Local coordinate frame.
For the optical flow based registration (demons, levelset motion, etc)
and perhaps in the FEM registration, these gradients should probably in the 
the Physical coordinate frame.

The more I keep writing "Local coordinate frame" the more I like the term.
Though I am still a bit concerned as to whether someone new to ITK would
recognize this frame as taking into account spacing or not.

Jim



-----Original Message-----
From: Luis Ibanez [mailto:luis.ibanez at kitware.com]
Sent: Thursday, March 23, 2006 6:14 PM
To: Miller, James V (GE, Research)
Cc: Peter Cech; Insight-developers (E-mail); Blezek, Daniel J (GE,
Research)
Subject: [Suspected Spam]Re: [Insight-developers] OrientedImage and
gradient calculations



Hi Jim,


Adding a TransformIndexGradientToPhysicalGradient() to the itk::Image
sounds like a good idea. It sounds better than correcting inside
the filters because Gradient operations can be useful in the local
coordinate system for the purpose of performing operations such as
segmentation and feature detection.


The Image gradient that is computed by the Image Metrics is
produced by the GradientRecursiveGaussianImageFilter. This
filter already takes into account the image spacing. So the
proposed Image method

       TransformIndexGradientToPhysicalGradient()

doesn't need to apply any spacing correction... However, other gradient
filters may not be taking the spacing into account (something to double
check...).


As Jim pointed out, this problem will apply to all other filters
that compute gradients in the image, or higher order derivatives.


In the most generic case we should use a Metric Tensor for mapping
the derivatives from one space into the other. The reason why the
Tensor notation may be more appropriate is that for higher order
derivatives the Tensor has to be applied multiple times.

Derivatives are covariant tensors, therefore a Gradient, is a
Tensor of covariant rank one, while a Hessian is a Tensor of
covariant rank two.

The MetricTensor that maps the image space into the world
coordinates space can be obtained as the inverse of the matrix
used in the itkImageTransformHelper.

In the particular case of the GradientRecursiveGaussianImageFilter,
we only need the correction of the orientation, so the matrix
that will be useful is the m_Direction matrix instead of the
m_IndexToPhysicalPoint matrix which also contains the scaling
corrections for pixel spacing.

One way of avoiding the ambiguity and potential bugs of use,
is to provide two new methods along the lines of what Jim
suggested:

   A) TransformIndexGradientToPhysicalGradient()
   B) TransformLocalGradientToPhysicalGradient()

The first one will use a matrix that involves both rotation
and pixel spacing, while the second use a matrix that involves
only rotations. We probably want to add also the inverse methods:

   C) TransformPhysicalGradientToIndexGradient()
   D) TransformPhysicalGradientToLocalGradient()

and eventually:

   E) TransformLocalGradientToIndexGradient()
   F) TransformIndexGradientToLocalGradient()

although for these four methods we could wait until we see
a real need for them...



I would suggest to add only the two methods (A) and (B)
as Jim proposed.



    Luis



Thanks Dan for



        "using your deductive powers to
         figure out an extremely vexing ITK bug".




==================
Peter Cech wrote:
> On Fri, Mar 17, 2006 at 09:49:17 -0500, Miller, James V (GE, Research) wrote:
> 
>>Dan Blezek stumbled across the following problem:
>>
>>In the registration framework, the metrics use the image gradient in the calculation of the derivative of the metric wrt to the parameters of the transformation.  The image gradients are calculated taking into account the image spacing to provide a gradient in physical space.
>>
>>When an OrientedImage is used, these gradient calculations are not truly in physical space since gradients are not reoriented into the physical coordinate frame.  This affects the registration because the mapping of positions does take into account orientation but the image gradients, and hence the derivative of the metric wrt the parameters, does not.  In Dan's test case, when the derivative of the metric should force the transformation to move the image up, it ends up moving the image down.
>>
>>One solution would be to add methods to image like:
>>
>>  template<class TCoordRep>
>>  void TransformIndexGradientToPhysicalGradient(
>>                      const IndexType & index,
>>                      CovariantVector<TCoordRep, VImageDimension>& gradient) const
>>
>>which would convert a gradient computed in index space (without taking into account spacing or orientation) to a gradient computed in physical space.
>>
>>In Image, this method would merely apply the spacing scale factors.  In OrientedImage, this method would apply the spacing scale factors and the orientation matrix.  Anywhere gradients are calculated, we would have to make an extra function call to convert the gradients into physical space.
>>
>>We may also need methods for transforming higher order derivatives (and cross-derivatives).
> 
> 
> Definitely. I just realized, that I recently witnessed the behavior when
> displaying eigenvectors of hessian. I thought it was off because the
> scale did not fit well for the structure...
> 
> 
>>Alternatively, wherever gradients are calculated, we could check whether we are operating on an Image or an OrientedImage and reorient the gradients as necessary there.
> 
> 
> If we stay consistent with current Index to Point behavior (1) of Image,
> a simple pass-through implementation will do.
> 
> (1) It assumes origin is in the same orientation as the Image.
> 
> Regards,
> Peter
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
> 
>