Reset of qform and sform codes on save

[matthew-brett]

Andrew Connolly pointed out the following misbehavior:

import numpy as np
import nibabel as nib

from nose.tools import assert_equal

img = nib.Nifti1Image(np.zeros((2,3,4)), None)
hdr = img.get_header()
hdr['qform_code'] = 3
hdr['sform_code'] = 4
nib.save(img, 'afile.nii')
img_back = nib.load('afile.nii')
hdr_back = img_back.get_header()
# These succeed
assert_equal(hdr_back['qform_code'], 3)
assert_equal(hdr_back['sform_code'], 4)

nib.save(img_back, 'afile2.nii')
# These fail - the save has modified the image in memory
assert_equal(hdr_back['qform_code'], 3)
assert_equal(hdr_back['sform_code'], 4)
img_back_again = nib.load('afile2.nii')
hdr_back_again = img_back_again.get_header()
# These fail - the saved image is also modified
assert_equal(hdr_back_again['qform_code'], 3)
assert_equal(hdr_back_again['sform_code'], 4)

[MrBago]

Hi Matthew,
Thanks for spending time on this.

It seems like you guys want to have set_sform and set_qform be private methods because calling those methods directly causes inconsistencies in the image and header affine.

In [1]: import numpy as np

In [2]: import nibabel as nib

In [3]: aff = np.eye(4)

In [4]: vol = np.zeros((5,5,5))

In [6]: img = nib.Nifti1Image(vol, aff)

In [7]: hdr = img.get_header()

In [8]: hdr.set_sform(-aff, 1)

In [9]: hdr.get_best_affine()
Out[9]: 
array([[-1., -0., -0., -0.],
       [-0., -1., -0., -0.],
       [-0., -0., -1., -0.],
       [ 0.,  0.,  0.,  1.]])

In [10]: img.update_header()

In [11]: hdr.get_best_affine()
Out[11]: 
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

I personally think any consistency checks should be done when the affine gets modified and not on write because the following seems really unintuitive to me and I think it'll be confusing to users.

In [18]: hdr.set_sform(-aff, 1)

In [19]: hdr.get_sform()
Out[13]: 
array([[-1., -0., -0., -0.],
       [-0., -1., -0., -0.],
       [-0., -0., -1., -0.],
       [ 0.,  0.,  0.,  1.]])

In [20]: nib.save(img, 'my_img.nii')

In [21]: hdr.get_sform()
Out[21]: 
array([[ 1.,  0.,  0.,  0.],
       [ 0.,  1.,  0.,  0.],
       [ 0.,  0.,  1.,  0.],
       [ 0.,  0.,  0.,  1.]])

I'm willing to spend some time on this, but it would help me to know how many ways the affine on an 'image' can get changed.

[matthew-brett]

Hi Bago,

On Mon, Dec 19, 2011 at 8:02 PM, MrBago
[email protected]
wrote:

Hi Matthew,
 Thanks for spending time on this.

It seems like you guys want to have set_sform and set_qform be private methods because calling those methods directly causes inconsistencies in the image and header affine.

In [1]: import numpy as np

In [2]: import nibabel as nib

In [3]: aff = np.eye(4)

In [4]: vol = np.zeros((5,5,5))

In [6]: img = nib.Nifti1Image(vol, aff)

In [7]: hdr = img.get_header()

In [8]: hdr.set_sform(-aff, 1)

In [9]: hdr.get_best_affine()
Out[9]:
array([[-1., -0., -0., -0.],
      [-0., -1., -0., -0.],
      [-0., -0., -1., -0.],
      [ 0.,  0.,  0.,  1.]])

In [10]: img.update_header()

In [11]: hdr.get_best_affine()
Out[11]:
array([[ 1.,  0.,  0.,  0.],
      [ 0.,  1.,  0.,  0.],
      [ 0.,  0.,  1.,  0.],
      [ 0.,  0.,  0.,  1.]])

I personally think any consistency checks should be done when the affine gets modified and not on write because the following seems really unintuitive to me and I think it'll be confusing to users.

In [18]: hdr.set_sform(-aff, 1)

In [19]: hdr.get_sform()
Out[13]:
array([[-1., -0., -0., -0.],
      [-0., -1., -0., -0.],
      [-0., -0., -1., -0.],
      [ 0.,  0.,  0.,  1.]])

In [20]: nib.save(img, 'my_img.nii')

In [21]: hdr.get_sform()
Out[21]:
array([[ 1.,  0.,  0.,  0.],
      [ 0.,  1.,  0.,  0.],
      [ 0.,  0.,  1.,  0.],
      [ 0.,  0.,  0.,  1.]])

I'm willing to spend some time on this, but it would help me to know how many ways the affine on an 'image' can

Thanks for thinking about this.

The problems are two-fold:

The affine is separate from the header in general. For example an
Analyze image can contain any affine, but can't store one in the
header. Even for nifti, we can't use the same storage for the image
affine and the sform, because the sform is only 3 x 4.

We can't tell in general when the affine is modified. Consider:

img = load('some_image.nii')
img.get_affine()[0, 0] = 6

See you,

Matthew

[MrBago]

