1# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
2# vi: set ft=python sts=4 ts=4 sw=4 et:
3### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
4#
5#   See COPYING file distributed along with the NiBabel package for the
6#   copyright and license terms.
7#
8### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
9''' Tests for nifti reading package '''
10from __future__ import with_statement
11import os
12
13from ..py3k import BytesIO, ZEROB, asbytes
14
15import numpy as np
16
17from ..casting import type_info, have_binary128
18from ..tmpdirs import InTemporaryDirectory
19from ..spatialimages import HeaderDataError
20from ..affines import from_matvec
21from .. import nifti1 as nifti1
22from ..nifti1 import (load, Nifti1Header, Nifti1PairHeader, Nifti1Image,
23                      Nifti1Pair, Nifti1Extension, Nifti1Extensions,
24                      data_type_codes, extension_codes, slice_order_codes)
25
26from .test_arraywriters import rt_err_estimate, IUINT_TYPES
27
28from numpy.testing import assert_array_equal, assert_array_almost_equal
29from nose.tools import (assert_true, assert_false, assert_equal,
30                        assert_raises)
31from nose import SkipTest
32
33from ..testing import data_path
34
35from . import test_analyze as tana
36
37header_file = os.path.join(data_path, 'nifti1.hdr')
38image_file = os.path.join(data_path, 'example4d.nii.gz')
39
40
41# Example transformation matrix
42R = [[0, -1, 0], [1, 0, 0], [0, 0, 1]] # rotation matrix
43Z = [2.0, 3.0, 4.0] # zooms
44T = [20, 30, 40] # translations
45A = np.eye(4)
46A[:3,:3] = np.array(R) * Z # broadcasting does the job
47A[:3,3] = T
48
49
50class TestNifti1PairHeader(tana.TestAnalyzeHeader):
51    header_class = Nifti1PairHeader
52    example_file = header_file
53
54    def test_empty(self):
55        tana.TestAnalyzeHeader.test_empty(self)
56        hdr = self.header_class()
57        assert_equal(hdr['magic'], asbytes('ni1'))
58        assert_equal(hdr['scl_slope'], 1)
59        assert_equal(hdr['vox_offset'], 0)
60
61    def test_from_eg_file(self):
62        hdr = Nifti1Header.from_fileobj(open(self.example_file, 'rb'))
63        assert_equal(hdr.endianness, '<')
64        assert_equal(hdr['magic'], asbytes('ni1'))
65        assert_equal(hdr['sizeof_hdr'], 348)
66
67    def test_big_scaling(self):
68        # Test that upcasting works for huge scalefactors
69        # See tests for apply_read_scaling in test_utils
70        hdr = self.header_class()
71        hdr.set_data_shape((2,1,1))
72        hdr.set_data_dtype(np.int16)
73        sio = BytesIO()
74        dtt = np.float32
75        # This will generate a huge scalefactor
76        finf = type_info(dtt)
77        data = np.array([finf['min'], finf['max']], dtype=dtt)[:,None, None]
78        hdr.data_to_fileobj(data, sio)
79        data_back = hdr.data_from_fileobj(sio)
80        assert_true(np.allclose(data, data_back))
81
82    def test_nifti_log_checks(self):
83        # in addition to analyze header checks
84        HC = self.header_class
85        # intercept and slope
86        hdr = HC()
87        # Slope of 0 is OK
88        hdr['scl_slope'] = 0
89        fhdr, message, raiser = self.log_chk(hdr, 0)
90        assert_equal((fhdr, message), (hdr, ''))
91        # But not with non-zero intercept
92        hdr['scl_inter'] = 3
93        fhdr, message, raiser = self.log_chk(hdr, 20)
94        assert_equal(fhdr['scl_inter'], 0)
95        assert_equal(message,
96                           'Unused "scl_inter" is 3.0; should be 0; '
97                           'setting "scl_inter" to 0')
98        # Or not-finite intercept
99        hdr['scl_inter'] = np.nan
100        # NaN string representation can be odd on windows
101        nan_str = '%s' % np.nan
102        fhdr, message, raiser = self.log_chk(hdr, 20)
103        assert_equal(fhdr['scl_inter'], 0)
104        assert_equal(message,
105                           'Unused "scl_inter" is %s; should be 0; '
106                           'setting "scl_inter" to 0' % nan_str)
107        # Reset to usable scale
108        hdr['scl_slope'] = 1
109        # not finite inter is more of a problem
110        hdr['scl_inter'] = np.nan # severity 30
111        fhdr, message, raiser = self.log_chk(hdr, 40)
112        assert_equal(fhdr['scl_inter'], 0)
113        assert_equal(message,
114                           '"scl_slope" is 1.0; but "scl_inter" is %s; '
115                           '"scl_inter" should be finite; setting '
116                           '"scl_inter" to 0' % nan_str)
117        assert_raises(*raiser)
118        # Not finite scale also bad, generates message for scale and offset
119        hdr['scl_slope'] = np.nan
120        fhdr, message, raiser = self.log_chk(hdr, 30)
121        assert_equal(fhdr['scl_slope'], 0)
122        assert_equal(fhdr['scl_inter'], 0)
123        assert_equal(message,
124                           '"scl_slope" is nan; should be finite; '
125                           'Unused "scl_inter" is nan; should be 0; '
126                           'setting "scl_slope" to 0 (no scaling); '
127                           'setting "scl_inter" to 0')
128        assert_raises(*raiser)
129        # Or just scale if inter is already 0
130        hdr['scl_inter'] = 0
131        fhdr, message, raiser = self.log_chk(hdr, 30)
132        assert_equal(fhdr['scl_slope'], 0)
133        assert_equal(fhdr['scl_inter'], 0)
134        assert_equal(message,
135                           '"scl_slope" is nan; should be finite; '
136                           'setting "scl_slope" to 0 (no scaling)')
137        assert_raises(*raiser)
138        # qfac
139        hdr = HC()
140        hdr['pixdim'][0] = 0
141        fhdr, message, raiser = self.log_chk(hdr, 20)
142        assert_equal(fhdr['pixdim'][0], 1)
143        assert_equal(message, 'pixdim[0] (qfac) should be 1 '
144                           '(default) or -1; setting qfac to 1')
145        # magic and offset
146        hdr = HC()
147        hdr['magic'] = 'ooh'
148        fhdr, message, raiser = self.log_chk(hdr, 45)
149        assert_equal(fhdr['magic'], asbytes('ooh'))
150        assert_equal(message, 'magic string "ooh" is not valid; '
151                           'leaving as is, but future errors are likely')
152        hdr['magic'] = 'n+1' # single file needs suitable offset
153        hdr['vox_offset'] = 0
154        fhdr, message, raiser = self.log_chk(hdr, 40)
155        assert_equal(fhdr['vox_offset'], 352)
156        assert_equal(message, 'vox offset 0 too low for single '
157                           'file nifti1; setting to minimum value '
158                           'of 352')
159        # qform, sform
160        hdr = HC()
161        hdr['qform_code'] = -1
162        fhdr, message, raiser = self.log_chk(hdr, 30)
163        assert_equal(fhdr['qform_code'], 0)
164        assert_equal(message, 'qform_code -1 not valid; '
165                           'setting to 0')
166        hdr = HC()
167        hdr['sform_code'] = -1
168        fhdr, message, raiser = self.log_chk(hdr, 30)
169        assert_equal(fhdr['sform_code'], 0)
170        assert_equal(message, 'sform_code -1 not valid; '
171                           'setting to 0')
172
173    def test_freesurfer_hack(self):
174        # For large vector images, Freesurfer appears to set dim[1] to -1 and
175        # then use glmin for the vector length (an i4)
176        HC = self.header_class
177        # The standard case
178        hdr = HC()
179        hdr.set_data_shape((2, 3, 4))
180        assert_equal(hdr.get_data_shape(), (2, 3, 4))
181        assert_equal(hdr['glmin'], 0)
182        # Just left of the freesurfer case
183        dim_type = hdr.template_dtype['dim'].base
184        glmin = hdr.template_dtype['glmin'].base
185        too_big = int(np.iinfo(dim_type).max) + 1
186        hdr.set_data_shape((too_big-1, 1, 1))
187        assert_equal(hdr.get_data_shape(), (too_big-1, 1, 1))
188        # The freesurfer case
189        hdr.set_data_shape((too_big, 1, 1))
190        assert_equal(hdr.get_data_shape(), (too_big, 1, 1))
191        assert_array_equal(hdr['dim'][:4], [3, -1, 1, 1])
192        assert_equal(hdr['glmin'], too_big)
193        # This only works for the case of a 3D with -1, 1, 1
194        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,))
195        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,1))
196        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,1,2))
197        assert_raises(HeaderDataError, hdr.set_data_shape, (too_big,2,1))
198        assert_raises(HeaderDataError, hdr.set_data_shape, (1, too_big))
199        assert_raises(HeaderDataError, hdr.set_data_shape, (1, too_big, 1))
200        assert_raises(HeaderDataError, hdr.set_data_shape, (1, 1, too_big))
201        # Outside range of glmin raises error
202        far_too_big = int(np.iinfo(glmin).max) + 1
203        hdr.set_data_shape((far_too_big-1, 1, 1))
204        assert_equal(hdr.get_data_shape(), (far_too_big-1, 1, 1))
205        assert_raises(HeaderDataError, hdr.set_data_shape, (far_too_big,1,1))
206        # glmin of zero raises error (implausible vector length)
207        hdr.set_data_shape((-1,1,1))
208        hdr['glmin'] = 0
209        assert_raises(HeaderDataError, hdr.get_data_shape)
210        # Lists or tuples or arrays will work for setting shape
211        for shape in ((too_big-1, 1, 1), (too_big, 1, 1)):
212            for constructor in (list, tuple, np.array):
213                hdr.set_data_shape(constructor(shape))
214                assert_equal(hdr.get_data_shape(), shape)
215
216
217    def test_qform_sform(self):
218        HC = self.header_class
219        hdr = HC()
220        assert_array_equal(hdr.get_qform(), np.eye(4))
221        empty_sform = np.zeros((4,4))
222        empty_sform[-1,-1] = 1
223        assert_array_equal(hdr.get_sform(), empty_sform)
224        assert_equal(hdr.get_qform(coded=True), (None, 0))
225        assert_equal(hdr.get_sform(coded=True), (None, 0))
226        # Affine with no shears
227        nice_aff = np.diag([2, 3, 4, 1])
228        # Affine with shears
229        nasty_aff = from_matvec(np.arange(9).reshape((3,3)), [9, 10, 11])
230        fixed_aff = unshear_44(nasty_aff)
231        for in_meth, out_meth in ((hdr.set_qform, hdr.get_qform),
232                                  (hdr.set_sform, hdr.get_sform)):
233            in_meth(nice_aff, 2)
234            aff, code = out_meth(coded=True)
235            assert_array_equal(aff, nice_aff)
236            assert_equal(code, 2)
237            assert_array_equal(out_meth(), nice_aff) # non coded
238            # Affine can also be passed if code == 0, affine will be suitably set
239            in_meth(nice_aff, 0)
240            assert_equal(out_meth(coded=True), (None, 0))
241            assert_array_almost_equal(out_meth(), nice_aff)
242            # Default qform code when previous == 0 is 2
243            in_meth(nice_aff)
244            aff, code = out_meth(coded=True)
245            assert_equal(code, 2)
246            # Unless code was non-zero before
247            in_meth(nice_aff, 1)
248            in_meth(nice_aff)
249            aff, code = out_meth(coded=True)
250            assert_equal(code, 1)
251            # Can set code without modifying affine, by passing affine=None
252            assert_array_equal(aff, nice_aff) # affine same as before
253            in_meth(None, 3)
254            aff, code = out_meth(coded=True)
255            assert_array_equal(aff, nice_aff) # affine same as before
256            assert_equal(code, 3)
257            # affine is None on its own, or with code==0, resets code to 0
258            in_meth(None, 0)
259            assert_equal(out_meth(coded=True), (None, 0))
260            in_meth(None)
261            assert_equal(out_meth(coded=True), (None, 0))
262            # List works as input
263            in_meth(nice_aff.tolist())
264            assert_array_equal(out_meth(), nice_aff)
265        # Qform specifics
266        # inexact set (with shears) is OK
267        hdr.set_qform(nasty_aff, 1)
268        assert_array_almost_equal(hdr.get_qform(), fixed_aff)
269        # Unless allow_shears is False
270        assert_raises(HeaderDataError, hdr.set_qform, nasty_aff, 1, False)
271        # Reset sform, give qform a code, to test sform
272        hdr.set_sform(None)
273        hdr.set_qform(nice_aff, 1)
274        # Check sform unchanged by setting qform
275        assert_equal(hdr.get_sform(coded=True), (None, 0))
276        # Setting does change the sform ouput
277        hdr.set_sform(nasty_aff, 1)
278        aff, code = hdr.get_sform(coded=True)
279        assert_array_equal(aff, nasty_aff)
280        assert_equal(code, 1)
281
282
283def unshear_44(affine):
284    RZS = affine[:3, :3]
285    zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
286    R = RZS / zooms
287    P, S, Qs = np.linalg.svd(R)
288    PR = np.dot(P, Qs)
289    return from_matvec(PR * zooms, affine[:3,3])
290
291
292class TestNifti1SingleHeader(TestNifti1PairHeader):
293
294    header_class = Nifti1Header
295
296    def test_empty(self):
297        tana.TestAnalyzeHeader.test_empty(self)
298        hdr = self.header_class()
299        assert_equal(hdr['magic'], asbytes('n+1'))
300        assert_equal(hdr['scl_slope'], 1)
301        assert_equal(hdr['vox_offset'], 352)
302
303    def test_binblock_is_file(self):
304        # Override test that binary string is the same as the file on disk; in
305        # the case of the single file version of the header, we need to append
306        # the extension string (4 0s)
307        hdr = self.header_class()
308        str_io = BytesIO()
309        hdr.write_to(str_io)
310        assert_equal(str_io.getvalue(), hdr.binaryblock + ZEROB * 4)
311
312    def test_float128(self):
313        hdr = self.header_class()
314        if have_binary128():
315            hdr.set_data_dtype(np.longdouble)
316            assert_equal(hdr.get_data_dtype().type, np.longdouble)
317        else:
318            assert_raises(HeaderDataError, hdr.set_data_dtype, np.longdouble)
319
320
321class TestNifti1Image(tana.TestAnalyzeImage):
322    # Run analyze-flavor spatialimage tests
323    image_class = Nifti1Image
324
325    def _qform_rt(self, img):
326        # Round trip image after setting qform, sform codes
327        hdr = img.get_header()
328        hdr['qform_code'] = 3
329        hdr['sform_code'] = 4
330        # Save / reload using bytes IO objects
331        for key, value in img.file_map.items():
332            value.fileobj = BytesIO()
333        img.to_file_map()
334        return img.from_file_map(img.file_map)
335
336    def test_qform_cycle(self):
337        # Qform load save cycle
338        img_klass = self.image_class
339        # None affine
340        img = img_klass(np.zeros((2,3,4)), None)
341        hdr_back = self._qform_rt(img).get_header()
342        assert_equal(hdr_back['qform_code'], 3)
343        assert_equal(hdr_back['sform_code'], 4)
344        # Try non-None affine
345        img = img_klass(np.zeros((2,3,4)), np.eye(4))
346        hdr_back = self._qform_rt(img).get_header()
347        assert_equal(hdr_back['qform_code'], 3)
348        assert_equal(hdr_back['sform_code'], 4)
349        # Modify affine in-place - does it hold?
350        img.get_affine()[0,0] = 9
351        img.to_file_map()
352        img_back = img.from_file_map(img.file_map)
353        exp_aff = np.diag([9,1,1,1])
354        assert_array_equal(img_back.get_affine(), exp_aff)
355        hdr_back = img.get_header()
356        assert_array_equal(hdr_back.get_sform(), exp_aff)
357        assert_array_equal(hdr_back.get_qform(), exp_aff)
358
359    def test_header_update_affine(self):
360        # Test that updating occurs only if affine is not allclose
361        img = self.image_class(np.zeros((2,3,4)), np.eye(4))
362        hdr = img.get_header()
363        aff = img.get_affine()
364        aff[:] = np.diag([1.1, 1.1, 1.1, 1]) # inexact floats
365        hdr.set_qform(aff, 2)
366        hdr.set_sform(aff, 2)
367        img.update_header()
368        assert_equal(hdr['sform_code'], 2)
369        assert_equal(hdr['qform_code'], 2)
370
371    def test_set_qform(self):
372        img = self.image_class(np.zeros((2,3,4)), np.diag([2.2, 3.3, 4.3, 1]))
373        hdr = img.get_header()
374        new_affine = np.diag([1.1, 1.1, 1.1, 1])
375        # Affine is same as sform (best affine)
376        assert_array_almost_equal(img.get_affine(), hdr.get_best_affine())
377        # Reset affine to something different again
378        aff_affine = np.diag([3.3, 4.5, 6.6, 1])
379        img.get_affine()[:] = aff_affine
380        assert_array_almost_equal(img.get_affine(), aff_affine)
381        # Set qform using new_affine
382        img.set_qform(new_affine, 1)
383        assert_array_almost_equal(img.get_qform(), new_affine)
384        assert_equal(hdr['qform_code'], 1)
385        # Image get is same as header get
386        assert_array_almost_equal(img.get_qform(), new_affine)
387        # Coded version of get gets same information
388        qaff, code = img.get_qform(coded=True)
389        assert_equal(code, 1)
390        assert_array_almost_equal(qaff, new_affine)
391        # Image affine now reset to best affine (which is sform)
392        assert_array_almost_equal(img.get_affine(), hdr.get_best_affine())
393        # Reset image affine and try update_affine == False
394        img.get_affine()[:] = aff_affine
395        img.set_qform(new_affine, 1, update_affine=False)
396        assert_array_almost_equal(img.get_affine(), aff_affine)
397        # Clear qform using None, zooms unchanged
398        assert_array_almost_equal(hdr.get_zooms(), [1.1, 1.1, 1.1])
399        img.set_qform(None)
400        qaff, code = img.get_qform(coded=True)
401        assert_equal((qaff, code), (None, 0))
402        assert_array_almost_equal(hdr.get_zooms(), [1.1, 1.1, 1.1])
403        # Best affine similarly
404        assert_array_almost_equal(img.get_affine(), hdr.get_best_affine())
405        # If sform is not set, qform should update affine
406        img.set_sform(None)
407        img.set_qform(new_affine, 1)
408        qaff, code = img.get_qform(coded=True)
409        assert_equal(code, 1)
410        assert_array_almost_equal(img.get_affine(), new_affine)
411        new_affine[0, 1] = 2
412        # If affine has has shear, should raise Error if strip_shears=False
413        img.set_qform(new_affine, 2)
414        assert_raises(HeaderDataError, img.set_qform, new_affine, 2, False)
415        # Unexpected keyword raises error
416        assert_raises(TypeError, img.get_qform, strange=True)
417
418    def test_set_sform(self):
419        orig_aff = np.diag([2.2, 3.3, 4.3, 1])
420        img = self.image_class(np.zeros((2,3,4)), orig_aff)
421        hdr = img.get_header()
422        new_affine = np.diag([1.1, 1.1, 1.1, 1])
423        qform_affine = np.diag([1.2, 1.2, 1.2, 1])
424        # Reset image affine to something different again
425        aff_affine = np.diag([3.3, 4.5, 6.6, 1])
426        img.get_affine()[:] = aff_affine
427        assert_array_almost_equal(img.get_affine(), aff_affine)
428        # Sform, Qform codes are 'aligned',  'unknown' by default
429        assert_equal((hdr['sform_code'], hdr['qform_code']), (2, 0))
430        # Set sform using new_affine when qform is 0
431        img.set_sform(new_affine, 1)
432        assert_equal(hdr['sform_code'], 1)
433        assert_array_almost_equal(hdr.get_sform(), new_affine)
434        # Image get is same as header get
435        assert_array_almost_equal(img.get_sform(), new_affine)
436        # Coded version gives same result
437        saff, code = img.get_sform(coded=True)
438        assert_equal(code, 1)
439        assert_array_almost_equal(saff, new_affine)
440        # Because we've reset the sform with update_affine, the affine changes
441        assert_array_almost_equal(img.get_affine(), hdr.get_best_affine())
442        # Reset image affine and try update_affine == False
443        img.get_affine()[:] = aff_affine
444        img.set_sform(new_affine, 1, update_affine=False)
445        assert_array_almost_equal(img.get_affine(), aff_affine)
446        # zooms do not get updated when qform is 0
447        assert_array_almost_equal(img.get_qform(), orig_aff)
448        assert_array_almost_equal(hdr.get_zooms(), [2.2, 3.3, 4.3])
449        img.set_qform(None)
450        assert_array_almost_equal(hdr.get_zooms(), [2.2, 3.3, 4.3])
451        # Set sform using new_affine when qform is set
452        img.set_qform(qform_affine, 1)
453        img.set_sform(new_affine, 1)
454        saff, code = img.get_sform(coded=True)
455        assert_equal(code, 1)
456        assert_array_almost_equal(saff, new_affine)
457        assert_array_almost_equal(img.get_affine(), new_affine)
458        # zooms follow qform
459        assert_array_almost_equal(hdr.get_zooms(), [1.2, 1.2, 1.2])
460        # Clear sform using None, best_affine should fall back on qform
461        img.set_sform(None)
462        assert_equal(hdr['sform_code'], 0)
463        assert_equal(hdr['qform_code'], 1)
464        # Sform holds previous affine from last set
465        assert_array_almost_equal(hdr.get_sform(), saff)
466        # Image affine follows qform
467        assert_array_almost_equal(img.get_affine(), qform_affine)
468        assert_array_almost_equal(hdr.get_best_affine(), img.get_affine())
469        # Unexpected keyword raises error
470        assert_raises(TypeError, img.get_sform, strange=True)
471
472    def test_hdr_diff(self):
473        # Check an offset beyond data does not raise an error
474        img = self.image_class(np.zeros((2,3,4)), np.eye(4))
475        ext = dict(img.files_types)['image']
476        hdr = img.get_header()
477        hdr['vox_offset'] = 400
478        with InTemporaryDirectory():
479            img.to_filename('another_file' + ext)
480
481
482class TestNifti1Pair(TestNifti1Image):
483    # Run analyze-flavor spatialimage tests
484    image_class = Nifti1Pair
485
486
487def test_datatypes():
488    hdr = Nifti1Header()
489    for code in data_type_codes.value_set():
490        dt = data_type_codes.type[code]
491        if dt == np.void:
492            continue
493        hdr.set_data_dtype(code)
494        (assert_equal,
495               hdr.get_data_dtype(),
496               data_type_codes.dtype[code])
497    # Check that checks also see new datatypes
498    hdr.set_data_dtype(np.complex128)
499    hdr.check_fix()
500
501
502def test_quaternion():
503    hdr = Nifti1Header()
504    hdr['quatern_b'] = 0
505    hdr['quatern_c'] = 0
506    hdr['quatern_d'] = 0
507    assert_true(np.allclose(hdr.get_qform_quaternion(), [1.0, 0, 0, 0]))
508    hdr['quatern_b'] = 1
509    hdr['quatern_c'] = 0
510    hdr['quatern_d'] = 0
511    assert_true(np.allclose(hdr.get_qform_quaternion(), [0, 1, 0, 0]))
512    # Check threshold set correctly for float32
513    hdr['quatern_b'] = 1+np.finfo(np.float32).eps
514    assert_array_almost_equal(hdr.get_qform_quaternion(), [0, 1, 0, 0])
515
516
517def test_qform():
518    # Test roundtrip case
519    ehdr = Nifti1Header()
520    ehdr.set_qform(A)
521    qA = ehdr.get_qform()
522    assert_true, np.allclose(A, qA, atol=1e-5)
523    assert_true, np.allclose(Z, ehdr['pixdim'][1:4])
524    xfas = nifti1.xform_codes
525    assert_true, ehdr['qform_code'] == xfas['aligned']
526    ehdr.set_qform(A, 'scanner')
527    assert_true, ehdr['qform_code'] == xfas['scanner']
528    ehdr.set_qform(A, xfas['aligned'])
529    assert_true, ehdr['qform_code'] == xfas['aligned']
530
531
532def test_sform():
533    # Test roundtrip case
534    ehdr = Nifti1Header()
535    ehdr.set_sform(A)
536    sA = ehdr.get_sform()
537    assert_true, np.allclose(A, sA, atol=1e-5)
538    xfas = nifti1.xform_codes
539    assert_true, ehdr['sform_code'] == xfas['aligned']
540    ehdr.set_sform(A, 'scanner')
541    assert_true, ehdr['sform_code'] == xfas['scanner']
542    ehdr.set_sform(A, xfas['aligned'])
543    assert_true, ehdr['sform_code'] == xfas['aligned']
544
545
546def test_dim_info():
547    ehdr = Nifti1Header()
548    assert_true(ehdr.get_dim_info() == (None, None, None))
549    for info in ((0,2,1),
550                 (None, None, None),
551                 (0,2,None),
552                 (0,None,None),
553                 (None,2,1),
554                 (None, None,1),
555                 ):
556        ehdr.set_dim_info(*info)
557        assert_true(ehdr.get_dim_info() == info)
558
559
560def test_slice_times():
561    hdr = Nifti1Header()
562    # error if slice dimension not specified
563    assert_raises(HeaderDataError, hdr.get_slice_times)
564    hdr.set_dim_info(slice=2)
565    # error if slice dimension outside shape
566    assert_raises(HeaderDataError, hdr.get_slice_times)
567    hdr.set_data_shape((1, 1, 7))
568    # error if slice duration not set
569    assert_raises(HeaderDataError, hdr.get_slice_times)
570    hdr.set_slice_duration(0.1)
571    # We need a function to print out the Nones and floating point
572    # values in a predictable way, for the tests below.
573    _stringer = lambda val: val is not None and '%2.1f' % val or None
574    _print_me = lambda s: map(_stringer, s)
575    #The following examples are from the nifti1.h documentation.
576    hdr['slice_code'] = slice_order_codes['sequential increasing']
577    assert_equal(_print_me(hdr.get_slice_times()),
578                       ['0.0', '0.1', '0.2', '0.3', '0.4',
579                        '0.5', '0.6'])
580    hdr['slice_start'] = 1
581    hdr['slice_end'] = 5
582    assert_equal(_print_me(hdr.get_slice_times()),
583        [None, '0.0', '0.1', '0.2', '0.3', '0.4', None])
584    hdr['slice_code'] = slice_order_codes['sequential decreasing']
585    assert_equal(_print_me(hdr.get_slice_times()),
586        [None, '0.4', '0.3', '0.2', '0.1', '0.0', None])
587    hdr['slice_code'] = slice_order_codes['alternating increasing']
588    assert_equal(_print_me(hdr.get_slice_times()),
589        [None, '0.0', '0.3', '0.1', '0.4', '0.2', None])
590    hdr['slice_code'] = slice_order_codes['alternating decreasing']
591    assert_equal(_print_me(hdr.get_slice_times()),
592        [None, '0.2', '0.4', '0.1', '0.3', '0.0', None])
593    hdr['slice_code'] = slice_order_codes['alternating increasing 2']
594    assert_equal(_print_me(hdr.get_slice_times()),
595        [None, '0.2', '0.0', '0.3', '0.1', '0.4', None])
596    hdr['slice_code'] = slice_order_codes['alternating decreasing 2']
597    assert_equal(_print_me(hdr.get_slice_times()),
598        [None, '0.4', '0.1', '0.3', '0.0', '0.2', None])
599    # test set
600    hdr = Nifti1Header()
601    hdr.set_dim_info(slice=2)
602    # need slice dim to correspond with shape
603    times = [None, 0.2, 0.4, 0.1, 0.3, 0.0, None]
604    assert_raises(HeaderDataError, hdr.set_slice_times, times)
605    hdr.set_data_shape([1, 1, 7])
606    assert_raises(HeaderDataError,
607                        hdr.set_slice_times,
608                        times[:-1]) # wrong length
609    assert_raises(HeaderDataError,
610                        hdr.set_slice_times,
611                        (None,) * len(times)) # all None
612    n_mid_times = times[:]
613    n_mid_times[3] = None
614    assert_raises(HeaderDataError,
615                        hdr.set_slice_times,
616                        n_mid_times) # None in middle
617    funny_times = times[:]
618    funny_times[3] = 0.05
619    assert_raises(HeaderDataError,
620                        hdr.set_slice_times,
621                        funny_times) # can't get single slice duration
622    hdr.set_slice_times(times)
623    assert_equal(hdr.get_value_label('slice_code'),
624                       'alternating decreasing')
625    assert_equal(hdr['slice_start'], 1)
626    assert_equal(hdr['slice_end'], 5)
627    assert_array_almost_equal(hdr['slice_duration'], 0.1)
628
629
630def test_intents():
631    ehdr = Nifti1Header()
632    ehdr.set_intent('t test', (10,), name='some score')
633    assert_equal(ehdr.get_intent(),
634                 ('t test', (10.0,), 'some score'))
635    # invalid intent name
636    assert_raises(KeyError,
637                  ehdr.set_intent, 'no intention')
638    # too many parameters
639    assert_raises(HeaderDataError,
640                  ehdr.set_intent,
641                  't test', (10,10))
642    # too few parameters
643    assert_raises(HeaderDataError,
644                  ehdr.set_intent,
645                  'f test', (10,))
646    # check unset parameters are set to 0, and name to ''
647    ehdr.set_intent('t test')
648    assert_equal((ehdr['intent_p1'],
649                  ehdr['intent_p2'],
650                  ehdr['intent_p3']), (0,0,0))
651    assert_equal(ehdr['intent_name'], asbytes(''))
652    ehdr.set_intent('t test', (10,))
653    assert_equal((ehdr['intent_p2'],
654                  ehdr['intent_p3']), (0,0))
655
656
657def test_set_slice_times():
658    hdr = Nifti1Header()
659    hdr.set_dim_info(slice=2)
660    hdr.set_data_shape([1, 1, 7])
661    hdr.set_slice_duration(0.1)
662    times = [0] * 6
663    assert_raises(HeaderDataError, hdr.set_slice_times, times)
664    times = [None] * 7
665    assert_raises(HeaderDataError, hdr.set_slice_times, times)
666    times = [None, 0, 1, None, 3, 4, None]
667    assert_raises(HeaderDataError, hdr.set_slice_times, times)
668    times = [None, 0, 1, 2.1, 3, 4, None]
669    assert_raises(HeaderDataError, hdr.set_slice_times, times)
670    times = [None, 0, 4, 3, 2, 1, None]
671    assert_raises(HeaderDataError, hdr.set_slice_times, times)
672    times = [0, 1, 2, 3, 4, 5, 6]
673    hdr.set_slice_times(times)
674    assert_equal(hdr['slice_code'], 1)
675    assert_equal(hdr['slice_start'], 0)
676    assert_equal(hdr['slice_end'], 6)
677    assert_equal(hdr['slice_duration'], 1.0)
678    times = [None, 0, 1, 2, 3, 4, None]
679    hdr.set_slice_times(times)
680    assert_equal(hdr['slice_code'], 1)
681    assert_equal(hdr['slice_start'], 1)
682    assert_equal(hdr['slice_end'], 5)
683    assert_equal(hdr['slice_duration'], 1.0)
684    times = [None, 0.4, 0.3, 0.2, 0.1, 0, None]
685    hdr.set_slice_times(times)
686    assert_true(np.allclose(hdr['slice_duration'], 0.1))
687    times = [None, 4, 3, 2, 1, 0, None]
688    hdr.set_slice_times(times)
689    assert_equal(hdr['slice_code'], 2)
690    times = [None, 0, 3, 1, 4, 2, None]
691    hdr.set_slice_times(times)
692    assert_equal(hdr['slice_code'], 3)
693    times = [None, 2, 4, 1, 3, 0, None]
694    hdr.set_slice_times(times)
695    assert_equal(hdr['slice_code'], 4)
696    times = [None, 2, 0, 3, 1, 4, None]
697    hdr.set_slice_times(times)
698    assert_equal(hdr['slice_code'], 5)
699    times = [None, 4, 1, 3, 0, 2, None]
700    hdr.set_slice_times(times)
701    assert_equal(hdr['slice_code'], 6)
702
703
704def test_nifti1_images():
705    shape = (2, 4, 6)
706    npt = np.float32
707    data = np.arange(np.prod(shape), dtype=npt).reshape(shape)
708    affine = np.diag([1, 2, 3, 1])
709    img = Nifti1Image(data, affine)
710    assert_equal(img.shape, shape)
711    img.set_data_dtype(npt)
712    stio = BytesIO()
713    img.file_map['image'].fileobj = stio
714    img.to_file_map()
715    img2 = Nifti1Image.from_file_map(img.file_map)
716    assert_array_equal(img2.get_data(), data)
717    with InTemporaryDirectory() as tmpdir:
718        for ext in ('.gz', '.bz2'):
719            fname = os.path.join(tmpdir, 'test.nii' + ext)
720            img.to_filename(fname)
721            img3 = Nifti1Image.load(fname)
722            assert_true(isinstance(img3, img.__class__))
723            assert_array_equal(img3.get_data(), data)
724            assert_equal(img3.get_header(), img.get_header())
725            # del to avoid windows errors of form 'The process cannot
726            # access the file because it is being used'
727            del img3
728
729
730def test_extension_basics():
731    raw = '123'
732    ext = Nifti1Extension('comment', raw)
733    assert_true(ext.get_sizeondisk() == 16)
734    assert_true(ext.get_content() == raw)
735    assert_true(ext.get_code() == 6)
736
737
738def test_ext_eq():
739    ext = Nifti1Extension('comment', '123')
740    assert_true(ext == ext)
741    assert_false(ext != ext)
742    ext2 = Nifti1Extension('comment', '124')
743    assert_false(ext == ext2)
744    assert_true(ext != ext2)
745
746
747def test_extension_codes():
748    for k in extension_codes.keys():
749        ext = Nifti1Extension(k, 'somevalue')
750
751
752def test_extension_list():
753    ext_c0 = Nifti1Extensions()
754    ext_c1 = Nifti1Extensions()
755    assert_equal(ext_c0, ext_c1)
756    ext = Nifti1Extension('comment', '123')
757    ext_c1.append(ext)
758    assert_false(ext_c0 == ext_c1)
759    ext_c0.append(ext)
760    assert_true(ext_c0 == ext_c1)
761
762
763def test_nifti_extensions():
764    nim = load(image_file)
765    # basic checks of the available extensions
766    hdr = nim.get_header()
767    exts_container = hdr.extensions
768    assert_equal(len(exts_container), 2)
769    assert_equal(exts_container.count('comment'), 2)
770    assert_equal(exts_container.count('afni'), 0)
771    assert_equal(exts_container.get_codes(), [6, 6])
772    assert_equal((exts_container.get_sizeondisk()) % 16, 0)
773    # first extension should be short one
774    assert_equal(exts_container[0].get_content(), asbytes('extcomment1'))
775    # add one
776    afniext = Nifti1Extension('afni', '<xml></xml>')
777    exts_container.append(afniext)
778    assert_true(exts_container.get_codes() == [6, 6, 4])
779    assert_true(exts_container.count('comment') == 2)
780    assert_true(exts_container.count('afni') == 1)
781    assert_true((exts_container.get_sizeondisk()) % 16 == 0)
782    # delete one
783    del exts_container[1]
784    assert_true(exts_container.get_codes() == [6, 4])
785    assert_true(exts_container.count('comment') == 1)
786    assert_true(exts_container.count('afni') == 1)
787
788
789def test_loadsave_cycle():
790    nim = load(image_file)
791    # ensure we have extensions
792    hdr = nim.get_header()
793    exts_container = hdr.extensions
794    assert_true(len(exts_container) > 0)
795    # write into the air ;-)
796    stio = BytesIO()
797    nim.file_map['image'].fileobj = stio
798    nim.to_file_map()
799    stio.seek(0)
800    # reload
801    lnim = Nifti1Image.from_file_map(nim.file_map)
802    hdr = lnim.get_header()
803    lexts_container = hdr.extensions
804    assert_equal(exts_container,
805                 lexts_container)
806    # build int16 image
807    data = np.ones((2,3,4,5), dtype='int16')
808    img = Nifti1Image(data, np.eye(4))
809    hdr = img.get_header()
810    assert_equal(hdr.get_data_dtype(), np.int16)
811    # default should have no scaling
812    assert_equal(hdr.get_slope_inter(), (1.0, 0.0))
813    # set scaling
814    hdr.set_slope_inter(2, 8)
815    assert_equal(hdr.get_slope_inter(), (2, 8))
816    # now build new image with updated header
817    wnim = Nifti1Image(data, np.eye(4), header=hdr)
818    assert_equal(wnim.get_data_dtype(), np.int16)
819    assert_equal(wnim.get_header().get_slope_inter(), (2, 8))
820    # write into the air again ;-)
821    stio = BytesIO()
822    wnim.file_map['image'].fileobj = stio
823    wnim.to_file_map()
824    stio.seek(0)
825    lnim = Nifti1Image.from_file_map(wnim.file_map)
826    assert_equal(lnim.get_data_dtype(), np.int16)
827    # the test below does not pass, because the slope and inter are
828    # always reset from the data, by the image write
829    raise SkipTest
830    assert_equal(lnim.get_header().get_slope_inter(), (2, 8))
831
832
833def test_slope_inter():
834    hdr = Nifti1Header()
835    assert_equal(hdr.get_slope_inter(), (1.0, 0.0))
836    for intup, outup in (((2.0,), (2.0, 0.0)),
837                         ((None,), (None, None)),
838                         ((3.0, None), (3.0, 0.0)),
839                         ((0.0, None), (None, None)),
840                         ((None, 0.0), (None, None)),
841                         ((None, 3.0), (None, None)),
842                         ((2.0, 3.0), (2.0, 3.0))):
843        hdr.set_slope_inter(*intup)
844        assert_equal(hdr.get_slope_inter(), outup)
845        # Check set survives through checking
846        hdr = Nifti1Header.from_header(hdr, check=True)
847        assert_equal(hdr.get_slope_inter(), outup)
848
849
850def test_xyzt_units():
851    hdr = Nifti1Header()
852    assert_equal(hdr.get_xyzt_units(), ('unknown', 'unknown'))
853    hdr.set_xyzt_units('mm', 'sec')
854    assert_equal(hdr.get_xyzt_units(), ('mm', 'sec'))
855    hdr.set_xyzt_units()
856    assert_equal(hdr.get_xyzt_units(), ('unknown', 'unknown'))
857
858
859def test_recoded_fields():
860    hdr = Nifti1Header()
861    assert_equal(hdr.get_value_label('qform_code'), 'unknown')
862    hdr['qform_code'] = 3
863    assert_equal(hdr.get_value_label('qform_code'), 'talairach')
864    assert_equal(hdr.get_value_label('sform_code'), 'unknown')
865    hdr['sform_code'] = 3
866    assert_equal(hdr.get_value_label('sform_code'), 'talairach')
867    assert_equal(hdr.get_value_label('intent_code'), 'none')
868    hdr.set_intent('t test', (10,), name='some score')
869    assert_equal(hdr.get_value_label('intent_code'), 't test')
870    assert_equal(hdr.get_value_label('slice_code'), 'unknown')
871    hdr['slice_code'] = 4 # alternating decreasing
872    assert_equal(hdr.get_value_label('slice_code'),
873                       'alternating decreasing')
874
875
876def test_load():
877    # test module level load.  We try to load a nii and an .img and a .hdr and
878    # expect to get a nifti back of single or pair type
879    arr = np.arange(24).reshape((2,3,4))
880    aff = np.diag([2, 3, 4, 1])
881    simg = Nifti1Image(arr, aff)
882    pimg = Nifti1Pair(arr, aff)
883    with InTemporaryDirectory():
884        nifti1.save(simg, 'test.nii')
885        assert_array_equal(arr, nifti1.load('test.nii').get_data())
886        nifti1.save(simg, 'test.img')
887        assert_array_equal(arr, nifti1.load('test.img').get_data())
888        nifti1.save(simg, 'test.hdr')
889        assert_array_equal(arr, nifti1.load('test.hdr').get_data())
890
891
892def test_load_pixdims():
893    # Make sure load preserves separate qform, pixdims, sform
894    arr = np.arange(24).reshape((2,3,4))
895    qaff = np.diag([2, 3, 4, 1])
896    saff = np.diag([5, 6, 7, 1])
897    hdr = Nifti1Header()
898    hdr.set_qform(qaff)
899    assert_array_equal(hdr.get_qform(), qaff)
900    hdr.set_sform(saff)
901    assert_array_equal(hdr.get_sform(), saff)
902    simg = Nifti1Image(arr, None, hdr)
903    img_hdr = simg.get_header()
904    # Check qform, sform, pixdims are the same
905    assert_array_equal(img_hdr.get_qform(), qaff)
906    assert_array_equal(img_hdr.get_sform(), saff)
907    assert_array_equal(img_hdr.get_zooms(), [2,3,4])
908    # Save to stringio
909    fm = Nifti1Image.make_file_map()
910    fm['image'].fileobj = BytesIO()
911    simg.to_file_map(fm)
912    # Load again
913    re_simg = Nifti1Image.from_file_map(fm)
914    assert_array_equal(re_simg.get_data(), arr)
915    # Check qform, sform, pixdims are the same
916    rimg_hdr = re_simg.get_header()
917    assert_array_equal(rimg_hdr.get_qform(), qaff)
918    assert_array_equal(rimg_hdr.get_sform(), saff)
919    assert_array_equal(rimg_hdr.get_zooms(), [2,3,4])
920
921
922def test_affines_init():
923    # Test we are doing vaguely spec-related qform things.  The 'spec' here is
924    # some thoughts by Mark Jenkinson:
925    # http://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/qsform_brief_usage
926    arr = np.arange(24).reshape((2,3,4))
927    aff = np.diag([2, 3, 4, 1])
928    # Default is sform set, qform not set
929    img = Nifti1Image(arr, aff)
930    hdr = img.get_header()
931    assert_equal(hdr['qform_code'], 0)
932    assert_equal(hdr['sform_code'], 2)
933    assert_array_equal(hdr.get_zooms(), [2, 3, 4])
934    # This is also true for affines with header passed
935    qaff = np.diag([3, 4, 5, 1])
936    saff = np.diag([6, 7, 8, 1])
937    hdr.set_qform(qaff, code='scanner')
938    hdr.set_sform(saff, code='talairach')
939    assert_array_equal(hdr.get_zooms(), [3, 4, 5])
940    img = Nifti1Image(arr, aff, hdr)
941    new_hdr = img.get_header()
942    # Again affine is sort of anonymous space
943    assert_equal(new_hdr['qform_code'], 0)
944    assert_equal(new_hdr['sform_code'], 2)
945    assert_array_equal(new_hdr.get_sform(), aff)
946    assert_array_equal(new_hdr.get_zooms(), [2, 3, 4])
947    # But if no affine passed, codes and matrices stay the same
948    img = Nifti1Image(arr, None, hdr)
949    new_hdr = img.get_header()
950    assert_equal(new_hdr['qform_code'], 1) # scanner
951    assert_array_equal(new_hdr.get_qform(), qaff)
952    assert_equal(new_hdr['sform_code'], 3) # Still talairach
953    assert_array_equal(new_hdr.get_sform(), saff)
954    # Pixdims as in the original header
955    assert_array_equal(new_hdr.get_zooms(), [3, 4, 5])
956
957
958def round_trip(img):
959    stio = BytesIO()
960    img.file_map['image'].fileobj = stio
961    img.to_file_map()
962    return Nifti1Image.from_file_map(img.file_map)
963
964
965def test_float_int_min_max():
966    # Conversion between float and int
967    # Parallel test to arraywriters
968    aff = np.eye(4)
969    for in_dt in (np.float32, np.float64):
970        finf = type_info(in_dt)
971        arr = np.array([finf['min'], finf['max']], dtype=in_dt)
972        for out_dt in IUINT_TYPES:
973            img = Nifti1Image(arr, aff)
974            img_back = round_trip(img)
975            arr_back_sc = img_back.get_data()
976            assert_true(np.allclose(arr, arr_back_sc))
977
978
979def test_float_int_spread():
980    # Test rounding error for spread of values
981    # Parallel test to arraywriters
982    powers = np.arange(-10, 10, 0.5)
983    arr = np.concatenate((-10**powers, 10**powers))
984    aff = np.eye(4)
985    for in_dt in (np.float32, np.float64):
986        arr_t = arr.astype(in_dt)
987        for out_dt in IUINT_TYPES:
988            img = Nifti1Image(arr_t, aff)
989            img_back = round_trip(img)
990            arr_back_sc = img_back.get_data()
991            slope, inter = img_back.get_header().get_slope_inter()
992            # Get estimate for error
993            max_miss = rt_err_estimate(arr_t, arr_back_sc.dtype, slope, inter)
994            # Simulate allclose test with large atol
995            diff = np.abs(arr_t - arr_back_sc)
996            rdiff = diff / np.abs(arr_t)
997            assert_true(np.all((diff <= max_miss) | (rdiff <= 1e-5)))
998
999
1000def test_rt_bias():
1001    # Check for bias in round trip
1002    # Parallel test to arraywriters
1003    rng = np.random.RandomState(20111214)
1004    mu, std, count = 100, 10, 100
1005    arr = rng.normal(mu, std, size=(count,))
1006    eps = np.finfo(np.float32).eps
1007    aff = np.eye(4)
1008    for in_dt in (np.float32, np.float64):
1009        arr_t = arr.astype(in_dt)
1010        for out_dt in IUINT_TYPES:
1011            img = Nifti1Image(arr_t, aff)
1012            img_back = round_trip(img)
1013            arr_back_sc = img_back.get_data()
1014            slope, inter = img_back.get_header().get_slope_inter()
1015            bias = np.mean(arr_t - arr_back_sc)
1016            # Get estimate for error
1017            max_miss = rt_err_estimate(arr_t, arr_back_sc.dtype, slope, inter)
1018            # Hokey use of max_miss as a std estimate
1019            bias_thresh = np.max([max_miss / np.sqrt(count), eps])
1020            assert_true(np.abs(bias) < bias_thresh)
1021