[Insight-developers] OrientedImage and gradient calculations

Miller, James V (GE, Research) millerjv at crd.ge.com
Fri Mar 24 09:17:42 EST 2006


Stephen,

I think of it as follows:  the PixelContainer class is really the 
data container. The Image places that container in context.

If we factored out these coordinate manipulations, 

1) would Image still store spacing, origin, and directions?
2) could you write an algorithm (say a registration) abstractly
    such that if an image did not have/use orientations (like 
    Image does not) the (new) external transformation functions 
    would take into account just spacing, and if an image did 
    have/use orientations(like OrientedImage does) the (new)
    external transformation functions would take
    into account spacing, origin, and orientation?


This actually brings up another issue.  Currently, the methods

    A) TransformPhysicalPointToContinuousIndex()
    B) TransformContinuousIndexToPhysicalPoint()

are templated over the coordinate representation.  As such they
are not virtual.  This means that if you have an OrientedImage
cast back to an Image, and you run a registration on that Image, 
the registration will NOT take into account orientation. This is
probably not good. It is understandable how this design came about
because at the time the decision was made to have the methods
TransformPhysicalPointToContinuousIndex, we did not have the 
subclass OrientedImage.

Unfortunately, I don't think we can have a virtual templated member
function. We also cannot "untemplate" this method and provide
overloaded versions (for double and float coordinate representations)
since that may break backward compatibility.  Someone may have code
like

  mImage->TransformPhysicalPointToContinuousIndex<float>(index, point);

which would not compile if we "untemplated" the code. I also don't
know whether we could now provide untempleted overloaded versions
like

  virtual void TransformPhysicalPointToContinuousIndex( const ContinuousIndex<float, VImageDimension>&,
                                                Point<float, VImageDimension>&) const;
  virtual void TransformPhysicalPointToContinuousIndex( const ContinuousIndex<double, VImageDimension>&,
                                                Point<double, VImageDimension>&) const;

since these overloaded versions probably conflict with the templated member
function versions.  We could test this.

We could perhaps add new methods with slightly different names that do what we want 
in term of polymorphism and change the implementations of the templated member function 
versions to call the new virtual function versions.

Jim



-----Original Message-----
From: Stephen R. Aylward [mailto:Stephen.Aylward at Kitware.com]
Sent: Thursday, March 23, 2006 7:48 PM
To: Luis Ibanez
Cc: Miller, James V (GE, Research); Blezek, Daniel J (GE, Research);
Insight-developers (E-mail)
Subject: Re: [Insight-developers] OrientedImage and gradient
calculations



I see the convenience of the transform in the image...and we've already 
headed down that road with point transforms....but doesn't "proper" c++ 
programming style suggest keeping itk::Image as a data container and 
putting functions in appropriate classes?

Stephen

Luis Ibanez wrote:
> 
> 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
>>
>>
> 
> _______________________________________________
> Insight-developers mailing list
> Insight-developers at itk.org
> http://www.itk.org/mailman/listinfo/insight-developers
> 

-- 
=============================================================
Stephen R. Aylward, Ph.D.
Chief Medical Scientist
Kitware, Inc.
http://www.kitware.com


More information about the Insight-developers mailing list