New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Function to reorient gradient directions according to moco parameters #669

Merged
merged 9 commits into from Dec 30, 2015

Conversation

Projects
None yet
6 participants
@arokem
Member

arokem commented Jun 24, 2015

Loosely based on @rfdougherty's implementation here: https://github.com/cni/nims/blob/master/scripts/pepolar_unwarp.py

TODO:

  • More thorough testing.
  • Allow inputs of shape (n, 4, 4) if users get this kind of output from their motion correction algorithm.
@@ -244,3 +244,55 @@ def gradient_table(bvals, bvecs=None, big_delta=None, small_delta=None,
small_delta=small_delta,
b0_threshold=b0_threshold,
atol=atol)
def reorient_bvecs(gtab, rots):

This comment has been minimized.

@omarocegueda

omarocegueda Jun 24, 2015

Contributor

Where is the user expected to obtain the rots input from? For example, they may have the rotation matrices themselves instead of Euler angles
Edit: sorry, I just found out that rfdougherty's implementation uses commands from FSL (topup and eddy).

@arokem

This comment has been minimized.

Member

arokem commented Jun 25, 2015

No worries. This is a good point for discussion! I think that we should at least allow two kinds of input:

  1. An (n, 3) with the x,y,z rotation angles associated with each volume.
  2. An (n, 4, 4) with the affine transformation associated with each volume.

Hopefully, we'll get some parameters from a motion-correction algorithm based on your current PR #654. I assume that the output from that would be more like the latter? Would maybe be good to write a little utility function that extracts the former from the latter.

----------
gtab : a GradientTable class instance.
The nominal gradient table with which the data were acquired.
rots : (n, 3) array.

This comment has been minimized.

@omarocegueda

omarocegueda Jun 25, 2015

Contributor

