BUG?: MGHHeader.set_zooms() does not persist on I/O roundtrip or affect affine

[larsoner]

I am trying to set the affine transform in order to scale an MRI along three dimensions (on latest master). I decided to try to do this by setting the zooms, but the scaling does not seem to persist on I/O round-trip, at least based on how I'm checking.

Reproducing code:

import os
import numpy as np
import nibabel

fname = os.path.join(os.getenv('FREESURFER_HOME'), 'subjects', 'fsaverage', 'mri', 'T1.mgz')

img_mgz = nibabel.load(fname)
nibabel.save(img_mgz, 'test.nii')
img_nii = nibabel.load('test.nii')
zooms = (0.5, 0.5, 0.5)
affine_zoomed = img_mgz.header.get_affine().copy()
affine_zoomed[:3] *= 0.5  # this is what I expect to persist
img_mgz.header.set_zooms(zooms)
img_nii.header.set_zooms(zooms)
# good (3)
assert np.allclose(img_mgz.header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(img_nii.header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(img_mgz.header.get_affine(), affine_zoomed, atol=1e-5)
# bad (1)
assert not np.allclose(img_nii.header.get_best_affine(), affine_zoomed, atol=1e-5)

#
# Do some I/O
#

nibabel.save(img_nii, 'test_nii.nii')
nibabel.save(img_mgz, 'test_mgz.nii')
nibabel.save(img_nii, 'test_nii.mgz')
nibabel.save(img_mgz, 'test_mgz.mgz')

# good (1)
assert np.allclose(nibabel.load('test_nii.nii').header.get_zooms(), zooms, atol=1e-5)
# bad (5)
assert not np.allclose(nibabel.load('test_mgz.nii').header.get_zooms(), zooms, atol=1e-5)
assert not np.allclose(nibabel.load('test_nii.mgz').header.get_zooms(), zooms, atol=1e-5)
assert not np.allclose(nibabel.load('test_mgz.mgz').header.get_zooms(), zooms, atol=1e-5)
# changed due to I/O steps!
assert not np.allclose(img_mgz.header.get_affine(), affine_zoomed, atol=1e-5)
assert not np.allclose(img_mgz.header.get_zooms(), zooms, atol=1e-5)

It looks like with this code 1) there are consistency problems, and 2) problems with img_mgz.header.get_zooms() being affected / reset by the call(s) to nibabel.save.

Is there a different way I should rescale the MRI? Am I doing something silly here?

[larsoner]

It looks like the problems are caused by MGHImage._affine, as eventually this leads to a call to _affine2header that overwrites the .header values.

If I do this:

img_mgz._affine = img_mgz.header.get_affine()

or

img_mgz._affine = None

Then all of the MGZ-format round-trip bugs go away (the .nii ones persist). I'm a bit out of my depth to know why this was there originally, though.

[effigies]

Assigning myself to remind me to look into this. I don't think I'll be able to get to this before a release, so someone else is welcome to have a go, as well.

[ellisdg]

I ran into this issue as well. After looking into it a bit, it seems that when saving an image, nibabel looks at the affine from the image, image._affine, and then updates the header according to the information stored in the affine. When you set the zooms in the header, it doesn't actually update image._affine. Thus, when you try and save the image to file, the adjusted zooms just get overwritten by the zooms as stored in image._affine. I'm guessing that for MGZ-format images, setting

img_mgz._affine = None

forces the image to look at the header and compute the adjusted affine with the adjusted zooms before saving.

Due to the way Images and Headers are designed, it's not possible for a method in the Header to update the affine in the Image, right? Also, it seems like there may be a bug in that the sfrom and qform matrices are not updated when the zooms are set in the AnalyzeHeader.

For my use case I was able to write functions that change the affine directly:

import os
import numpy as np
import nibabel


def get_zooms_from_image(image):
    return get_zooms_from_affine(image.affine)


def get_zooms_from_affine(affine):
    RZS = affine[:3, :3]
    return np.sqrt(np.sum(RZS * RZS, axis=0))


def set_image_zooms(image, spacing):
    scale = np.divide(spacing, get_zooms_from_image(image))
    image.affine[:3, :3] *= scale
    image._affine2header()


fname = os.path.join(os.getenv('FREESURFER_HOME'), 'subjects', 'fsaverage', 'mri', 'T1.mgz')

img_mgz = nibabel.load(fname)
nibabel.save(img_mgz, 'test.nii')
img_nii = nibabel.load('test.nii')

zooms = (0.5, 0.5, 0.5)
set_image_zooms(img_mgz, zooms)
set_image_zooms(img_nii, zooms)

assert np.allclose(img_mgz.header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(img_nii.header.get_zooms(), zooms, atol=1e-5)

#
# Do some I/O
#

nibabel.save(img_nii, 'test_nii.nii')
nibabel.save(img_mgz, 'test_mgz.nii')
nibabel.save(img_nii, 'test_nii.mgz')
nibabel.save(img_mgz, 'test_mgz.mgz')

assert np.allclose(nibabel.load('test_nii.nii').header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(nibabel.load('test_mgz.nii').header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(nibabel.load('test_nii.mgz').header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(nibabel.load('test_mgz.mgz').header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(img_mgz.header.get_zooms(), zooms, atol=1e-5)

These functions don't adjust the origin of the affine when you adjust the zooms. I'm not sure how the zooms should affect the origin.

[effigies]

Hi @larsoner @ellisdg, thanks for this conversation, and sorry for being so slow in getting back to it. On the plus side, I don't think it's a bug, per se.

Now that I think about it, I'm not sure of the usefulness of set_zooms(), except for adjusting the TR. In NIFTI, it allows you to decouple the pixdim values from the sform affine, but I would probably consider a file where the zooms and affine do not match to be invalid.

I think the correct way would be to update your affine directly, as @ellisdg suggests.

I'm not sure how the zooms should affect the origin.

Assuming you want the RAS origin to reside within the same IJK voxel, scale the translation parameters by the same factors as the rotation/zoom/shear matrix.

Given an affine, decompose into M : 3x3 and B : 3x1:

RAS = A * IJK
RAS = M * IJK + B

The RAS origin (0, 0, 0) has IJK coordinate:

0 = M * IJK - B
IJK = -M^-1 * B

If we want the same IJK to point to 0 for M' as M, then we need a B' that satisfies:

-M'^-1 * B' = -M^-1 * B
B' = M' * M^-1 * B

So, in general, new_vec = new_mat.dot(np.linalg.inv(orig_mat)).dot(orig_vec).

If M = R * Z * S and M' = R * (0.5 * Z) * S, then M' = 0.5 * M, so this simplifies trivially:

B' = 0.5 * M * M^-1 * B
   = 0.5 * B

And if your affine has no shears, applying non-uniform scaling will be easy enough to do without inverting your matrix.

[larsoner]

I don't think it's a bug, per se.

This suggests that the minimal example I pasted above is doing the right thing and what users should expect. In this case that means that there can be some Image-class (NIfTI vs MGHImage) and I/O-dependent (round-trip failure) behavior.

I think the correct way would be to update your affine directly, as @ellisdg suggests.

Looking at that solution, it requires using the private _affine2header function. Without that line, these two lines fail in that example script:

assert np.allclose(img_mgz.header.get_zooms(), zooms, atol=1e-5)
assert np.allclose(img_nii.header.get_zooms(), zooms, atol=1e-5)

Changing private attributes and using private methods does not seem like a good suggested approach (though it's a fine workaround!). Fortunately, it looks like you can get away with using the public image.header.set_zooms(spacing) in place of the private image._affine2header() call in that script, though.

Now that I think about it, I'm not sure of the usefulness of set_zooms(), except for adjusting the TR

Naively I expected that changing the first three elements would change the affine. This line in my minimal example passes, suggesting that at least for MGHImage it does indeed change:

assert np.allclose(img_mgz.header.get_affine(), affine_zoomed, atol=1e-5)

If changing the zooms should not affect the affine and zooms[:3] can be overwritten by I/O round-trip, then this should probably be documented. Ideally the behavior or MGHImage changed not to change the affine for consistency with other classes, too.

But to avoid changing any code, this could also be fixed by documentation (probably in header.set_zooms) stating that if you intend to change zooms[:3], you should both modify image.affine and call image.header.set_zooms to ensure consistency across classes and I/O round-trips.

[effigies]

Ah, this isn't to say that I would recommend that exact approach. I agree that users shouldn't be expected or encouraged to use _affine2header(). Perhaps I should have said that "Using the target affine is the way to go", as what I would actually recommend is creating a new image with an updated affine:

orig = nb.load(fname)
affine = scale_affine(orig.affine, new_zooms)  # left as an exercise to the reader
img_nii = nb.Nifti1Image(orig.dataobj, affine)
img_mgz = nb.MGHImage(orig.dataobj, affine)

But you do raise some good points.

In this case that means that there can be some Image-class (NIfTI vs MGHImage) and I/O-dependent (round-trip failure) behavior.

There shouldn't be image-class differences, here. I would say, in this case, that it's actually a failure of Nifti1Image to update pixdim* to match the affine in update_header().

My understanding is that SpatialImage interface treats img.dataobj and img.affine as authoritative. img.header is a bunch of format-specific metadata, but to the extent that affine information is encoded there, it should be overridden by img.affine when img.update_header() is called. To set a new affine, you actually should create a new image.

It might make this more obvious to the user if img.affine were either marked as not writable or returned copies such that trying to change it would have no effect on the image.

Ideally the behavior or MGHImage changed not to change the affine for consistency with other classes, too.

The problem here is that the affine isn't stored in the MGH header; it is calculated from the zooms and the directional cosine matrix. (Similar to NIfTI's qform, actually. If you set_zooms() and then get_qform(), you'll see an updated affine matrix.)

Nonetheless, I would basically consider these file format quirks, where you can get the header out of sync with the affine matrix in somewhat different ways. We should mostly just take care of making sure they're in sync at creation/loading and saving. Perhaps we should add a warning when they are out-of-sync at save?

Does that all seem reasonable? I think this leaves us with the following potential fixes/updates:

1. Update pixdim[1..3] in Nifti1Image.update_header() (or the appropriate superclass)
2. Warn if the header's affine information is being changed when writing the image out
3. Make img.affine return a copy of the affine OR
4. Set img.affine.writable = False
5. Function to update affine with new zooms
6. Clearly document in SpatialImage that affine-related header info may be clobbered during write

Let me know what you think.

[larsoner]

To set a new affine, you actually should create a new image.

I didn't realize this was the preferred approach. Doesn't this mean anything stored in header other than the affine is lost, including things like the TR (which would live in zooms[3] and not in affine)? I assume there are other important fields that would be lost, too. So this is a bit surprising to me.

Perhaps we should add a warning when they are out-of-sync at save?

This would at least make it more obvious that potential problems exist, yeah.

I think this leaves us with the following potential fixes/updates:

I don't mess with MRI file formats (or data!) enough to have a good informed opinion here. But it sounds like (1), (2), (5), and (6) would all be useful. (3) and (4) sound risky because they could break existing code that relies on the current in-place-modification behaviors (4 would at least emit an error, but 3 could break things silently!), and I assume breakages would includes the workaround code that @ellisdg proposed (as it relies on image.affine[:3, :3] *= scale working properly).

[effigies]

Doesn't this mean anything stored in header other than the affine is lost, including things like the TR (which would live in zooms[3] and not in affine)?

Yes, though you can pass in the header as well. It just won't necessarily go smoothly in the general case where you don't know whether the image class you're targeting has the same header format as your source image.

I often use this idiom:

img = orig.__class__(data, affine, orig.header)  # Usually affine is just orig.affine, too

I think that if you pass an MGHHeader to Nifti1Image, it'll be about the same as passing no header, but I didn't try it, so didn't want to post it.

Anyway, thanks for the feedback. I'll try to fork these out to separate issues in the next day or so, so that others can weigh in or try their hand at addressing them.