Ok, now I understand the scope of the problem better. I'm sure a more complete and elegant solution is possible, but maybe the simplest thing we can do is to create a keyword for the save function which will turn on/off update_header(). It's not that pretty but it'll allow a user to manually manipulate the header and make sure those values get used durring the save. Let me know what you think.

[matthew-brett]

On Wed, Dec 21, 2011 at 1:06 AM, MrBago
[email protected]
wrote:

Ok, now I understand the scope of the problem better. I'm sure a more complete and elegant solution is possible, but maybe the simplest thing we can do is to create a keyword for the save function which will turn on/off update_header(). It's not that pretty but it'll allow a user to manually manipulate the header and make sure those values get used durring the save. Let me know what you think.

I'm a bit reluctant to add keywords to 'save' just on the basis that
it could start getting rather cluttered if we carry on down that
track. Not completely against.

It would not be hard as things stand to prevent the affine updating
the header if the affine is the same as it was at image instantiation.
How far would that go towards solving the problem do you think?

[MrBago]

On 12/22/2011 05:43 AM, Matthew Brett wrote:

I'm a bit reluctant to add keywords to 'save' just on the basis that
it could start getting rather cluttered if we carry on down that
track. Not completely against.
I agree
It would not be hard as things stand to prevent the affine updating
the header if the affine is the same as it was at image instantiation.
How far would that go towards solving the problem do you think?

I feel like that's unnecessarily complex without addressing the root of
the problem. I think the fix you've proposed already takes care of my
use case, and maybe I'm trying to solve a problem that doesn't exist.
Let me give this some more thought of-line (and give others smarter than
me a chance to think about it too). Also if others run into this issue
it'll give us a better idea of the use cases where this matters.

Bago

[moloney]

Shouldn't the Nifti1 update_header just harmonize the affine to the sform and not touch the qform or the codes? Let the users set these header specific bits through the header object.

[moloney]

After thinking about this some more, my above suggestion does not work if the sform_code is 0 (since it would never get updated).

The approach in the master branch (do nothing if the affine matches the header.get_best_affine()) seems like a good solution. I do wonder if the qform should really be getting overwritten when the affines don't match, but at least this can be avoided.

[matthew-brett]

On Tue, Feb 28, 2012 at 8:03 PM, moloney
[email protected]
wrote:

After thinking about this some more, my above suggestion does not work if the sform_code is 0 (since it would never get updated).

The approach in the master branch (do nothing if the affine matches the header.get_best_affine()) seems like a good solution. I do wonder if the qform should really be getting overwritten when the affines don't match, but at least this can be avoided.

Thanks for thinking about this.

The reason it works as it does is that I imagine some situation like this:

some_affine = np.eye(4)
data = np.zeros((2,3,4))
img = Nifti1Image(data, some_affine)

I'm imagining then that you'd want to set both of sform and qform -
because neither are set at the moment. Then imagine I have done some
registration that result in a new affine for the img:

reg_aff = some_registration(img, another_img)

and I want to set this new affine into the image:

img.get_affine()[:] = reg_aff

Now, I could mean - 'only set the sform to this affine leave the qform
as it is'. But I probably don't mean that, I probably mean 'the
affine for this image has changed, please set all the header fields
accordingly when I save the image'.

You can set just the sform (not the qform) with:

img.get_affine()[:] = reg_aff
img.get_header().set_sform(reg_aff)

because then the update_header check will see that the sform and
affine match. I realize this is rather complicated to keep in mind.

[MrBago]

Unfortunately I'm still running into issues with this qform/sform stuff. Here is the issue:

>>> aff = img.get_affine()
>>> aff[:] = np.diag([1.1, 1.1, 1.1, 1])
>>> hdr.set_qform(aff, 2)
>>> hdr.set_sform(aff, 2)
>>> img.update_header()
>>> hdr['qform_code']
array(0, dtype=int16)

I think it's happening because the affine is float64 and the sform is float32:

>>> aff == hdr.get_best_affine()
array([[False,  True,  True,  True],
       [ True, False,  True,  True],
       [ True,  True, False,  True],
       [ True,  True,  True,  True]], dtype=bool)
>>> aff[:] = np.diag([1.1, 1.1, 1.1, 1]).astype('float32')
>>> aff == hdr.get_best_affine()
array([[ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True],
       [ True,  True,  True,  True]], dtype=bool)

I've resorted to the following in order to be able to set the qform to match (as much as possible) the sform:

aff = img.get_affine()
aff[:] = np.diag([1.1, 1.1, 1.1, 1])
img.update_header()
aff[:] = hdr.get_best_affine()
hdr.set_qform(aff, 2)

[moloney]

Seems like the code should be using something like numpy.allclose instead of comparing floats directly.

[matthew-brett]

On Thu, Mar 8, 2012 at 11:47 AM, moloney
[email protected]
wrote:

Seems like the code should be using something like numpy.allclose instead of comparing floats directly.

Sounds reasonable - care to review : #90 ?

Thanks a lot

[matthew-brett]

See discussion at #90 about how to solve sform / qform problem in a more general way.

[matthew-brett]

Closed by the pull request above