Skip to content
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

applying affine is creating some distortion or an illusion #1076

Open
satra opened this issue Jan 28, 2022 · 5 comments
Open

applying affine is creating some distortion or an illusion #1076

satra opened this issue Jan 28, 2022 · 5 comments

Comments

@satra
Copy link
Member

satra commented Jan 28, 2022

i'm taking an image and simply changing the zoom with this:

epi_2mm_nii = resample_to_output(epi0, (2,2,2))

and then plotting the before and after with the niworkflows display module. this results in this image:

https://dl.dropbox.com/s/tvx6c1wsbatdhgv/test.svg?dl=0 (couldn't figure out a way to display this in github - you will have to download and open with a browser).

i wasn't expecting any visual differences, so now i'm contemplating if this is just an illusion.

@satra
Copy link
Member Author

satra commented Jan 28, 2022

orig affine

[[  3.5940001    0.           0.         -65.4303894 ]
 [  0.           3.5940001    0.         -96.33889008]
 [  0.           0.           4.         -86.82946777]
 [  0.           0.           0.           1.        ]]

new affine

[[  2.           0.           0.         -65.4303894 ]
 [  0.           2.           0.         -96.33889008]
 [  0.           0.           2.         -86.82946777]
 [  0.           0.           0.           1.        ]]

@effigies
Copy link
Member

@satra Did you figure anything out? I think I looked at this back then, but also am not sure if there's an illusion.

@romainVala
Copy link

Hi
it would say there is an induced shift due to the resampling. It is difficult to be sure because change in resolution also lead to different partial volume, and with high folding of cortical gray matter reslicing may give different shape (for a given slice).

If I compare with mrgrid resampling function of mrtrix I get different result for the affine for instance if I started with

3.600000  0.000000  0.000000  -122.649391
0.000000  3.600000  0.000000  -128.450073
0.000000  0.000000  4.000000  -96.661674
0.000000  0.000000  0.000000  1.000000

if I resample to 2 mm iso, resample_to_output lead to the exact same affine matrix (as shown by satra) but mrgrid is giving me a small difference in the transaltion :

 2.000000  0.000000  0.000000  -123.849396
 0.000000  2.000000  0.000000  -129.650070
 0.000000  0.000000  2.000000  -97.661674
 0.000000  0.000000  0.000000  1.000000

I am very confident with mrtrix handeling of the affine, but may be it is just due to a different convention. Is it related to the position of the corner of the voxel or of the center (as a reference ?... or I am wrong and this is pertinent only for the viewer ... ?) having a shift less than a voxel size is small and may be we do not care, but if one can be exact ...

Since I see difference I compared with a third tool (I tested Resample transform from trochio, which rely on sitk) and .... unfurtunatly I get a third result ...

2.000000  0.000000  0.000000  -123.449394
0.000000  2.000000  0.000000  -129.250076
0.000000  0.000000  2.000000  -97.661674
0.000000  0.000000  0.000000  1.000000

It also changing the affine (which I think make sense) but end up with a different one ... who's right ???

I also do this a clear difference of resample_to_output when you have a non strictly axial orientation

For instance if I start with this affine

-1.632297  0.065883  0.102517  114.033714
0.017452  1.494715  -0.652199  -74.317558
0.122627  0.664246  1.457438  -96.661674
0.000000  0.000000  0.000000  1.000000

and resample to 1mm iso, nibabel (resample_to_output) give me a strictly axial output

1.000000  0.000000  0.000000  -122.649391
0.000000  1.000000  0.000000  -128.450073
0.000000  0.000000  1.000000  -96.661674
0.000000  0.000000  0.000000  1.000000

whereas mrtrix gives

-0.997133  0.040246  0.064073  114.331818
0.010661  0.913087  -0.407624  -74.570763
0.074910  0.405773  0.910899  -96.905731
0.000000  0.000000  0.000000  1.000000

and torchio (sitk) give a similar one

 -0.997133  0.040246  0.064073  114.319260
 0.010661  0.913087  -0.407624  -74.489479
 0.074910  0.405773  0.910899  -97.088036
 0.000000  0.000000  0.000000  1.000000

