[Insight-developers] [Insight-registration] vector transformation / reorientation

[Luis Ibanez]

Hi Brian,

Please see comments below,

-----------------------
On Sat, Nov 20, 2010 at 3:35 PM, brian avants <stnava at gmail.com> wrote:

> hi everyone
>
> it's time to start thinking about how we want to transform
> vectors/tensors --- important for tensor registration, etc.  does
> anyone have strong opinions about how this should be implemented in
> itk-v4?   the important question is how to incorporate tensor/vector
> reorientation for affine & deformable transformations.
>
> some options:
>
> 1.  a filter that takes in an already resampled (by the resample image
> filter)  vector/tensor image, a transformation and a reorientation
> option (if necessary) and outputs the reoriented image.   the
> disadvantage is that the user would need to apply this filter after
> resampling and would have to know how and why this needs to be done.
>
>
Documentation & education
are always a good thing...       :-)

Note that users may not always need
the reorientation, so, they should not be
forced to pay always for that computation time.



> 2.  an implementation that "hides" the reorientation from the user.
> that is, the reorientation is applied to the tensor/vector after the
> resampling is performed but, perhaps, within the resample filter.



This implementation will require to use a Traits-like
set of helper classes to delegate the reorientation to
the pixel type itself.  That is, only the pixel type trait
class will know how is that this pixel should be oriented
under a given transformation.



>  the
> issue, here, is that we'd have to assume that a vector pixel
> represents a geometric vector / tensor.   e.g. the CovariantVector or
> SymmetricSecondRankTensor.   if someone uses the geometric vector type
> incorrectly, e.g. to represent RGB, then they will get baffling
> results.



Good point.


one nice thing about this approach is we use only one loop
> over the data for the resampling and reorientation.
>
>
Yeap.
However, how significant that benefit is,
should be measured with a profiling experiment.



> generally, spatial jacobian maps are required for these operations to
> be available.
>
>
Yeap, this will depend on the Transform.
(and in some cases, the orientation depend on the position
of the pixel in the image as well).


> Michael Stauffer pointed out that the transform base class currently
> has unimplemented TransformVector and TransformCovariantVector
> functions.


The itk::TransformBase class doesn't implement any transformation.
It is only intended to facilitate polymorphism in the process of reading
Transforms from files.

The itk::Transform class provides pure virtual functions for:


   - TransformPoint(const InputPointType  &) const = 0;
   - TransformVector(const InputVectorType &) const = 0;
   - TransformVector(const InputVnlVectorType &) const = 0;
   - TransformCovariantVector(const InputCovariantVectorType &) const = 0;


This could have been simply


   - Transform(const InputPointType  &) const = 0;
   - Transform(const InputVectorType &) const = 0;
   - Transform(const InputVnlVectorType &) const = 0;
   - Transform(const InputCovariantVectorType &) const = 0;

and we could have let the polymorphism do its job.

The reason for the redundant naming is that historically,
this class was based on VTK, where there is no type
distinction between: Point, Vector and CovariantVector
(normals), as types. (they are all represented by an array
of doubles).



> does anyone know what these were intended to implement,
> specifically?



The transformation that must be applied to a Point, is different
from what is applied to a Vector, and to a CovariantVector.


> e.g. if the transformation is an affine transformation,
> what would these functions do to the input vector?


A Point is affected by the Rotation, Scaling and Translation
components of an Affine Transform.

A Vector is not affected by the Translation part of an Affine
Transform, only the scaling and rotational part of the affine
transform. Because a vector is a difference between two
points.

A CovariantVector (e.g. the cross product of two vectors) is not
affected by the Translation part of an affine Transform either,
it is affected by the rotational part and by the reciprocal of the
scaling part. (in practice, the inverse of the transpose matrix).



> for instance, the
> functions may only resample the vector into the new space and not
> reorient them.   alternatively, the vectors may be both transformed
> and reoriented under rotation.  a third option is that they would be
> transformed, reoriented and scaled under the full affine
> transformation.  all of these use cases may come up.  and would the
> behavior for TransformVector and TransformCovariantVector be the same?
>
>
Nope,

the behavior of a covariant vector is different from a vector under
affine transformations.

