# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
#
# See COPYING file distributed along with the NiBabel package for the
# copyright and license terms.
#
### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
''' Test Analyze headers
See test_wrapstruct.py for tests of the wrapped structarr-ness of the Analyze
header
'''
import os
import re
import logging
import pickle
import numpy as np
from ..py3k import BytesIO, StringIO, asbytes
from ..volumeutils import array_to_file
from ..spatialimages import (HeaderDataError, HeaderTypeError)
from ..analyze import AnalyzeHeader, AnalyzeImage
from ..nifti1 import Nifti1Header
from ..loadsave import read_img_data
from .. import imageglobals
from ..casting import as_int
from numpy.testing import (assert_array_equal,
assert_array_almost_equal)
from ..testing import (assert_equal, assert_not_equal, assert_true,
assert_false, assert_raises, data_path)
from .test_wrapstruct import _TestWrapStructBase
from . import test_spatialimages as tsi
header_file = os.path.join(data_path, 'analyze.hdr')
PIXDIM0_MSG = 'pixdim[1,2,3] should be non-zero; setting 0 dims to 1'
def _write_data(hdr, data, fileobj):
# auxilary function to write data
out_dtype = hdr.get_data_dtype()
offset = hdr.get_data_offset()
array_to_file(data, fileobj, out_dtype, offset)
class TestAnalyzeHeader(_TestWrapStructBase):
header_class = AnalyzeHeader
example_file = header_file
def test_general_init(self):
super(TestAnalyzeHeader, self).test_general_init()
hdr = self.header_class()
# an empty header has shape (0,) - like an empty array
# (np.array([]))
assert_equal(hdr.get_data_shape(), (0,))
# The affine is always homogenous 3D regardless of shape. The
# default affine will have -1 as the X zoom iff default_x_flip
# is True (which it is by default). We have to be careful of the
# translations though - these arise from SPM's use of the origin
# field, and the center of the image.
assert_array_equal(np.diag(hdr.get_base_affine()),
[-1,1,1,1])
# But zooms only go with number of dimensions
assert_equal(hdr.get_zooms(), (1.0,))
def test_header_size(self):
assert_equal(self.header_class.template_dtype.itemsize, 348)
def test_empty(self):
hdr = self.header_class()
assert_true(len(hdr.binaryblock) == 348)
assert_true(hdr['sizeof_hdr'] == 348)
assert_true(np.all(hdr['dim'][1:] == 1))
assert_true(hdr['dim'][0] == 0 )
assert_true(np.all(hdr['pixdim'] == 1))
assert_true(hdr['datatype'] == 16) # float32
assert_true(hdr['bitpix'] == 32)
def _set_something_into_hdr(self, hdr):
# Called from test_bytes test method. Specific to the header data type
hdr.set_data_shape((1, 2, 3))
def test_checks(self):
# Test header checks
hdr_t = self.header_class()
# _dxer just returns the diagnostics as a string
assert_equal(self._dxer(hdr_t), '')
hdr = hdr_t.copy()
hdr['sizeof_hdr'] = 1
assert_equal(self._dxer(hdr), 'sizeof_hdr should be 348')
hdr = hdr_t.copy()
hdr['datatype'] = 0
assert_equal(self._dxer(hdr), 'data code 0 not supported\n'
'bitpix does not match datatype')
hdr = hdr_t.copy()
hdr['bitpix'] = 0
assert_equal(self._dxer(hdr), 'bitpix does not match datatype')
for i in (1,2,3):
hdr = hdr_t.copy()
hdr['pixdim'][i] = -1
assert_equal(self._dxer(hdr), 'pixdim[1,2,3] should be positive')
def test_log_checks(self):
# Test logging, fixing, errors for header checking
HC = self.header_class
# magic
hdr = HC()
hdr['sizeof_hdr'] = 350 # severity 30
fhdr, message, raiser = self.log_chk(hdr, 30)
assert_equal(fhdr['sizeof_hdr'], 348)
assert_equal(message, 'sizeof_hdr should be 348; '
'set sizeof_hdr to 348')
assert_raises(*raiser)
# RGB datatype does not raise error
hdr = HC()
hdr.set_data_dtype('RGB')
fhdr, message, raiser = self.log_chk(hdr, 0)
# datatype not recognized
hdr = HC()
hdr['datatype'] = -1 # severity 40
fhdr, message, raiser = self.log_chk(hdr, 40)
assert_equal(message, 'data code -1 not recognized; '
'not attempting fix')
assert_raises(*raiser)
# datatype not supported
hdr['datatype'] = 255 # severity 40
fhdr, message, raiser = self.log_chk(hdr, 40)
assert_equal(message, 'data code 255 not supported; '
'not attempting fix')
assert_raises(*raiser)
# bitpix
hdr = HC()
hdr['datatype'] = 16 # float32
hdr['bitpix'] = 16 # severity 10
fhdr, message, raiser = self.log_chk(hdr, 10)
assert_equal(fhdr['bitpix'], 32)
assert_equal(message, 'bitpix does not match datatype; '
'setting bitpix to match datatype')
assert_raises(*raiser)
# pixdim positive
hdr = HC()
hdr['pixdim'][1] = -2 # severity 35
fhdr, message, raiser = self.log_chk(hdr, 35)
assert_equal(fhdr['pixdim'][1], 2)
assert_equal(message, 'pixdim[1,2,3] should be positive; '
'setting to abs of pixdim values')
assert_raises(*raiser)
hdr = HC()
hdr['pixdim'][1] = 0 # severity 30
fhdr, message, raiser = self.log_chk(hdr, 30)
assert_equal(fhdr['pixdim'][1], 1)
assert_equal(message, PIXDIM0_MSG)
assert_raises(*raiser)
# both
hdr = HC()
hdr['pixdim'][1] = 0 # severity 30
hdr['pixdim'][2] = -2 # severity 35
fhdr, message, raiser = self.log_chk(hdr, 35)
assert_equal(fhdr['pixdim'][1], 1)
assert_equal(fhdr['pixdim'][2], 2)
assert_equal(message, 'pixdim[1,2,3] should be '
'non-zero and pixdim[1,2,3] should '
'be positive; setting 0 dims to 1 '
'and setting to abs of pixdim values')
assert_raises(*raiser)
def test_logger_error(self):
# Check that we can reset the logger and error level
HC = self.header_class
hdr = HC()
# Make a new logger
str_io = StringIO()
logger = logging.getLogger('test.logger')
logger.setLevel(30) # defaultish level
logger.addHandler(logging.StreamHandler(str_io))
# Prepare an error
hdr['pixdim'][1] = 0 # severity 30
log_cache = imageglobals.logger, imageglobals.error_level
try:
# Check log message appears in new logger
imageglobals.logger = logger
hdr.copy().check_fix()
assert_equal(str_io.getvalue(), PIXDIM0_MSG + '\n')
# Check that error_level in fact causes error to be raised
imageglobals.error_level = 30
assert_raises(HeaderDataError, hdr.copy().check_fix)
finally:
imageglobals.logger, imageglobals.error_level = log_cache
def test_data_dtype(self):
# check getting and setting of data type
# codes / types supported by all binary headers
supported_types = ((2, np.uint8),
(4, np.int16),
(8, np.int32),
(16, np.float32),
(32, np.complex64),
(64, np.float64),
(128, np.dtype([('R','u1'),
('G', 'u1'),
('B', 'u1')])))
# and unsupported - here using some labels instead
unsupported_types = (np.void, 'none', 'all', 0)
hdr = self.header_class()
for code, npt in supported_types:
# Can set with code value, or numpy dtype, both return the
# dtype as output on get
hdr.set_data_dtype(code)
assert_equal(hdr.get_data_dtype(), npt)
hdr.set_data_dtype(npt)
assert_equal(hdr.get_data_dtype(), npt)
for inp in unsupported_types:
assert_raises(HeaderDataError,
hdr.set_data_dtype,
inp)
def test_shapes(self):
# Test that shape checks work
hdr = self.header_class()
for shape in ((2, 3, 4), (2, 3, 4, 5), (2, 3), (2,)):
hdr.set_data_shape(shape)
assert_equal(hdr.get_data_shape(), shape)
# Check max works, but max+1 raises error
dim_dtype = hdr.structarr['dim'].dtype
# as_int for safety to deal with numpy 1.4.1 int conversion errors
mx = as_int(np.iinfo(dim_dtype).max)
shape = (mx,)
hdr.set_data_shape(shape)
assert_equal(hdr.get_data_shape(), shape)
shape = (mx+1,)
assert_raises(HeaderDataError, hdr.set_data_shape, shape)
# Lists or tuples or arrays will work for setting shape
shape = (2, 3, 4)
for constructor in (list, tuple, np.array):
hdr.set_data_shape(constructor(shape))
assert_equal(hdr.get_data_shape(), shape)
def test_read_write_data(self):
# Check reading and writing of data
hdr = self.header_class()
# Trying to read data from an empty header gives no data
bytes = hdr.data_from_fileobj(BytesIO())
assert_equal(len(bytes), 0)
# Setting no data into an empty header results in - no data
str_io = BytesIO()
hdr.data_to_fileobj([], str_io)
assert_equal(str_io.getvalue(), asbytes(''))
# Setting more data then there should be gives an error
assert_raises(HeaderDataError,
hdr.data_to_fileobj,
np.zeros(3),
str_io)
# Test valid write
hdr.set_data_shape((1,2,3))
hdr.set_data_dtype(np.float32)
S = BytesIO()
data = np.arange(6, dtype=np.float64)
# data have to be the right shape
assert_raises(HeaderDataError, hdr.data_to_fileobj, data, S)
data = data.reshape((1,2,3))
# and size
assert_raises(HeaderDataError, hdr.data_to_fileobj, data[:,:,:-1], S)
assert_raises(HeaderDataError, hdr.data_to_fileobj, data[:,:-1,:], S)
# OK if so
hdr.data_to_fileobj(data, S)
# Read it back
data_back = hdr.data_from_fileobj(S)
# Should be about the same
assert_array_almost_equal(data, data_back)
# but with the header dtype, not the data dtype
assert_equal(hdr.get_data_dtype(), data_back.dtype)
# this is with native endian, not so for swapped
S2 = BytesIO()
hdr2 = hdr.as_byteswapped()
hdr2.set_data_dtype(np.float32)
hdr2.set_data_shape((1,2,3))
hdr2.data_to_fileobj(data, S2)
data_back2 = hdr2.data_from_fileobj(S2)
# Compares the same
assert_array_almost_equal(data_back, data_back2)
# Same dtype names
assert_equal(data_back.dtype.name, data_back2.dtype.name)
# But not the same endianness
assert_not_equal(data.dtype.byteorder, data_back2.dtype.byteorder)
# Try scaling down to integer
hdr.set_data_dtype(np.uint8)
S3 = BytesIO()
# Analyze header cannot do scaling, but, if not scaling, AnalyzeHeader
# is OK
_write_data(hdr, data, S3)
data_back = hdr.data_from_fileobj(S3)
assert_array_almost_equal(data, data_back)
# But, the data won't always be same as input if not scaling
data = np.arange(6, dtype=np.float64).reshape((1,2,3)) + 0.5
_write_data(hdr, data, S3)
data_back = hdr.data_from_fileobj(S3)
assert_false(np.allclose(data, data_back))
# Test RGB image
dtype = np.dtype([('R', 'uint8'), ('G', 'uint8'), ('B', 'uint8')])
data = np.ones((1, 2, 3), dtype)
hdr.set_data_dtype(dtype)
S4 = BytesIO()
hdr.data_to_fileobj(data, S4)
data_back = hdr.data_from_fileobj(S4)
assert_array_equal(data, data_back)
def test_datatype(self):
ehdr = self.header_class()
codes = self.header_class._data_type_codes
for code in codes.value_set():
npt = codes.type[code]
if npt is np.void:
assert_raises(
HeaderDataError,
ehdr.set_data_dtype,
code)
continue
dt = codes.dtype[code]
ehdr.set_data_dtype(npt)
assert_true(ehdr['datatype'] == code)
assert_true(ehdr['bitpix'] == dt.itemsize*8)
ehdr.set_data_dtype(code)
assert_true(ehdr['datatype'] == code)
ehdr.set_data_dtype(dt)
assert_true(ehdr['datatype'] == code)
def test_data_shape_zooms_affine(self):
hdr = self.header_class()
for shape in ((1,2,3),(0,),(1,),(1,2),(1,2,3,4)):
L = len(shape)
hdr.set_data_shape(shape)
if L:
assert_equal(hdr.get_data_shape(), shape)
else:
assert_equal(hdr.get_data_shape(), (0,))
# Default zoom - for 3D - is 1(())
assert_equal(hdr.get_zooms(), (1,) * L)
# errors if zooms do not match shape
if len(shape):
assert_raises(HeaderDataError,
hdr.set_zooms,
(1,) * (L-1))
# Errors for negative zooms
assert_raises(HeaderDataError,
hdr.set_zooms,
(-1,) + (1,)*(L-1))
assert_raises(HeaderDataError,
hdr.set_zooms,
(1,) * (L+1))
# Errors for negative zooms
assert_raises(HeaderDataError,
hdr.set_zooms,
(-1,) * L)
# reducing the dimensionality of the array and then increasing
# it again reverts the previously set zoom values to 1.0
hdr = self.header_class()
hdr.set_data_shape((1,2,3))
hdr.set_zooms((4,5,6))
assert_array_equal(hdr.get_zooms(), (4,5,6))
hdr.set_data_shape((1,2))
assert_array_equal(hdr.get_zooms(), (4,5))
hdr.set_data_shape((1,2,3))
assert_array_equal(hdr.get_zooms(), (4,5,1))
# Setting zooms changes affine
assert_array_equal(np.diag(hdr.get_base_affine()),
[-4,5,1,1])
hdr.set_zooms((1,1,1))
assert_array_equal(np.diag(hdr.get_base_affine()),
[-1,1,1,1])
def test_default_x_flip(self):
hdr = self.header_class()
hdr.default_x_flip = True
hdr.set_data_shape((1,2,3))
hdr.set_zooms((1,1,1))
assert_array_equal(np.diag(hdr.get_base_affine()),
[-1,1,1,1])
hdr.default_x_flip = False
# Check avoids translations
assert_array_equal(np.diag(hdr.get_base_affine()),
[1,1,1,1])
def test_from_eg_file(self):
fileobj = open(self.example_file, 'rb')
hdr = self.header_class.from_fileobj(fileobj, check=False)
assert_equal(hdr.endianness, '>')
assert_equal(hdr['sizeof_hdr'], 348)
def test_orientation(self):
# Test flips
hdr = self.header_class()
assert_true(hdr.default_x_flip)
hdr.set_data_shape((3,5,7))
hdr.set_zooms((4,5,6))
aff = np.diag((-4,5,6,1))
aff[:3,3] = np.array([1,2,3]) * np.array([-4,5,6]) * -1
assert_array_equal(hdr.get_base_affine(), aff)
hdr.default_x_flip = False
assert_false(hdr.default_x_flip)
aff[0]*=-1
assert_array_equal(hdr.get_base_affine(), aff)
def test_str(self):
super(TestAnalyzeHeader, self).test_str()
hdr = self.header_class()
s1 = str(hdr)
# check the datacode recoding
rexp = re.compile('^datatype +: float32', re.MULTILINE)
assert_true(rexp.search(s1) is not None)
def test_from_header(self):
# check from header class method.
klass = self.header_class
empty = klass.from_header()
assert_equal(klass(), empty)
empty = klass.from_header(None)
assert_equal(klass(), empty)
hdr = klass()
hdr.set_data_dtype(np.float64)
hdr.set_data_shape((1,2,3))
hdr.set_zooms((3.0, 2.0, 1.0))
copy = klass.from_header(hdr)
assert_equal(hdr, copy)
assert_false(hdr is copy)
class C(object):
def get_data_dtype(self): return np.dtype('i2')
def get_data_shape(self): return (5,4,3)
def get_zooms(self): return (10.0, 9.0, 8.0)
converted = klass.from_header(C())
assert_true(isinstance(converted, klass))
assert_equal(converted.get_data_dtype(), np.dtype('i2'))
assert_equal(converted.get_data_shape(), (5,4,3))
assert_equal(converted.get_zooms(), (10.0,9.0,8.0))
def test_base_affine(self):
klass = self.header_class
hdr = klass()
hdr.set_data_shape((3, 5, 7))
hdr.set_zooms((3, 2, 1))
assert_true(hdr.default_x_flip)
assert_array_almost_equal(
hdr.get_base_affine(),
[[-3., 0., 0., 3.],
[ 0., 2., 0., -4.],
[ 0., 0., 1., -3.],
[ 0., 0., 0., 1.]])
hdr.set_data_shape((3, 5))
assert_array_almost_equal(
hdr.get_base_affine(),
[[-3., 0., 0., 3.],
[ 0., 2., 0., -4.],
[ 0., 0., 1., -0.],
[ 0., 0., 0., 1.]])
hdr.set_data_shape((3, 5, 7))
assert_array_almost_equal(
hdr.get_base_affine(),
[[-3., 0., 0., 3.],
[ 0., 2., 0., -4.],
[ 0., 0., 1., -3.],
[ 0., 0., 0., 1.]])
def test_best_affine():
hdr = AnalyzeHeader()
hdr.set_data_shape((3,5,7))
hdr.set_zooms((4,5,6))
assert_array_equal(hdr.get_base_affine(), hdr.get_best_affine())
def test_scaling():
# Test integer scaling from float
# Analyze headers cannot do float-integer scaling '''
hdr = AnalyzeHeader()
assert_true(hdr.default_x_flip)
shape = (1,2,3)
hdr.set_data_shape(shape)
hdr.set_data_dtype(np.float32)
data = np.ones(shape, dtype=np.float64)
S = BytesIO()
# Writing to float datatype doesn't need scaling
hdr.data_to_fileobj(data, S)
rdata = hdr.data_from_fileobj(S)
assert_array_almost_equal(data, rdata)
# Writing to integer datatype does, and raises an error
hdr.set_data_dtype(np.int32)
assert_raises(HeaderTypeError, hdr.data_to_fileobj, data, BytesIO())
# unless we aren't scaling, in which case we convert the floats to
# integers and write
_write_data(hdr, data, S)
rdata = hdr.data_from_fileobj(S)
assert_true(np.allclose(data, rdata))
# This won't work for floats that aren't close to integers
data_p5 = data + 0.5
_write_data(hdr, data_p5, S)
rdata = hdr.data_from_fileobj(S)
assert_false(np.allclose(data_p5, rdata))
def test_slope_inter():
hdr = AnalyzeHeader()
assert_equal(hdr.get_slope_inter(), (None, None))
for slinter in ((None,),
(None, None),
(1.0,),
(1.0, None),
(None, 0),
(1.0, 0)):
hdr.set_slope_inter(*slinter)
assert_equal(hdr.get_slope_inter(), (None, None))
assert_raises(HeaderTypeError, hdr.set_slope_inter, 1.1)
assert_raises(HeaderTypeError, hdr.set_slope_inter, 1.0, 0.1)
def test_data_code_error():
# test analyze raising error for unsupported codes
hdr = Nifti1Header()
hdr['datatype'] = 256
assert_raises(HeaderDataError, AnalyzeHeader.from_header, hdr)
class TestAnalyzeImage(tsi.TestSpatialImage):
image_class = AnalyzeImage
def test_data_hdr_cache(self):
# test the API for loaded images, such that the data returned
# from img.get_data() is not affected by subsequent changes to
# the header.
IC = self.image_class
# save an image to a file map
fm = IC.make_file_map()
for key, value in fm.items():
fm[key].fileobj = BytesIO()
shape = (2,3,4)
data = np.arange(24, dtype=np.int8).reshape(shape)
affine = np.eye(4)
hdr = IC.header_class()
hdr.set_data_dtype(np.int16)
img = IC(data, affine, hdr)
img.to_file_map(fm)
img2 = IC.from_file_map(fm)
assert_equal(img2.shape, shape)
assert_equal(img2.get_data_dtype().type, np.int16)
hdr = img2.get_header()
hdr.set_data_shape((3,2,2))
assert_equal(hdr.get_data_shape(), (3,2,2))
hdr.set_data_dtype(np.uint8)
assert_equal(hdr.get_data_dtype(), np.dtype(np.uint8))
assert_array_equal(img2.get_data(), data)
# now check read_img_data function - here we do see the changed
# header
sc_data = read_img_data(img2)
assert_equal(sc_data.shape, (3,2,2))
us_data = read_img_data(img2, prefer='unscaled')
assert_equal(us_data.shape, (3,2,2))
def test_affine_44(self):
IC = self.image_class
shape = (2,3,4)
data = np.arange(24, dtype=np.int16).reshape(shape)
affine = np.diag([2, 3, 4, 1])
# OK - affine correct shape
img = IC(data, affine)
assert_array_equal(affine, img.get_affine())
# OK - affine can be array-like
img = IC(data, affine.tolist())
assert_array_equal(affine, img.get_affine())
# Not OK - affine wrong shape
assert_raises(ValueError, IC, data, np.diag([2, 3, 4]))
def test_header_updating(self):
# Only update on changes
img_klass = self.image_class
# With a None affine - don't overwrite zooms
img = img_klass(np.zeros((2,3,4)), None)
hdr = img.get_header()
hdr.set_zooms((4,5,6))
# Save / reload using bytes IO objects
for key, value in img.file_map.items():
value.fileobj = BytesIO()
img.to_file_map()
hdr_back = img.from_file_map(img.file_map).get_header()
assert_array_equal(hdr_back.get_zooms(), (4,5,6))
# With a real affine, update zooms
img = img_klass(np.zeros((2,3,4)), np.diag([2,3,4,1]), hdr)
hdr = img.get_header()
assert_array_equal(hdr.get_zooms(), (2, 3, 4))
# Modify affine in-place? Update on save.
img.get_affine()[0,0] = 9
for key, value in img.file_map.items():
value.fileobj = BytesIO()
img.to_file_map()
hdr_back = img.from_file_map(img.file_map).get_header()
assert_array_equal(hdr.get_zooms(), (9, 3, 4))
# Modify data in-place? Update on save
data = img.get_data()
data.shape = (3, 2, 4)
img.to_file_map()
img_back = img.from_file_map(img.file_map)
assert_array_equal(img_back.shape, (3, 2, 4))
def test_pickle(self):
# Test that images pickle
# Image that is not proxied can pickle
img_klass = self.image_class
img = img_klass(np.zeros((2,3,4)), None)
img_str = pickle.dumps(img)
img2 = pickle.loads(img_str)
assert_array_equal(img.get_data(), img2.get_data())
assert_equal(img.get_header(), img2.get_header())
# Save / reload using bytes IO objects
for key, value in img.file_map.items():
value.fileobj = BytesIO()
img.to_file_map()
img_prox = img.from_file_map(img.file_map)
img_str = pickle.dumps(img_prox)
img2_prox = pickle.loads(img_str)
assert_array_equal(img.get_data(), img2_prox.get_data())
def test_unsupported():
# analyze does not support uint32
data = np.arange(24, dtype=np.int32).reshape((2,3,4))
affine = np.eye(4)
data = np.arange(24, dtype=np.uint32).reshape((2,3,4))
assert_raises(HeaderDataError, AnalyzeImage, data, affine)
