Discussion:
[Insight-developers] vector transformation / reorientation
brian avants
2010-11-20 20:35:13 UTC
Permalink
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.

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. 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. one nice thing about this approach is we use only one loop
over the data for the resampling and reorientation.

generally, spatial jacobian maps are required for these operations to
be available.

Michael Stauffer pointed out that the transform base class currently
has unimplemented TransformVector and TransformCovariantVector
functions. does anyone know what these were intended to implement,
specifically? e.g. if the transformation is an affine transformation,
what would these functions do to the input vector? 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?

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.

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.

As always,

Brian


On Fri, Nov 19, 2010 at 3:59 PM, Michael Stauffer (Circular Logic)
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
--
ß®∫∆π
Matthew McCormick (thewtex)
2010-11-20 21:13:45 UTC
Permalink
Post by brian avants
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.
Hi Brian,

I prefer option 2 because it is too easy to forget to apply the
transformation afterwards. There are many specific pixel types in
ITK, and we can take advantage of this. If the wrong pixel type was
chosen, then the use of the correct pixel type should be advocated
when the choice is made.

Regards,
Matt
Luis Ibanez
2010-11-20 23:18:43 UTC
Permalink
Hi Brian,

Please see comments below,

-----------------------
Post by brian avants
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.
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.
Post by brian avants
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.
Post by brian avants
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
Post by brian avants
over the data for the resampling and reorientation.
Yeap.
However, how significant that benefit is,
should be measured with a profiling experiment.
Post by brian avants
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).
Post by brian avants
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).
Post by brian avants
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.
Post by brian avants
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).
Post by brian avants
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.
Post by brian avants
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
Post by brian avants
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


----------------------------
Post by brian avants
As always,
Brian
On Fri, Nov 19, 2010 at 3:59 PM, Michael Stauffer (Circular Logic)
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
http://public.kitware.com/cgi-bin/mailman/listinfo/insight-registration
Matthew McCormick (thewtex)
2010-11-21 00:00:34 UTC
Permalink
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
after I get another paper out and the following patch gets merged:

http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z

Matt
Luis Ibanez
2010-11-21 00:07:05 UTC
Permalink
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Matt
Excellent !

Yes,
that what I meant by the stress component.

Do you see a way of deriving as well
the rotational component of the Jacobian ?


Luis
Matthew McCormick (thewtex)
2010-11-21 00:14:48 UTC
Permalink
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help.  I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains.  I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian  ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details

And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?

Matt
Luis Ibanez
2010-11-21 00:36:03 UTC
Permalink
On Sat, Nov 20, 2010 at 7:14 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details
And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?
Matt
Yes,
That sound very reasonable to me.


Luis
brian avants
2010-11-21 01:43:59 UTC
Permalink
Right. The eigendecomposition approach is part of the "preservation of
principal directions" algorithm used when reorienting diffusion tensors.
The idea is that, as long as your deformation is a diffeomorphism, you can
estimate it, locally, as an affine transformation. Then you use the rigid
part of the affine transformation to reorient the tensor's eigenvectors.
There are a few other steps, though.

My question, in part, is whether ITK's pixel types are specific enough to
uniquely specify the required reorientation. Some examples:

A diffusion tensor may be reoriented by either PPD (above) or finite strain
estimates. So this is a non-unique mapping from pixel type to
reorientation strategy.

When reorienting a vector, one may want to change its relative orientation
but not change its length, even under an affine transformation. Also a
non-unique mapping from pixel type to reorientation strategy.

One approach (what we've done in the past) is to implement filters for each
of these strategies for DTI. While having filters would be nice, I can see
this becoming confusing, in particular if some metrics/interpolators are
performing reorientation as well.
Matt, I look forward to seeing the strain filter. Thans much!

Re: Tcon , when would you like to discuss? Previously, we used 10 am
Tuesday for registration discussions.

Brian
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:14 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?
Matt
Yes,
That sound very reasonable to me.
Luis
Tom Vercauteren
2010-11-22 08:41:30 UTC
Permalink
Hi all,

As Brian just mentioned, there may be several ways of performing
re-orientation for a given type. PPD or FS being the most widely
encountered choices for diffusion tensors reorientation. Besides, what
can be of interest, especially for image registration, is the the
derivative of this reorientation strategy. It should indeed come into
the computation of the derivative of the objective function.

[Shameless add] For the specific case of FS in dimension 3, the
analytical derivative can be found in Thomas Yeo's paper (which I
co-authored):
http://people.csail.mit.edu/ythomas/publications/2009DTI-TMI.pdf

Hence, I believe that we need a generic approach where one can plug
its own reorientation strategy. In that view, I would prefer a 2bis
option where resampling and reorientation are done in one pass but
where the user can choose its reorientation strategy much like he can
now choose its interpolator.