Think of the normal of a Sphere (they are covariant vectors)
If you scale the sphere along Z by a factor of 0.5 (to make it
like an ellipsoid whose Z axis is half of its X,Y axes, then the
Z components of the normals will be divided by 0.5 (multiplied
by 2), while the Z components of the Vectors that connect
points in the surface of the sphere with the center, will be
multiplied by 0.5.

This is probably better explained here:
http://www.cgafaq.info/wiki/Transforming_Normals


Note that covariant vectors are very common:

  The gradient of a scalar image
  is an image whose pixels are covariant vectors.




> my feeling is that vectors should be resampled by exploiting the same
> functionality as one would for a scalar.  then reorientation should be
> applied afterward via a separate operation either exposed to the user
> (option 1 above) or applied implicitly (option 2).   if we choose
> option 2, then the datatype would have to encode the type of
> reorientation that would be expected.  i do not think current itk data
> types are sufficiently specific to support this.
>
>

Performing first the resampling of the {vector, point, covariantvector}
values and then following it by the potential correction of orientations
seems to be the right thing to do.

I agree that in case (2), it should be up to the pixel type (or its Traits)
to define how is that it should be affected under a given Transformation.

Unfortunately, the pixel type doesn't know (and it shouldn't know)
about the family of Transformations (affine, rigid, bspline...), so
it is therefore up to the Transforms to define how is that they affect
each one of the potential pixel types.

We can probably specify this behavior for


* Vectors
* Points
* CovariantVectors
* Tensors (of different ranks and combination
         of covariant and contravariant indices)

(after alll,
Vectors are tensors of rank 1 with a single contravariant index,
and CovariantVectors are tensors of rank 1 with a single
covariant index).

This reorientation probably should be done in a specific
filter whose goal is to reorient the pixel content of
the image, instead of reorienting the grid itself.


Note that something similar to this has already been
done implicitly in the Metrics and the interpolators:

See:

ITK/Code/Algorithms:

itkImageToImageMetric.txx:434:
m_DerivativeCalculator->UseImageDirectionOn();
itkImageToImageMetric.txx:445:
m_BSplineInterpolator->UseImageDirectionOn();
itkImageToImageMetric.txx:808:  gradientFilter->SetUseImageDirection(true);
itkMutualInformationImageToImageMetric.txx:51:
m_DerivativeCalculator->UseImageDirectionOn();

Which is the equivalent of reorienting vectors
as part of the interpolation of the images.


if there are no opinions about this, then we will probably lean
> towards option 1 which will force people to think about what type of
> reorientation they should use for their data.   on the other hand, for
> some cases where data type is specific enough (e.g. DTI) we would
> prefer option 2.
>
>

I will vote for option (2): Creating a separate filter for
reorienting the internal (vectors, covariant vectors) of
an image.

Based on the following arguments:

A) Users should only pay for computation
     time of things they need.

B) The generic implementation of option 1 is much more
     difficult than option 2.



Also, coming back to Michaels original question:

How should Vectors and CovariantVectors be transformed
under deformation fields (whether they are dense or sparse
and interpolated with BSplines or Kernel transforms) ?

Here are two potential options:

(a) Consider that a Vector is the relative position between
     two points. Therefore, we use the Transform to map
     the tail and head of the vector, and then rebuild the vector
     in the output space by subtracting the images of the tail
     and head points. A similar geometrical construct could
     be done for the covariant vector (but involving N-1 points)
     where N is the dimension of space.

or

(b) Estimating the Rotational and Stress component of the
      deformation field at that pixel location and use it to correct
      the Vector (or covariant vector) accordingly.  This is
      essentially based on the Jacobian matrix of the deformation
      field, estimated at the location of the pixel.


Method (b) is probably more computationally expensive,
but seems to be the right mathematical approach.


How about scheduling a tcon early next week to
brainstorm about the pros and cons of these options ?



      Luis


----------------------------

> As always,
>
> Brian
>
>
> On Fri, Nov 19, 2010 at 3:59 PM, Michael Stauffer (Circular Logic)
> <mstauffer at circular-logic.com> wrote:
> > Hi Brian,
> >
> > The latest itkTransform class has new (pure virtual) transform methods,
> > see below. What does it mean to transform a vector or covariantvector
> > with a deformation field? for now, I'll just implement them and throw an
> > exception so I can get the rest of the code up for starting review.
> >
> >  virtual OutputVectorType    TransformVector(const InputVectorType &)
> > const = 0;
> >  virtual OutputVnlVectorType TransformVector(const InputVnlVectorType
> > &) const = 0;
> >  virtual OutputCovariantVectorType TransformCovariantVector(const
> > InputCovariantVectorType &) const = 0;
> >  virtual void SetParameters(const ParametersType &) = 0;
> >  virtual void SetParametersByValue(const ParametersType & p)
> >  { this->SetParameters (p); }
> >  virtual void SetFixedParameters(const ParametersType &) = 0;
> >
> > -M
> >
> >
>
>
>
> --
> ß®∫∆π
> _______________________________________________
> Insight-registration mailing list
> Insight-registration at public.kitware.com
> http://public.kitware.com/cgi-bin/mailman/listinfo/insight-registration
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/mailman/private/insight-developers/attachments/20101120/af66b580/attachment.htm>