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 mghformat reading writing""" 10 11import os 12import io 13 14import numpy as np 15 16from .. import load, save 17from ...openers import ImageOpener 18from ..mghformat import MGHHeader, MGHError, MGHImage 19from ...tmpdirs import InTemporaryDirectory 20from ...fileholders import FileHolder 21from ...spatialimages import HeaderDataError 22from ...volumeutils import sys_is_le 23from ...wrapstruct import WrapStructError 24from ... import imageglobals 25 26 27import pytest 28 29from numpy.testing import assert_array_equal, assert_array_almost_equal, assert_almost_equal 30 31from ...testing import data_path 32 33from ...tests import test_spatialimages as tsi 34from ...tests.test_wrapstruct import _TestLabeledWrapStruct 35 36MGZ_FNAME = os.path.join(data_path, 'test.mgz') 37 38# sample voxel to ras matrix (mri_info --vox2ras) 39v2r = np.array([[1, 2, 3, -13], [2, 3, 1, -11.5], 40 [3, 1, 2, -11.5], [0, 0, 0, 1]], dtype=np.float32) 41# sample voxel to ras - tkr matrix (mri_info --vox2ras-tkr) 42v2rtkr = np.array([[-1.0, 0.0, 0.0, 1.5], 43 [0.0, 0.0, 1.0, -2.5], 44 [0.0, -1.0, 0.0, 2.0], 45 [0.0, 0.0, 0.0, 1.0]], dtype=np.float32) 46 47BIG_CODES = ('>', 'big', 'BIG', 'b', 'be', 'B', 'BE') 48LITTLE_CODES = ('<', 'little', 'l', 'le', 'L', 'LE') 49 50if sys_is_le: 51 BIG_CODES += ('swapped', 's', 'S', '!') 52 LITTLE_CODES += ('native', 'n', 'N', '=', '|', 'i', 'I') 53else: 54 BIG_CODES += ('native', 'n', 'N', '=', '|', 'i', 'I') 55 LITTLE_CODES += ('swapped', 's', 'S', '!') 56 57 58 59def test_read_mgh(): 60 # test.mgz was generated by the following command 61 # mri_volsynth --dim 3 4 5 2 --vol test.mgz 62 # --cdircos 1 2 3 --rdircos 2 3 1 --sdircos 3 1 2 63 # mri_volsynth is a FreeSurfer command 64 mgz = load(MGZ_FNAME) 65 66 # header 67 h = mgz.header 68 assert h['version'] == 1 69 assert h['type'] == 3 70 assert h['dof'] == 0 71 assert h['goodRASFlag'] == 1 72 assert_array_equal(h['dims'], [3, 4, 5, 2]) 73 assert_almost_equal(h['tr'], 2.0) 74 assert_almost_equal(h['flip_angle'], 0.0) 75 assert_almost_equal(h['te'], 0.0) 76 assert_almost_equal(h['ti'], 0.0) 77 assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 2]) 78 assert_array_almost_equal(h.get_vox2ras(), v2r) 79 assert_array_almost_equal(h.get_vox2ras_tkr(), v2rtkr) 80 81 # data. will be different for your own mri_volsynth invocation 82 v = mgz.get_fdata() 83 assert_almost_equal(v[1, 2, 3, 0], -0.3047, 4) 84 assert_almost_equal(v[1, 2, 3, 1], 0.0018, 4) 85 86 87def test_write_mgh(): 88 # write our data to a tmp file 89 v = np.arange(120) 90 v = v.reshape((5, 4, 3, 2)).astype(np.float32) 91 # form a MGHImage object using data and vox2ras matrix 92 img = MGHImage(v, v2r) 93 with InTemporaryDirectory(): 94 save(img, 'tmpsave.mgz') 95 # read from the tmp file and see if it checks out 96 mgz = load('tmpsave.mgz') 97 h = mgz.header 98 dat = mgz.get_fdata() 99 # Delete loaded image to allow file deletion by windows 100 del mgz 101 # header 102 assert h['version'] == 1 103 assert h['type'] == 3 104 assert h['dof'] == 0 105 assert h['goodRASFlag'] == 1 106 assert np.array_equal(h['dims'], [5, 4, 3, 2]) 107 assert_almost_equal(h['tr'], 0.0) 108 assert_almost_equal(h['flip_angle'], 0.0) 109 assert_almost_equal(h['te'], 0.0) 110 assert_almost_equal(h['ti'], 0.0) 111 assert_almost_equal(h['fov'], 0.0) 112 assert_array_almost_equal(h.get_vox2ras(), v2r) 113 # data 114 assert_almost_equal(dat, v, 7) 115 116 117def test_write_noaffine_mgh(): 118 # now just save the image without the vox2ras transform 119 # and see if it uses the default values to save 120 v = np.ones((7, 13, 3, 22)).astype(np.uint8) 121 # form a MGHImage object using data 122 # and the default affine matrix (Note the "None") 123 img = MGHImage(v, None) 124 with InTemporaryDirectory(): 125 save(img, 'tmpsave.mgz') 126 # read from the tmp file and see if it checks out 127 mgz = load('tmpsave.mgz') 128 h = mgz.header 129 # Delete loaded image to allow file deletion by windows 130 del mgz 131 # header 132 assert h['version'] == 1 133 assert h['type'] == 0 # uint8 for mgh 134 assert h['dof'] == 0 135 assert h['goodRASFlag'] == 1 136 assert np.array_equal(h['dims'], [7, 13, 3, 22]) 137 assert_almost_equal(h['tr'], 0.0) 138 assert_almost_equal(h['flip_angle'], 0.0) 139 assert_almost_equal(h['te'], 0.0) 140 assert_almost_equal(h['ti'], 0.0) 141 assert_almost_equal(h['fov'], 0.0) 142 # important part -- whether default affine info is stored 143 assert_array_almost_equal(h['Mdc'], [[-1, 0, 0], [0, 0, 1], [0, -1, 0]]) 144 assert_array_almost_equal(h['Pxyz_c'], [0, 0, 0]) 145 146 147def test_set_zooms(): 148 mgz = load(MGZ_FNAME) 149 h = mgz.header 150 assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 2]) 151 h.set_zooms([1, 1, 1, 3]) 152 assert_array_almost_equal(h.get_zooms(), [1, 1, 1, 3]) 153 for zooms in ((-1, 1, 1, 1), 154 (1, -1, 1, 1), 155 (1, 1, -1, 1), 156 (1, 1, 1, -1), 157 (1, 1, 1, 1, 5)): 158 with pytest.raises(HeaderDataError): 159 h.set_zooms(zooms) 160 # smoke test for tr=0 161 h.set_zooms((1, 1, 1, 0)) 162 163 164def bad_dtype_mgh(): 165 """ This function raises an MGHError exception because 166 uint16 is not a valid MGH datatype. 167 """ 168 # try to write an unsigned short and make sure it 169 # raises MGHError 170 v = np.ones((7, 13, 3, 22)).astype(np.uint16) 171 # form a MGHImage object using data 172 # and the default affine matrix (Note the "None") 173 MGHImage(v, None) 174 175 176def test_bad_dtype_mgh(): 177 # Now test the above function 178 with pytest.raises(MGHError): 179 bad_dtype_mgh() 180 181 182def test_filename_exts(): 183 # Test acceptable filename extensions 184 v = np.ones((7, 13, 3, 22)).astype(np.uint8) 185 # form a MGHImage object using data 186 # and the default affine matrix (Note the "None") 187 img = MGHImage(v, None) 188 # Check if these extensions allow round trip 189 for ext in ('.mgh', '.mgz'): 190 with InTemporaryDirectory(): 191 fname = 'tmpname' + ext 192 save(img, fname) 193 # read from the tmp file and see if it checks out 194 img_back = load(fname) 195 assert_array_equal(img_back.get_fdata(), v) 196 del img_back 197 198 199def _mgh_rt(img, fobj): 200 file_map = {'image': FileHolder(fileobj=fobj)} 201 img.to_file_map(file_map) 202 return MGHImage.from_file_map(file_map) 203 204 205def test_header_updating(): 206 # Don't update the header information if the affine doesn't change. 207 # Luckily the test.mgz dataset had a bad set of cosine vectors, so these 208 # will be changed if the affine gets updated 209 mgz = load(MGZ_FNAME) 210 hdr = mgz.header 211 # Test against mri_info output 212 exp_aff = np.loadtxt(io.BytesIO(b""" 213 1.0000 2.0000 3.0000 -13.0000 214 2.0000 3.0000 1.0000 -11.5000 215 3.0000 1.0000 2.0000 -11.5000 216 0.0000 0.0000 0.0000 1.0000""")) 217 assert_almost_equal(mgz.affine, exp_aff, 6) 218 assert_almost_equal(hdr.get_affine(), exp_aff, 6) 219 # Test that initial wonky header elements have not changed 220 assert np.all(hdr['delta'] == 1) 221 assert_almost_equal(hdr['Mdc'].T, exp_aff[:3, :3]) 222 # Save, reload, same thing 223 img_fobj = io.BytesIO() 224 mgz2 = _mgh_rt(mgz, img_fobj) 225 hdr2 = mgz2.header 226 assert_almost_equal(hdr2.get_affine(), exp_aff, 6) 227 assert_array_equal(hdr2['delta'],1) 228 # Change affine, change underlying header info 229 exp_aff_d = exp_aff.copy() 230 exp_aff_d[0, -1] = -14 231 # This will (probably) become part of the official API 232 mgz2._affine[:] = exp_aff_d 233 mgz2.update_header() 234 assert_almost_equal(hdr2.get_affine(), exp_aff_d, 6) 235 RZS = exp_aff_d[:3, :3] 236 assert_almost_equal(hdr2['delta'], np.sqrt(np.sum(RZS ** 2, axis=0))) 237 assert_almost_equal(hdr2['Mdc'].T, RZS / hdr2['delta']) 238 239 240def test_cosine_order(): 241 # Test we are interpreting the cosine order right 242 data = np.arange(60).reshape((3, 4, 5)).astype(np.int32) 243 aff = np.diag([2., 3, 4, 1]) 244 aff[0] = [2, 1, 0, 10] 245 img = MGHImage(data, aff) 246 assert_almost_equal(img.affine, aff, 6) 247 img_fobj = io.BytesIO() 248 img2 = _mgh_rt(img, img_fobj) 249 hdr2 = img2.header 250 RZS = aff[:3, :3] 251 zooms = np.sqrt(np.sum(RZS ** 2, axis=0)) 252 assert_almost_equal(hdr2['Mdc'].T, RZS / zooms) 253 assert_almost_equal(hdr2['delta'], zooms) 254 255 256def test_eq(): 257 # Test headers compare properly 258 hdr = MGHHeader() 259 hdr2 = MGHHeader() 260 assert hdr == hdr2 261 hdr.set_data_shape((2, 3, 4)) 262 assert(hdr != hdr2) 263 hdr2.set_data_shape((2, 3, 4)) 264 assert hdr == hdr2 265 266 267def test_header_slope_inter(): 268 # Test placeholder slope / inter method 269 hdr = MGHHeader() 270 assert hdr.get_slope_inter() == (None, None) 271 272 273def test_mgh_load_fileobj(): 274 # Checks the filename gets passed to array proxy 275 # 276 # This is a bit of an implementation detail, but the test is to make sure 277 # that we aren't passing ImageOpener objects to the array proxy, as these 278 # were confusing mmap on Python 3. If there's some sensible reason not to 279 # pass the filename to the array proxy, please feel free to change this 280 # test. 281 img = MGHImage.load(MGZ_FNAME) 282 assert img.dataobj.file_like == MGZ_FNAME 283 # Check fileobj also passed into dataobj 284 with ImageOpener(MGZ_FNAME) as fobj: 285 contents = fobj.read() 286 bio = io.BytesIO(contents) 287 fm = MGHImage.make_file_map(mapping=dict(image=bio)) 288 img2 = MGHImage.from_file_map(fm) 289 assert(img2.dataobj.file_like is bio) 290 assert_array_equal(img.get_fdata(), img2.get_fdata()) 291 292 293def test_mgh_affine_default(): 294 hdr = MGHHeader() 295 hdr['goodRASFlag'] = 0 296 hdr2 = MGHHeader(hdr.binaryblock) 297 assert hdr2['goodRASFlag'] == 1 298 assert_array_equal(hdr['Mdc'], hdr2['Mdc']) 299 assert_array_equal(hdr['Pxyz_c'], hdr2['Pxyz_c']) 300 301 302def test_mgh_set_data_shape(): 303 hdr = MGHHeader() 304 hdr.set_data_shape((5,)) 305 assert_array_equal(hdr.get_data_shape(), (5, 1, 1)) 306 hdr.set_data_shape((5, 4)) 307 assert_array_equal(hdr.get_data_shape(), (5, 4, 1)) 308 hdr.set_data_shape((5, 4, 3)) 309 assert_array_equal(hdr.get_data_shape(), (5, 4, 3)) 310 hdr.set_data_shape((5, 4, 3, 2)) 311 assert_array_equal(hdr.get_data_shape(), (5, 4, 3, 2)) 312 with pytest.raises(ValueError): 313 hdr.set_data_shape((5, 4, 3, 2, 1)) 314 315 316def test_mghheader_default_structarr(): 317 hdr = MGHHeader.default_structarr() 318 assert hdr['version'] == 1 319 assert_array_equal(hdr['dims'], 1) 320 assert hdr['type'] == 3 321 assert hdr['dof'] == 0 322 assert hdr['goodRASFlag'] == 1 323 assert_array_equal(hdr['delta'], 1) 324 assert_array_equal(hdr['Mdc'], [[-1, 0, 0], [0, 0, 1], [0, -1, 0]]) 325 assert_array_equal(hdr['Pxyz_c'], 0) 326 assert hdr['tr'] == 0 327 assert hdr['flip_angle'] == 0 328 assert hdr['te'] == 0 329 assert hdr['ti'] == 0 330 assert hdr['fov'] == 0 331 332 for endianness in (None,) + BIG_CODES: 333 hdr2 = MGHHeader.default_structarr(endianness=endianness) 334 assert hdr2 == hdr 335 assert hdr2.newbyteorder('>') == hdr 336 337 for endianness in LITTLE_CODES: 338 with pytest.raises(ValueError): 339 MGHHeader.default_structarr(endianness=endianness) 340 341 342def test_deprecated_fields(): 343 hdr = MGHHeader() 344 hdr_data = MGHHeader._HeaderData(hdr.structarr) 345 346 # mrparams is the only deprecated field at the moment 347 # Accessing hdr_data is equivalent to accessing hdr, so double all checks 348 with pytest.deprecated_call(match="from version: 2.3"): 349 assert_array_equal(hdr['mrparams'], 0) 350 assert_array_equal(hdr_data['mrparams'], 0) 351 352 with pytest.deprecated_call(match="from version: 2.3"): 353 hdr['mrparams'] = [1, 2, 3, 4] 354 with pytest.deprecated_call(match="from version: 2.3"): 355 assert_array_almost_equal(hdr['mrparams'], [1, 2, 3, 4]) 356 assert hdr['tr'] == 1 357 assert hdr['flip_angle'] == 2 358 assert hdr['te'] == 3 359 assert hdr['ti'] == 4 360 assert hdr['fov'] == 0 361 assert_array_almost_equal(hdr_data['mrparams'], [1, 2, 3, 4]) 362 assert hdr_data['tr'] == 1 363 assert hdr_data['flip_angle'] == 2 364 assert hdr_data['te'] == 3 365 assert hdr_data['ti'] == 4 366 assert hdr_data['fov'] == 0 367 368 hdr['tr'] = 5 369 hdr['flip_angle'] = 6 370 hdr['te'] = 7 371 hdr['ti'] = 8 372 with pytest.deprecated_call(match="from version: 2.3"): 373 assert_array_almost_equal(hdr['mrparams'], [5, 6, 7, 8]) 374 assert_array_almost_equal(hdr_data['mrparams'], [5, 6, 7, 8]) 375 376 hdr_data['tr'] = 9 377 hdr_data['flip_angle'] = 10 378 hdr_data['te'] = 11 379 hdr_data['ti'] = 12 380 with pytest.deprecated_call(match="from version: 2.3"): 381 assert_array_almost_equal(hdr['mrparams'], [9, 10, 11, 12]) 382 assert_array_almost_equal(hdr_data['mrparams'], [9, 10, 11, 12]) 383 384 385class TestMGHImage(tsi.TestSpatialImage, tsi.MmapImageMixin): 386 """ Apply general image tests to MGHImage 387 """ 388 image_class = MGHImage 389 can_save = True 390 391 def check_dtypes(self, expected, actual): 392 # Some images will want dtypes to be equal including endianness, 393 # others may only require the same type 394 # MGH requires the actual to be a big endian version of expected 395 assert expected.newbyteorder('>') == actual 396 397 398class TestMGHHeader(_TestLabeledWrapStruct): 399 header_class = MGHHeader 400 401 def _set_something_into_hdr(self, hdr): 402 hdr['dims'] = [4, 3, 2, 1] 403 404 def get_bad_bb(self): 405 return b'\xff' + b'\x00' * self.header_class._hdrdtype.itemsize 406 407 # Update tests to account for big-endian requirement 408 def test_general_init(self): 409 hdr = self.header_class() 410 # binaryblock has length given by header data dtype 411 binblock = hdr.binaryblock 412 assert len(binblock) == hdr.structarr.dtype.itemsize 413 # Endianness will always be big, and cannot be set 414 assert hdr.endianness == '>' 415 # You can also pass in a check flag, without data this has no 416 # effect 417 hdr = self.header_class(check=False) 418 419 def test__eq__(self): 420 # Test equal and not equal 421 hdr1 = self.header_class() 422 hdr2 = self.header_class() 423 assert hdr1 == hdr2 424 self._set_something_into_hdr(hdr1) 425 assert hdr1 != hdr2 426 self._set_something_into_hdr(hdr2) 427 assert hdr1 == hdr2 428 # REMOVED as_byteswapped() test 429 # Check comparing to funny thing says no 430 assert hdr1 != None 431 assert hdr1 != 1 432 433 def test_to_from_fileobj(self): 434 # Successful write using write_to 435 hdr = self.header_class() 436 str_io = io.BytesIO() 437 hdr.write_to(str_io) 438 str_io.seek(0) 439 hdr2 = self.header_class.from_fileobj(str_io) 440 assert hdr2.endianness == '>' 441 assert hdr2.binaryblock == hdr.binaryblock 442 443 def test_endian_guess(self): 444 # Check guesses of endian 445 eh = self.header_class() 446 assert eh.endianness == '>' 447 assert self.header_class.guessed_endian(eh) == '>' 448 449 def test_bytes(self): 450 # Test get of bytes 451 hdr1 = self.header_class() 452 bb = hdr1.binaryblock 453 hdr2 = self.header_class(hdr1.binaryblock) 454 assert hdr1 == hdr2 455 assert hdr1.binaryblock == hdr2.binaryblock 456 # Do a set into the header, and try again. The specifics of 'setting 457 # something' will depend on the nature of the bytes object 458 self._set_something_into_hdr(hdr1) 459 hdr2 = self.header_class(hdr1.binaryblock) 460 assert hdr1 == hdr2 461 assert hdr1.binaryblock == hdr2.binaryblock 462 # Short binaryblocks give errors (here set through init) 463 # Long binaryblocks are truncated 464 with pytest.raises(WrapStructError): 465 self.header_class(bb[:self.header_class._hdrdtype.itemsize - 1]) 466 467 # Checking set to true by default, and prevents nonsense being 468 # set into the header. 469 bb_bad = self.get_bad_bb() 470 if bb_bad is None: 471 return 472 with imageglobals.LoggingOutputSuppressor(): 473 with pytest.raises(HeaderDataError): 474 self.header_class(bb_bad) 475 476 # now slips past without check 477 _ = self.header_class(bb_bad, check=False) 478 479 def test_as_byteswapped(self): 480 # Check byte swapping 481 hdr = self.header_class() 482 assert hdr.endianness == '>' 483 # same code just returns a copy 484 for endianness in BIG_CODES: 485 hdr2 = hdr.as_byteswapped(endianness) 486 assert(hdr2 is not hdr) 487 assert hdr2 == hdr 488 489 # Different code raises error 490 for endianness in (None,) + LITTLE_CODES: 491 with pytest.raises(ValueError): 492 hdr.as_byteswapped(endianness) 493 # Note that contents is not rechecked on swap / copy 494 class DC(self.header_class): 495 def check_fix(self, *args, **kwargs): 496 raise Exception 497 498 # Assumes check=True default 499 with pytest.raises(Exception): 500 DC(hdr.binaryblock) 501 502 hdr = DC(hdr.binaryblock, check=False) 503 hdr2 = hdr.as_byteswapped('>') 504 505 def test_checks(self): 506 # Test header checks 507 hdr_t = self.header_class() 508 # _dxer just returns the diagnostics as a string 509 # Default hdr is OK 510 assert self._dxer(hdr_t) == '' 511 # Version should be 1 512 hdr = hdr_t.copy() 513 hdr['version'] = 2 514 assert self._dxer(hdr) == 'Unknown MGH format version' 515