The tricky part is how to handle the derivative of the reorientation
strategy in a generic manner in the registration framework and not
kill performance.

If we add a good linear algebra package (say eigen 3*), we could have
a no-op reorientation strategy that simply multiplies by identity
matrices but does not cost anything at run-time thanks to expression
templates.

My two cents,
Tom

*Sadly, several people in ITK would like not to rely on eigen 3
because of its license:
http://www.itk.org/mailman/private/insight-developers/2010-March/014318.html
Right.  The eigendecomposition approach is part of the "preservation of
principal directions" algorithm used when reorienting diffusion tensors.
 The idea is that, as long as your deformation is a diffeomorphism, you can
estimate it, locally, as an affine transformation.   Then you use the rigid
part of the affine transformation to reorient the tensor's eigenvectors.
 There are a few other steps, though.
My question, in part, is whether ITK's pixel types are specific enough to
A diffusion tensor may be reoriented by either PPD (above) or finite strain
estimates.   So this is a non-unique mapping from pixel type to
reorientation strategy.
When reorienting a vector, one may want to change its relative orientation
but not change its length, even under an affine transformation.  Also a
non-unique mapping from pixel type to reorientation strategy.
One approach (what we've done in the past) is to implement filters for each
of these strategies for DTI.  While having filters would be nice, I can see
this becoming confusing, in particular if some metrics/interpolators are
performing reorientation as well.
Matt, I look forward to seeing the strain filter.  Thans much!
Re: Tcon , when would you like to discuss?  Previously, we used 10 am
Tuesday for registration discussions.
Brian
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:14 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details
And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?
Matt
Yes,
That sound very reasonable to me.
Luis
_______________________________________________
Insight-registration mailing list
http://public.kitware.com/cgi-bin/mailman/listinfo/insight-registration
Francois Budin
2010-11-22 15:47:36 UTC
Permalink
Hi everybody,

For tensor resampling, I had written a whole set of filters that could
be reused. They are available in this Insight Journal paper:
An ITK Implementation of a Diffusion Tensor Images Resampling Filter
http://hdl.handle.net/10380/3189
The resampling filter reorients the tensors. I guess if you actually
don't want to reorient the tensors, you can use the basic
itk::ResampleImageFilter.
Sincerely,

Francois

On Mon, Nov 22, 2010 at 3:41 AM, Tom Vercauteren
Post by Tom Vercauteren
Hi all,
As Brian just mentioned, there may be several ways of performing
re-orientation for a given type. PPD or FS being the most widely
encountered choices for diffusion tensors reorientation. Besides, what
can be of interest, especially for image registration, is the the
derivative of this reorientation strategy. It should indeed come into
the computation of the derivative of the objective function.
[Shameless add] For the specific case of FS in dimension 3, the
analytical derivative can be found in Thomas Yeo's paper (which I
http://people.csail.mit.edu/ythomas/publications/2009DTI-TMI.pdf
Hence, I believe that we need a generic approach where one can plug
its own reorientation strategy. In that view, I would prefer a 2bis
option where resampling and reorientation are done in one pass but
where the user can choose its reorientation strategy much like he can
now choose its interpolator.
The tricky part is how to handle the derivative of the reorientation
strategy in a generic manner in the registration framework and not
kill performance.
If we add a good linear algebra package (say eigen 3*), we could have
a no-op reorientation strategy that simply multiplies by identity
matrices but does not cost anything at run-time thanks to expression
templates.
My two cents,
Tom
*Sadly, several people in ITK would like not to rely on eigen 3
http://www.itk.org/mailman/private/insight-developers/2010-March/014318.html
Right.  The eigendecomposition approach is part of the "preservation of
principal directions" algorithm used when reorienting diffusion tensors.
 The idea is that, as long as your deformation is a diffeomorphism, you can
estimate it, locally, as an affine transformation.   Then you use the rigid
part of the affine transformation to reorient the tensor's eigenvectors.
 There are a few other steps, though.
My question, in part, is whether ITK's pixel types are specific enough to
A diffusion tensor may be reoriented by either PPD (above) or finite strain
estimates.   So this is a non-unique mapping from pixel type to
reorientation strategy.
When reorienting a vector, one may want to change its relative orientation
but not change its length, even under an affine transformation.  Also a
non-unique mapping from pixel type to reorientation strategy.
One approach (what we've done in the past) is to implement filters for each
of these strategies for DTI.  While having filters would be nice, I can see
this becoming confusing, in particular if some metrics/interpolators are
performing reorientation as well.
Matt, I look forward to seeing the strain filter.  Thans much!
Re: Tcon , when would you like to discuss?  Previously, we used 10 am
Tuesday for registration discussions.
Brian
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:14 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details
And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?
Matt
Yes,
That sound very reasonable to me.
Luis
_______________________________________________
Insight-registration mailing list
http://public.kitware.com/cgi-bin/mailman/listinfo/insight-registration
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://kitware.com/products/protraining.html
http://www.itk.org/Wiki/ITK_FAQ
http://www.itk.org/mailman/listinfo/insight-developers
Francois Budin
2010-11-22 16:41:14 UTC
Permalink
Hi everybody,

For tensor resampling, I had written a whole set of filters that could
be reused. They are available in this Insight Journal paper:
An ITK Implementation of a Diffusion Tensor Images Resampling Filter
http://hdl.handle.net/10380/3189
The resampling filter reorients the tensors. I guess if you actually
don't want to reorient the tensors, you can use the basic
itk::ResampleImageFilter.

Sincerely,
Francois
Post by Francois Budin
On Mon, Nov 22, 2010 at 3:41 AM, Tom Vercauteren
Post by Tom Vercauteren
Hi all,
As Brian just mentioned, there may be several ways of performing
re-orientation for a given type. PPD or FS being the most widely
encountered choices for diffusion tensors reorientation. Besides, what
can be of interest, especially for image registration, is the the
derivative of this reorientation strategy. It should indeed come into
the computation of the derivative of the objective function.
[Shameless add] For the specific case of FS in dimension 3, the
analytical derivative can be found in Thomas Yeo's paper (which I
http://people.csail.mit.edu/ythomas/publications/2009DTI-TMI.pdf
Hence, I believe that we need a generic approach where one can plug
its own reorientation strategy. In that view, I would prefer a 2bis
option where resampling and reorientation are done in one pass but
where the user can choose its reorientation strategy much like he can
now choose its interpolator.
The tricky part is how to handle the derivative of the reorientation
strategy in a generic manner in the registration framework and not
kill performance.
If we add a good linear algebra package (say eigen 3*), we could have
a no-op reorientation strategy that simply multiplies by identity
matrices but does not cost anything at run-time thanks to expression
templates.
My two cents,
Tom
*Sadly, several people in ITK would like not to rely on eigen 3
http://www.itk.org/mailman/private/insight-developers/2010-March/014318.html
Right.  The eigendecomposition approach is part of the "preservation of
principal directions" algorithm used when reorienting diffusion tensors.
 The idea is that, as long as your deformation is a diffeomorphism, you can
estimate it, locally, as an affine transformation.   Then you use the rigid
part of the affine transformation to reorient the tensor's eigenvectors.
 There are a few other steps, though.
My question, in part, is whether ITK's pixel types are specific enough to
A diffusion tensor may be reoriented by either PPD (above) or finite strain
estimates.   So this is a non-unique mapping from pixel type to
reorientation strategy.
When reorienting a vector, one may want to change its relative orientation
but not change its length, even under an affine transformation.  Also a
non-unique mapping from pixel type to reorientation strategy.
One approach (what we've done in the past) is to implement filters for each
of these strategies for DTI.  While having filters would be nice, I can see
this becoming confusing, in particular if some metrics/interpolators are
performing reorientation as well.
Matt, I look forward to seeing the strain filter.  Thans much!
Re: Tcon , when would you like to discuss?  Previously, we used 10 am
Tuesday for registration discussions.
Brian
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:14 PM, Matthew McCormick (thewtex) <
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
On Sat, Nov 20, 2010 at 7:00 PM, Matthew McCormick (thewtex)
Post by Matthew McCormick (thewtex)
Post by Luis Ibanez
(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.
I have a StrainImageFilter that calculates the strain tensor from a
deformation field if it would be of any help. I assume what is meant
by the 'Rotational and Stress component' are the same as the angles to
the principle strains and the principle strains. I intend on
submitting this filter to the Insight Journal in the next week or two
http://review.source.kitware.com/#q,status:open+project:ITK+branch:master+topic:BUG0010725_VTKTensors3,n,z
Post by Luis Ibanez
Post by Matthew McCormick (thewtex)
Matt
Excellent !
Yes,
that what I meant by the stress component.
Do you see a way of deriving as well
the rotational component of the Jacobian ?
I was thinking just using
http://www.itk.org/Doxygen320/html/classitk_1_1SymmetricEigenAnalysis.html#_details
And the eigenvectors determine the rotation, while the eigenvalues
(principle strains) determine the compression/expansion of the
components? Correct?
Matt
Yes,
That sound very reasonable to me.
Luis
_______________________________________________
Insight-registration mailing list
http://public.kitware.com/cgi-bin/mailman/listinfo/insight-registration
_______________________________________________
Powered by www.kitware.com
Visit other Kitware open-source projects at
http://www.kitware.com/opensource/opensource.html
http://kitware.com/products/protraining.html
http://www.itk.org/Wiki/ITK_FAQ
http://www.itk.org/mailman/listinfo/insight-developers
brian avants
2010-11-22 22:11:15 UTC
Permalink
I think Tom's idea is quite nice. We add a

resampler->SetReorientationStrategy( reorienter );

function to the ResampleImageFilter .... obviously, for most cases, it
would not be needed.

Tom, you mentioned "what can be of interest, especially for image
registration, is the the derivative of this reorientation strategy. It
should indeed come into the computation of the derivative of the
objective function." and also your FS paper.

I am not sure that we can implement the derivative wrt reorientation
in a general way, though it's possible to implement for some cases,
e.g. FS and PPD. Did you have more specific thoughts about this?


Also, Luis, should we schedule the tcon for tomorrow?

Brian
Tom Vercauteren
2010-11-23 07:41:31 UTC
Permalink
Hi Brian,
I think Tom's idea is quite nice.   We add a
resampler->SetReorientationStrategy( reorienter );
function to the ResampleImageFilter .... obviously, for most cases, it
would not be needed.
By the way, would there be any more appropriate name than
ReorientationStrategy? In some specific cases, I guess that the user
might be interested in not only rotating but also re-scaling the
diffusion tensor. That is, what might be wanted is to apply the full
Jacobian matrix: J^-1 . T . J

In this case, reorientation is not exactly what is done. Anyhow, with
some documentation, this might be clear enough.
Tom, you mentioned "what can be of interest, especially for image
registration, is the the derivative of this reorientation strategy. It
should indeed come into the computation of the derivative of the
objective function."  and also your FS paper.
I am not sure that we can implement the derivative wrt reorientation
in a general way, though it's possible to implement for some cases,
e.g. FS and PPD.   Did you have more specific thoughts about this?
Not really, I just thought it would be nice to support it if we can
(much like some computation of a pseudo-Hessian on objective
functions).

If no analytical derivation can be made, we could either return an
identity matrix, compute if from finite difference or automatic
differentiation (
http://en.wikipedia.org/wiki/Automatic_differentiation ) or provide a
way to inform the user that this function is not implemented. For
performance reasons, I wouldn't rely on exceptions but maybe on
something close to boost optional:
http://www.boost.org/doc/libs/release/libs/optional/index.html
http://www.boost.org/doc/libs/release/libs/optional/doc/html/boost_optional/examples.html#boost_optional.examples.optional_return_values


Tom
Luis Ibanez
2010-11-23 15:33:10 UTC
Permalink
Hi Tom,
Hi Brian,
Post by brian avants
I think Tom's idea is quite nice. We add a
resampler->SetReorientationStrategy( reorienter );
function to the ResampleImageFilter .... obviously, for most cases, it
would not be needed.
By the way, would there be any more appropriate name than
ReorientationStrategy? In some specific cases, I guess that the user
might be interested in not only rotating but also re-scaling the
diffusion tensor. That is, what might be wanted is to apply the full
Jacobian matrix: J^-1 . T . J
In this case, reorientation is not exactly what is done. Anyhow, with
some documentation, this might be clear enough.
How about:

"PixelTransformationFunction"
Post by brian avants
Tom, you mentioned "what can be of interest, especially for image
registration, is the the derivative of this reorientation strategy. It
should indeed come into the computation of the derivative of the
objective function." and also your FS paper.
I am not sure that we can implement the derivative wrt reorientation
in a general way, though it's possible to implement for some cases,
e.g. FS and PPD. Did you have more specific thoughts about this?
Not really, I just thought it would be nice to support it if we can
(much like some computation of a pseudo-Hessian on objective
functions).
If no analytical derivation can be made, we could either return an
identity matrix, compute if from finite difference or automatic
differentiation (
http://en.wikipedia.org/wiki/Automatic_differentiation ) or provide a
way to inform the user that this function is not implemented. For
performance reasons, I wouldn't rely on exceptions but maybe on
http://www.boost.org/doc/libs/release/libs/optional/index.html
http://www.boost.org/doc/libs/release/libs/optional/doc/html/boost_optional/examples.html#boost_optional.examples.optional_return_values
Tom
We may have to profile the computation time, to be sure,
but, looking at the code of the itkResampleImageFilter.txx,
we already perform enough computation per pixel, that
the additional call to a function per pixel may not have a
significant impact in the computation time.

If that's the case, we can simply have the pixel-rotation-helper
class be set by default to a class that performs a null-operation,
and that can be replaced, by the user, with any of the other
flavors of pixel rotation(scaling) that we may want to made
available.



Luis

Luis Ibanez
2010-11-23 15:24:29 UTC
Permalink
Hi Brian,
Post by brian avants
Also, Luis, should we schedule the tcon for tomorrow?
Brian
Hi Brian,

I'm reading this email too late.

Would you like to try a tcon at 1pm EST ?

or maybe tomorrow Wednesday at 10 am EST ?

(we have the SimpleITK tcon from 9:30am to 10am).


Thanks


Luis
Loading...