[Insight-developers] OrientedImage and gradient calculations

[Luis Ibanez]

Jim, Stephen,

We should avoid using virtual functions for this
conversion because virtual functions do not get inlined.

The conversions between Indexes / Local Points / Physical Points
are critical for the performance of the registration framework.

These functions are called *millions of times* per iteration of
the registration optimizer.

Polymorphism at this point may be too costly and may hinder
what has been gained by the use of template meta programming.

I don't see too much of a problem in adding this operations
to a helper class that is friend of the Image. The helper class
would play about the same role as the ImageCalculators or the
TransformInitializers.


    Luis



--------------------------------------
Miller, James V (GE, Research) wrote:
> 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
>>
> 
>