I am worried that some people may use this as a black box (I'm not sure if we should add a warning here or be more specific about what the rots are) because there may be an issue regarding how the rotational parameters are given by different software. For example, flirt (and mcflirt, and I assume all FSL software, like eddy) returns a matrix M that maps points in physical space from the moving image towards the static image (that is the most intuitive way for most people, because it is the moving image which is being "sent" to the static image). ANTS, on the other hand, returns a matrix that maps points in the static image towards points in the moving image, but this choice makes sense too: when we transform the moving image towards the static image, we actually send each point x in the static image towards its corresponding y in the moving image (contrary to the common intuition) so we can interpolate the moving image at y to obtain the intensity value that must be written at each x to form the transformed image. So, whenever Flirt needs to transform an image, it actually computes the inverse of the matrix it returned, while ants simply uses the returned matrix itself. However, if people used Flirt, they must apply the returned matrix unchanged to rotate the bvecs, but if for some reason they used ANTS they need to apply the inverse of the returned matrix. Another issue is whether the transformation was computed in RAS or LPS (or any other) coordinate system, and if the voxels in the image are given in C order of Fortran order. The dangerous thing is that the transforms are often very close to the identity, which makes it very hard to notice when we are making a mistake.

This comment has been minimized.

@MrBago

MrBago Jun 26, 2015

Contributor

My understanding is that rotation angles are not well defined, different rotation angles can be used to represent the same transformation and different conventions exist for different applications. For example Euler angles vs Tait-Bryan angles.

Is there a way to simplify this, for example would it be easier to use an affine matrix or a "transformation object" to help standardize this?

Subject Motion in DTI Data. Leemans, A. and Jones, D.K. (2009).
MRM, 61: 1336-1349
"""
new_bvecs = gtab.bvecs[~gtab.b0s_mask]

This comment has been minimized.

@MrBago

MrBago Jun 26, 2015

Contributor

Are you sure you wan to just apply this to only ~gtab.b0s_mask?

Either way, there should error checking here to make sure the number rotations matches the number of gradients. When the data has b0s mixed in with gradients (for example hc data) it's easy to occidentally pass in too many or too few rotation angles.

@arokem

This comment has been minimized.

Member

arokem commented Jun 29, 2015

Thanks for the comments!

@omarocegueda: I think that the most important would be to be internally consistent within dipy. That is, best would be to be compatible with your registration work, because I think that we will eventually want to use your work for motion correction as well (?). I see that you are discussing implementation of a Transformation class, over on #654. Maybe the input here should be a list of these objects? I saw that you mentioned a Transformation class that is implemented in SymmetricDiffeomorphicRegistration and in StreamlineLinearRegistration. Could you please point out to me where I could find the implementation of this Transformation class? Might be useful to think of this use-case as we implement that.

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jun 29, 2015

Agreed, @arokem. I think each software should provide its own method for rotating the b-vectors to avoid any confusion. However, I don't know much about the state of the art on registration-based motion correction for Diffusion: as far as I know, the easiest approach is simply register all DW volumes to the B0, but I'm not convinced that this is a good approach. I have been studying FSL's eddy code, which is based on a more sound approach by attempting to model the diffusion signal under the eddy-current and subject motion artifacts (it models the signal as a Gaussian process). I think we should go for something more in that direction, but for now, with the current framework it is very easy to implement the simple approach of registering all DW images to the B0 (maybe we can put that in the tutorial, like "Write your own motion correction method with b-vec rotation").

I saw that you mentioned a Transformation class that is implemented in SymmetricDiffeomorphicRegistration and in StreamlineLinearRegistration

Sure!, here is the one for streamlines:
https://github.com/nipy/dipy/blob/master/dipy/align/streamlinear.py#L395

And here is the one for images:
https://github.com/nipy/dipy/blob/master/dipy/align/imwarp.py#L403

Actually, @Garyfallidis, I think we can use one single AffineMap for images and streamlines, the only difference is that in the case of images it needs to provide the grid2world transform for each image, while for streamlines, it's not necessary because the streamlines are already in physical space (wouldn't it be nice to transform an image with the same transformation you got from streamline registration? =) ), actually, there's nothing that prevents us from defining a transform method for streamlines in DiffeomorphicMap as well (i.e., how the streamlines are warped by the displacement field?). Let me define the AffineMap with both capabilities and we'll see if it is convenient.

@arokem

This comment has been minimized.

Member

arokem commented Jun 29, 2015

On Mon, Jun 29, 2015 at 11:07 AM, Omar Ocegueda notifications@github.com
wrote:

Agreed, @arokem https://github.com/arokem. I think each software should
provide its own method for rotating the b-vectors to avoid any confusion.
However, I don't know much about the state of the art on registration-based
motion correction for Diffusion: as far as I know, the easiest approach is
simply register all DW volumes to the B0, but I'm not convinced that this
is a good approach. I have been studying FSL's eddy code, which is based
on a more sound approach by attempting to model the diffusion signal under
the eddy-current and subject motion artifacts (it models the signal as a
Gaussian process). I think we should go for something more in that
direction, but for now, with the current framework it is very easy to
implement the simple approach of registering all DW images to the B0 (maybe
we can put that in the tutorial, like "Write your own motion correction
method with b-vec toration").

I agree that more advanced methods would be good to have in the future. For
now, I would like to have this functionality (moco by registration to mean
b0) implemented and easy to use in Dipy. My colleagues and I have used this
approach on many data-sets in the past, and it seems to work well for a
lot
of scenarios. As you point out, the implementation of this approach
flows naturally from this current work.

Ideally, this would end up not as a tutorial, but an actual scripting API
to do this, possibly as part of a workflows module that we will have for
users at various levels to use in their scripting of common tasks. But, for
the time being, the tutorial you propose sounds like an excellent idea. I
think that it can be a separate PR, once this is merged.

I saw that you mentioned a Transformation class that is implemented in
SymmetricDiffeomorphicRegistration and in StreamlineLinearRegistration

Sure!, here is the one for streamlines:
https://github.com/nipy/dipy/blob/master/dipy/align/streamlinear.py#L395

And here is the one for Diffeomorphisms:
https://github.com/nipy/dipy/blob/master/dipy/align/imwarp.py#L403

Thanks.

Upon further thought and reading through these things, I see two options:

  1. Use a list of matrices that are the output of affreg.optimize in the
    example you have here.
  2. Implement a uniform API for affine transformations, that has transorm,
    transform_inverse and invert, similar to the API you have in imwarp.
    This could take as input one of the matrices mentioned in 1.

Actually, @Garyfallidis https://github.com/Garyfallidis, I think we can
use one single AffineMap for images and streamlines, the only difference
is that in the case of images it needs to provide the grid2world
transform for each image, while for streamlines, it's not necessary because
the streamlines are already in physical space (wouldn't it be nice to
transform an image with the same transformation you got from streamline
registration? =) ), actually, there's nothing that prevents us from
defining a transform method for streamlines in DiffeomorphicMap as well
(i.e., how the streamlines are warped by the displacement field?). Let me
define the AffineMap with both capabilities and we'll see if it is
convenient.


Reply to this email directly or view it on GitHub
#669 (comment).

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jun 29, 2015

y colleagues and I have used this approach on many data-sets in the past, and it seems to work well for a lot of scenarios

That's great @arokem!, by any chance can you share a dataset with noticeable motion that is successfully corrected with this approach? So far I have been using this year's ISMRM tractography challenge data set: http://www.tractometer.org/ismrm_2015_challenge/data (unfortunately it seems that it's going to take a long time to have access to the ground truth, though) but I have not been able to obtain a correction that looks as good as reported in FSL's eddy page: http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/EDDY

@matthew-brett

This comment has been minimized.

Member

matthew-brett commented Jun 29, 2015

Forgive me for not really following, but if you are going to go for a simple rotation parameter set, I agree it would be best to avoid Euler angles. I would personally go for 3x3 rotation matrices, but you could also go for a rotation axis and rotation angle.

@alexis-roche - any views?

http://matthew-brett.github.io/transforms3d/reference/transforms3d.euler.html

@jchoude

This comment has been minimized.

Contributor

jchoude commented Jul 1, 2015

@omarocegueda Just a quick note: contrarily to your beliefs, the ground truth for the Challenge are online: http://www.tractometer.org/ismrm_2015_challenge/data

@omarocegueda

This comment has been minimized.

Contributor

omarocegueda commented Jul 1, 2015

This is fantastic @jchoude! may be we can now quantitatively compare different motion correction approaches under controlled conditions, I'll take a look, thank you very much!

@arokem arokem force-pushed the arokem:reorient-bvecs branch from d7fb62b to 0916a22 Dec 10, 2015

@arokem arokem changed the title from WIP: Function to reorient gradient directions according to moco parameters to Function to reorient gradient directions according to moco parameters Dec 10, 2015

@arokem

This comment has been minimized.

Member

arokem commented Dec 10, 2015

I finally need this myself, so I've revisited this and made fixes according to the comments I got in the previous round. In particular, this now assumes a list or array of affine 4 by 4 transformation matrices as input, and tests for a match between the input and the number of non-zero gradients in the input gradient table.

@arokem

This comment has been minimized.

Member

arokem commented Dec 11, 2015

Apparently the polar transformation was only recently added to scipy, so I added this in a fixes file that I called scipy. Is this the right way to support use in older scipy?

Parameters
----------
gtab : a GradientTable class instance.

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member

GradientTable is sufficient no need for a...

----------
gtab : a GradientTable class instance.
The nominal gradient table with which the data were acquired.
affines : list or ndarray of shape (n, 4, 4) or (n, 3, 3)

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member
affines : list of ndarrays of shape (4, 4) or (3, 3) or ndarray of shape (n, 4, 4) or (n, 3, 3)

This comment has been minimized.

@arokem

arokem Dec 29, 2015

Member

That doesn't fit within 80 characters.

new_gt = reorient_bvecs(gt, affs)
npt.assert_equal(gt.bvecs, new_gt.bvecs)
# Now apply some rotations to these suckers:

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member

language please...

rotated_bvecs[i] = np.dot(rotation_affines[-1][:3, :3],
bvecs[i])
# Copy over the rotation affines:

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member

no need for : in the comments.

def reorient_bvecs(gtab, affines):
"""
Reorient the directions in a GradientTable according

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member

Reduce docstring size by starting from the first line of the docstring

Parameters
----------
a : (m, n) array_like

This comment has been minimized.

@Garyfallidis

Garyfallidis Dec 29, 2015

Member

We usually write the shape after the object like in here
http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.spatial.Delaunay.html
It would be nice if we could follow a standard for this.

This comment has been minimized.

@arokem

arokem Dec 29, 2015

Member

This is copied straight from the scipy code-base. But will do.

@arokem arokem force-pushed the arokem:reorient-bvecs branch from 51a604d to ca154fa Dec 29, 2015

@arokem

This comment has been minimized.

Member

arokem commented Dec 29, 2015

Rebased, and addressed all the comments by @Garyfallidis.

Garyfallidis added a commit that referenced this pull request Dec 30, 2015

Merge pull request #669 from arokem/reorient-bvecs
Function to reorient gradient directions according to moco parameters

@Garyfallidis Garyfallidis merged commit 161b9f7 into nipy:master Dec 30, 2015

1 check passed

continuous-integration/travis-ci/pr The Travis CI build passed
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment