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