For the last example I found nibable output weird: With a change in the voxel size I expect the global orientation to be unchanged (numerical values of the affine parameters change because of a different total number of voxel, but for both (mrtix and torchio) I get the same global orientation as the input volume

As for the strick axial case, I do see small differences in the translation part between sitk and mrtrix but at least the rotation part are exactly identical !

@effigies
Copy link
Member

effigies commented Apr 4, 2023

If we're creating a bad affine, it will be in:

def vox2out_vox(mapped_voxels, voxel_sizes=None):
"""output-aligned shape, affine for input implied by `mapped_voxels`
The input (voxel) space, and the affine mapping to output space, are given
in `mapped_voxels`.
The output space is implied by the affine, we don't need to know what that
is, we just return something with the same (implied) output space.
Our job is to work out another voxel space where the voxel array axes and
the output axes are aligned (top left 3 x 3 of affine is diagonal with all
positive entries) and which contains all the voxels of the implied input
image at their correct output space positions, once resampled into the
output voxel space.
Parameters
----------
mapped_voxels : object or length 2 sequence
If object, has attributes ``shape`` giving input voxel shape, and
``affine`` giving mapping of input voxels to output space. If length 2
sequence, elements are (shape, affine) with same meaning as above. The
affine is a (4, 4) array-like.
voxel_sizes : None or sequence
Gives the diagonal entries of `output_affine` (except the trailing 1
for the homogeneous coordinates) (``output_affine == np.diag(voxel_sizes
+ [1])``). If None, return identity `output_affine`.
Returns
-------
output_shape : sequence
Shape of output image that has voxel axes aligned to original image
output space axes, and encloses all the voxel data from the original
image implied by input shape.
output_affine : (4, 4) array
Affine of output image that has voxel axes aligned to the output axes
implied by input affine. Top-left 3 x 3 part of affine is diagonal with
all positive entries. The entries come from `voxel_sizes` if
specified, or are all 1. If the image is < 3D, then the missing
dimensions will have a 1 in the matching diagonal.
"""
try:
in_shape, in_affine = mapped_voxels.shape, mapped_voxels.affine
except AttributeError:
in_shape, in_affine = mapped_voxels
n_axes = len(in_shape)
if n_axes > 3:
raise ValueError('This function can only deal with 3D images')
if n_axes < 3:
in_shape += (1,) * (3 - n_axes)
out_vox = np.ones((3,))
if voxel_sizes is not None:
if not len(voxel_sizes) == n_axes:
raise ValueError('voxel sizes length should match shape')
if not np.all(np.array(voxel_sizes) > 0):
raise ValueError('voxel sizes should all be positive')
out_vox[:n_axes] = voxel_sizes
in_mn_mx = zip([0, 0, 0], np.array(in_shape) - 1)
in_corners = list(product(*in_mn_mx))
out_corners = apply_affine(in_affine, in_corners)
out_mn = out_corners.min(axis=0)
out_mx = out_corners.max(axis=0)
out_shape = np.ceil((out_mx - out_mn) / out_vox) + 1
out_affine = np.diag(list(out_vox) + [1])
out_affine[:3, 3] = out_mn
return tuple(int(i) for i in out_shape[:n_axes]), out_affine

I suspect the issue is that we're keeping constant the location of the corner of the first voxel rather than the center.

@romainVala
Copy link

romainVala commented Apr 4, 2023

Yes the small differences I pointed are indeed a matter of the exact chosen reference. This choice should have been defined in the format so that every one take the same ... but ok we have to live with imperfect standard definition ...

The important is to stay consistent and also to use the visualization software that goes with the same choice (corner versus center )

So we may forget those small differences in the translation part but for the rotation information I would say there is a bug in the nibabel implementation ... it construct a diagonal 3*3 matrix for the affine although it should start from the input affine.

I would build a zoom affine (identity with zooming factor in the diagonal) and then just divide or multiply with the input affine in order to get the output affine ...
I hope it make sense

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

3 participants