How does the inverseWarp work?

[devDonnn]

Hi. I'm trying to understand how ANTs internally manages to apply the inverse warp transformation. Say, I run SyN of two images: fixed and moving. I get the tx_0GenericAffine.mat, tx_1Warp.nii.gz and tx_1InverseWarp.nii.gz.

Now, both images - fixed and moving - are in different physical space (their RAS coordinates are different as well as dimensions). I want to apply tx_1InverseWarp.nii.gz and tx_0GenericAffine.mat to map the fixed image to the moving image space. I would do so by:

antsApplyTransforms -d 3 -r moving.nii.gz -i fixed.nii.gz -o fixed_warped.nii.gz -t [ tx_0GenericAffine.mat,1 ] -t tx_1InverseWarp.nii.gz

I found that RAS[fixed] = inverse_affine(inverse_warp(RAS[moving])). What I don't understand is how ANTs knows how to apply inverse warps to the RAS[moving] if inverse warps are in the fixed space and RAS[moving] corresponds to another space? I mean there is no voxel-wise correspondence between moving and 1InverseWarp.nii.gz neither their RAS coordinates are the same that it could interpolate. I feel like I'm missing something here...

[cookpa]

The warps are applied to point sets in physical space. The header transforms are only used when looking up intensities at the point's final location. For the purposes of transforms, it's all points in physical space.

The most confusing aspect is that the transformation of points is the opposite to how the "moving" image is transformed. Thus, if applying the forward warps, we start with a point set on a regular grid in the fixed space. For a point p in the center of a voxel v in the fixed space, we apply the warp and then the affine transform to get to a continuous point in moving space. There, we interpolate a value that becomes the value of v in the output image. If the Warp is W and the affine A, we do A(W(p)).

For inverse warps, we start with a point q at the center of voxel x in moving space, and need an interpolated value from the fixed space. We first apply the inverse affine A' to get approximately into the fixed space, and then find the interpolated inverse warp vector from W'. In other words, we do W'(A'(q)) to find the final position to interpolate for x in the output image.

At least, that's my reading of the transform logic

https://itk.org/Doxygen/html/classitk_1_1CompositeTransform.html

[gdevenyi]

The's a good breakdown in SimpleITK's fundamental concepts, https://simpleitk.readthedocs.io/en/release/fundamentalConcepts.html

[devDonnn]

Hmmm... Even more confusing. Let's consider the forward warps then. Based on the documentation, it says that first the affine and then the warp is applied: W(A(p)). In the forward case I was thinking that in order to get a moving image warped to the fixed image, we take the grid points in the fixed space, we apply the affine registration which brings us to the moving space. If we interpolate at this point, we get an affine registered moving to the fixed image.

To my understanding, since 1Warp.nii.gz is in the same space as the fixed image, it directly tells how each voxel of the fixed image is displaced, no?

[cookpa]

The 1Warp.nii.gz is the displacement of voxel centers in the fixed image towards the moving image for the SyN stage of registration, yes. If you just ran SyN on two images without any affine component, you'd get 0Warp.nii.gz which would be the entire transform.

But typically, registration is a multi-stage process, with the affine part computed first (or more usually, multiple affine / rigid / translation steps merged into one final 0GenericAffine.mat). At the SyN stage, the affine transform is done, so effectively, for fixed image F and moving image M, we're trying to optimize the similarity of F and A(M) during the SyN stage.

But in terms of how the resampling is actually done, the images are "pulled" into the output space. For forward warps, we start at p in the fixed space and find A(W(p)), interpolate a value, and so on to create the deformed moving image.

[devDonnn]

"<...> we're trying to optimize the similarity of F and A(M) during the SyN stage." - yes, this part makes sense. The obtained forward warps then directly tells how each voxel of the fixed image is displaced to the moving image, i.e. we find locations in the moving image which map to the fixed image. We only need to interpolate the moving image at these locations.

However, I don't understand what 1InverseWarp.nii.gz represents. If we now take the moving image as the reference space, how does it know how these moving image locations shall be displaced if 1InverseWarp and moving.nii.gz are in different spaces.

[cookpa]

1InverseWarp doesn't know how to displace points from moving space, only from A'(moving space)

[devDonnn]

I don't understand what goes first then. Consider this:
antsApplyTransforms -d 3 -r moving.nii.gz -i fixed.nii.gz -o fixedWarped.nii.gz -t [ 0GenericAffine.mat,1 ] -t 1InverseWarp.nii.gz

According to the documentation, it first applies the warps - moving voxel locations are displaced, and then affine-transformed to the fixed image, and finally interpolated. But now I see that it is not the case?

[cookpa]

Can you please link me to the documentation you're looking at?

[devDonnn]

My goal is to register the moving image to the fixed image and then apply the warps (inverse probably) to the fixed image only without the affine registration. I only want to see how the fixed image has displaced in its own space without reslicing to the moving space.

Does it make any sense to
antsApplyTransforms -d 3 -r fixed.nii.gz -i fixed.nii.gz -t tx_1InverseWarp.nii.gz -o fixedWarped.nii.gz ?

[cookpa]

Yes, you can do that to see how the fixed image is displaced by the warp without the affine component.

[devDonnn]

Ah, I see. I think I do understand now what's going 'under the hood'. The only concern I have now is with the following (just for curiosity):
antsApplyTransforms -d 3 fixed.nii.gz -i moving.nii.gz -t 1Warp.nii.gz

I only apply the warps to the points in the fixed space without the affine registration afterwards. The resulting image looks a bit weird and skewed but still fairly well interpolates the intensities of the moving image. I wonder what exactly happens at this particular case? We get W(p) but how does it relate to the moving image then without the affine transformation?

[cookpa]

Its relationship to the moving image is going to depend on the nature of the affine transform. If the affine transform is close to identity, then the results of dropping it will be small. If it's larger, then the sample points will be further off.

The size of the offset of a point from the affine transform A is going to depend on the correspondence of the images in physical space before registration. Sometimes this is very close (eg, registering T2w to T1w from the same imaging session).

[devDonnn]

I see, it all makes sense now. Thanks a lot, @cookpa!