Commit 90fcc75d authored by Paul McCarthy's avatar Paul McCarthy 🚵
Browse files

Merge branch 'rf/nibabel-get_data' into 'master'

Rf/nibabel get data

See merge request fsl/fslpy!184
parents 3c4c5bdd 70216109
Pipeline #4840 failed with stages
in 19 minutes and 31 seconds
......@@ -17,6 +17,16 @@ Added
* New :func:`.affine.rescale` function, for adjusting a scaling matrix.
Changed
^^^^^^^
* The :class:`.ImageWrapper` now maintains its own image data cache, rather
than depending on ``nibabel``.
* Internal changes to avoid using the deprecated
``nibabel.dataobj_images.DataobjImage.get_data`` method.
Fixed
^^^^^
......
......@@ -36,6 +36,7 @@ import shlex
import logging
import binascii
import numpy as np
import nibabel as nib
import fsl.utils.tempdir as tempdir
......@@ -279,9 +280,11 @@ def loadSeries(series):
# copy images so nibabel no longer
# refs to the files (as they will
# be deleted), and use get_data()
# to force-load the image data.
images = [nib.Nifti1Image(i.get_data(), None, i.header)
# be deleted), and force-load the
# the image data into memory (to
# avoid any disk accesses due to
# e.g. memmap)
images = [nib.Nifti1Image(np.asanyarray(i.dataobj), None, i.header)
for i in images]
return [DicomImage(i, series, dcmdir, name=desc) for i in images]
......@@ -1204,6 +1204,9 @@ class Image(Nifti):
@property
def nibImage(self):
"""Returns a reference to the ``nibabel`` NIFTI image instance.
Note that if the image data has been modified through this ``Image``,
it will be out of sync with what is returned by the ``nibabel`` object,
until a call to :meth:`save` is made.
"""
return self.__nibImage
......@@ -1387,6 +1390,21 @@ class Image(Nifti):
os.close(tmphd)
try:
# First of all, the nibabel object won't know
# about any image data modifications, so if
# any have occurred, we need to create a new
# nibabel image using the data managed by the
# imagewrapper, and the old header.
#
# Assuming here that analyze/nifti1/nifti2
# nibabel classes have an __init__ which
# expects (data, affine, header)
if not self.saveState:
self.__nibImage = type(self.__nibImage)(self[:],
None,
self.header)
self.header = self.__nibImage.header
nib.save(self.__nibImage, tmpfname)
# nibabel should close any old
......@@ -1395,6 +1413,8 @@ class Image(Nifti):
self.__nibImage = None
self.header = None
# Copy to final destination,
# and reload from there
imcp.imcp(tmpfname, filename, overwrite=True)
self.__nibImage = nib.load(filename)
......
......@@ -76,8 +76,7 @@ class ImageWrapper(notifier.Notifier):
- The image data is not modified (via :meth:`__setitem__`.
If any of these conditions do not hold, the image data will be loaded into
memory and accessed directly, via the ``nibabel.Nifti1Image.get_data``
method.
memory and accessed directly.
*Image dimensionality*
......@@ -218,7 +217,11 @@ class ImageWrapper(notifier.Notifier):
self.reset(dataRange)
if loadData:
# We keep an internal ref to
# the data numpy array if/when
# it is loaded in memory
self.__data = None
if loadData or image.in_memory:
self.loadData()
if threaded:
......@@ -232,6 +235,7 @@ class ImageWrapper(notifier.Notifier):
the :class:`.TaskThread` is stopped.
"""
self.__image = None
self.__data = None
if self.__taskThread is not None:
self.__taskThread.stop()
self.__taskThread = None
......@@ -377,16 +381,11 @@ class ImageWrapper(notifier.Notifier):
"""Forces all of the image data to be loaded into memory.
.. note:: This method will be called by :meth:`__init__` if its
``loadData`` parameter is ``True``.
``loadData`` parameter is ``True``. It will also be called
on all write operations (see :meth:`__setitem__`).
"""
# If the data is not already loaded, this will
# cause nibabel to load it. By default, nibabel
# will cache the numpy array that contains the
# image data, so subsequent calls to this
# method will not overwrite any changes that
# have been made to the data array.
self.__image.get_data()
if self.__data is None:
self.__data = np.asanyarray(self.__image.dataobj)
def __getData(self, sliceobj, isTuple=False):
......@@ -407,14 +406,7 @@ class ImageWrapper(notifier.Notifier):
# ArrayProxy. Otheriwse if it is in
# memory, we can access it directly.
#
# Furthermore, if it is in memory and
# has been modified, the ArrayProxy
# will give us out-of-date values (as
# the ArrayProxy reads from disk). So
# we have to read from the in-memory
# array to get changed values.
#
# Finally, note that if the caller has
# Note also that if the caller has
# given us a 'fancy' slice object (a
# boolean numpy array), but the image
# data is not in memory, we can't access
......@@ -422,8 +414,8 @@ class ImageWrapper(notifier.Notifier):
# (the dataobj attribute) cannot handle
# fancy indexing. In this case an error
# will be raised.
if self.__image.in_memory: return self.__image.get_data()[sliceobj]
else: return self.__image.dataobj[ sliceobj]
if self.__data is not None: return self.__data[ sliceobj]
else: return self.__image.dataobj[sliceobj]
def __imageIsCovered(self):
......@@ -444,18 +436,18 @@ class ImageWrapper(notifier.Notifier):
_, expansions = calcExpansion(slices, self.__coverage)
expansions = collapseExpansions(expansions, self.__numRealDims - 1)
log.debug('Updating image {} data range [slice: {}] '
'(current range: [{}, {}]; '
'number of expansions: {}; '
'current coverage: {}; '
'volume ranges: {})'.format(
self.__name,
slices,
self.__range[0],
self.__range[1],
len(expansions),
self.__coverage,
self.__volRanges))
log.debug('Updating image %s data range [slice: %s] '
'(current range: [%s, %s]; '
'number of expansions: %s; '
'current coverage: %s; '
'volume ranges: %s)',
self.__name,
slices,
self.__range[0],
self.__range[1],
len(expansions),
self.__coverage,
self.__volRanges)
# As we access the data for each expansions,
# we want it to have the same dimensionality
......@@ -507,12 +499,12 @@ class ImageWrapper(notifier.Notifier):
if any((oldmin is None, oldmax is None)) or \
not np.all(np.isclose([oldmin, oldmax], [newmin, newmax])):
log.debug('Image {} range changed: [{}, {}] -> [{}, {}]'.format(
self.__name,
oldmin,
oldmax,
newmin,
newmax))
log.debug('Image %s range changed: [%s, %s] -> [%s, %s]',
self.__name,
oldmin,
oldmax,
newmin,
newmax)
self.notify()
......@@ -591,15 +583,15 @@ class ImageWrapper(notifier.Notifier):
slices = np.array(slices.T, dtype=np.uint32)
slices = tuple(it.chain(map(tuple, slices), [(lowVol, highVol)]))
log.debug('Image {} data written - clearing known data '
'range on volumes {} - {} (write slice: {}; '
'coverage: {}; volRanges: {})'.format(
self.__name,
lowVol,
highVol,
slices,
self.__coverage[:, :, lowVol:highVol],
self.__volRanges[lowVol:highVol, :]))
log.debug('Image %s data written - clearing known data '
'range on volumes %s - %s (write slice: %s; '
'coverage: %s; volRanges: %s)',
self.__name,
lowVol,
highVol,
slices,
self.__coverage[:, :, lowVol:highVol],
self.__volRanges[lowVol:highVol, :])
for vol in range(lowVol, highVol):
self.__coverage[:, :, vol] = np.nan
......@@ -622,7 +614,7 @@ class ImageWrapper(notifier.Notifier):
:arg sliceobj: Something which can slice the image data.
"""
log.debug('Getting image data: {}'.format(sliceobj))
log.debug('Getting image data: %s', sliceobj)
shape = self.__canonicalShape
realShape = self.__image.shape
......@@ -721,7 +713,7 @@ class ImageWrapper(notifier.Notifier):
# have any effect.
self.loadData()
self.__image.get_data()[sliceobj] = values
self.__data[sliceobj] = values
self.__updateDataRangeOnWrite(slices, values)
......
......@@ -12,6 +12,7 @@ Freesurfer ``mgh``/``mgz`` image files.
import os.path as op
import six
import numpy as np
import nibabel as nib
import fsl.utils.path as fslpath
......@@ -54,7 +55,7 @@ class MGHImage(fslimage.Image):
name = 'MGH image'
filename = None
data = image.get_data()
data = np.asanyarray(image.dataobj)
xform = image.affine
vox2surf = image.header.get_vox2ras_tkr()
......
......@@ -952,16 +952,17 @@ def fileOrImage(*args, **kwargs):
# create an independent in-memory
# copy of the image file
img = nib.load(path, mmap=False)
img = nib.load(path, mmap=False)
data = np.asanyarray(img.dataobj)
# if any arguments were fsl images,
# that takes precedence.
if fslimage.Image in intypes:
return fslimage.Image(img.get_data(), header=img.header)
return fslimage.Image(data, header=img.header)
# but if all inputs were file names,
# nibabel takes precedence
elif nib.nifti1.Nifti1Image in intypes or len(intypes) == 0:
return nib.nifti1.Nifti1Image(img.get_data(), None, img.header)
return nib.nifti1.Nifti1Image(data, None, img.header)
# this function should not be called
# under any other circumstances
......
......@@ -169,12 +169,12 @@ def test_DTIFitTensor():
l2file = op.join(testdir, 'dti_L2.nii')
l3file = op.join(testdir, 'dti_L3.nii')
v1 = tests.make_random_image(v1file, (5, 5, 5, 3)).get_data()
v2 = tests.make_random_image(v2file, (5, 5, 5, 3)).get_data()
v3 = tests.make_random_image(v3file, (5, 5, 5, 3)).get_data()
l1 = tests.make_random_image(l1file, (5, 5, 5)) .get_data()
l2 = tests.make_random_image(l2file, (5, 5, 5)) .get_data()
l3 = tests.make_random_image(l3file, (5, 5, 5)) .get_data()
v1 = np.asanyarray(tests.make_random_image(v1file, (5, 5, 5, 3)).dataobj)
v2 = np.asanyarray(tests.make_random_image(v2file, (5, 5, 5, 3)).dataobj)
v3 = np.asanyarray(tests.make_random_image(v3file, (5, 5, 5, 3)).dataobj)
l1 = np.asanyarray(tests.make_random_image(l1file, (5, 5, 5)) .dataobj)
l2 = np.asanyarray(tests.make_random_image(l2file, (5, 5, 5)) .dataobj)
l3 = np.asanyarray(tests.make_random_image(l3file, (5, 5, 5)) .dataobj)
dtiobj = dtifit.DTIFitTensor(testdir)
......
......@@ -27,7 +27,7 @@ def test_ensureIsImage():
for l in loaded:
assert isinstance(l, nib.nifti1.Nifti1Image)
assert np.all(img.get_data() == l.get_data())
assert np.all(np.asanyarray(img.dataobj) == np.asanyarray(l.dataobj))
l = None
loaded = None
......
......@@ -823,7 +823,7 @@ def _test_Image_changeData(imgtype):
img.register('name3', onDataRange, 'dataRange')
# Calculate the actual data range
data = img.nibImage.get_data()
data = np.asanyarray(img.nibImage.dataobj)
dmin = data.min()
dmax = data.max()
drange = dmax - dmin
......@@ -1098,14 +1098,16 @@ def _test_Image_save(imgtype):
# Load the image back in
img2 = fslimage.Image(img.dataSource)
assert img.saveState
assert img.dataSource == expDataSource
assert img2.saveState
assert img2.dataSource == expDataSource
if imgtype > 0:
assert np.all(np.isclose(img.voxToWorldMat, xform))
for i in (img, img2):
if imgtype > 0:
assert np.all(np.isclose(i.voxToWorldMat, xform))
for (x, y, z), v in zip(rvoxs, rvals):
assert np.isclose(img[x, y, z], v)
for (x, y, z), v in zip(rvoxs, rvals):
assert np.isclose(i[x, y, z], v)
img = None
img2 = None
......@@ -1135,7 +1137,7 @@ def _test_Image_init_xform(imgtype):
# an image created off a
# header should have
# identical sform/qform
fimg = fslimage.Image(img.get_data(), header=img.header)
fimg = fslimage.Image(np.asanyarray(img.dataobj), header=img.header)
fsform, fsform_code = fimg.header.get_sform(True)
fqform, fqform_code = fimg.header.get_qform(True)
......@@ -1153,7 +1155,7 @@ def _test_Image_init_xform(imgtype):
# set to that xform,
# qform to None, and
# and codes set to (s2, q0)
fimg = fslimage.Image(img.get_data(), xform=sform)
fimg = fslimage.Image(np.asanyarray(img.dataobj), xform=sform)
fsform, fsform_code = fimg.header.get_sform(True)
fqform, fqform_code = fimg.header.get_qform(True)
......@@ -1174,7 +1176,7 @@ def _test_Image_init_xform(imgtype):
rxform = affine.compose(np.random.random(3),
np.random.random(3),
np.random.random(3))
fimg = fslimage.Image(img.get_data(),
fimg = fslimage.Image(np.asanyarray(img.dataobj),
header=img.header,
xform=rxform)
......
......@@ -790,7 +790,7 @@ def _test_ImageWrapper_write_out(niters, seed, threaded):
wrapper[tuple(sliceobjs)] = newData
_ImageWraper_busy_wait(wrapper)
expLo, expHi = coverageDataRange(img.get_data(), cov, slices)
expLo, expHi = coverageDataRange(np.asanyarray(img.dataobj), cov, slices)
newLo, newHi = wrapper.dataRange
# print('Old range: {} - {}'.format(clo, chi))
......@@ -917,10 +917,10 @@ def _test_ImageWrapper_write_in_overlap(niters, seed, threaded):
print('Expected range: {} - {}'.format(expLo, expHi))
print('New range: {} - {}'.format(newLo, newHi))
print('Slice min/max: {} - {}'.format(img.get_data()[tuple(sliceobjs)].min(),
img.get_data()[tuple(sliceobjs)].max()))
print('Data min/max: {} - {}'.format(img.get_data().min(),
img.get_data().max()))
print('Slice min/max: {} - {}'.format(np.asanyarray(img.dataobj)[tuple(sliceobjs)].min(),
np.asanyarray(img.dataobj)[tuple(sliceobjs)].max()))
print('Data min/max: {} - {}'.format(np.asanyarray(img.dataobj).min(),
np.asanyarray(img.dataobj).max()))
assert np.all(newCov == expCov)
......
......@@ -14,6 +14,7 @@ import os
import shutil
import tempfile
import numpy as np
import nibabel as nib
import fsl.utils.imcp as imcp
......@@ -31,7 +32,7 @@ def print(*args, **kwargs):
def makeImage(filename):
return hash(make_random_image(filename).get_data().tobytes())
return hash(np.asanyarray(make_random_image(filename).dataobj).tobytes())
def checkImageHash(filename, datahash):
......@@ -39,7 +40,7 @@ def checkImageHash(filename, datahash):
"""
img = nib.load(filename)
assert hash(img.get_data().tobytes()) == datahash
assert hash(np.asanyarray(img.dataobj).tobytes()) == datahash
def checkFilesToExpect(files, outdir, outputType, datahashes):
......
......@@ -31,7 +31,7 @@ def test_MGHImage():
v2s = nbimg.header.get_vox2ras_tkr()
w2s = affine.concat(v2s, affine.invert(nbimg.affine))
assert np.all(np.isclose(img[:], nbimg.get_data()))
assert np.all(np.isclose(img[:], np.asanyarray(nbimg.dataobj)))
assert np.all(np.isclose(img.voxToWorldMat, nbimg.affine))
assert np.all(np.isclose(img.voxToSurfMat, v2s))
assert np.all(np.isclose(img.surfToVoxMat, affine.invert(v2s)))
......@@ -44,7 +44,7 @@ def test_MGHImage():
# Load from an in-memory nibabel object
img = fslmgh.MGHImage(nbimg)
assert np.all(np.isclose(img[:], nbimg.get_data()))
assert np.all(np.isclose(img[:], np.asanyarray(nbimg.dataobj)))
assert np.all(np.isclose(img.voxToWorldMat, nbimg.affine))
assert img.dataSource is None
assert img.mghImageFile is None
......
......@@ -44,7 +44,7 @@ def test_help():
def test_linear():
def test_linear(seed):
with tempdir.tempdir():
src2ref = _random_affine()
......@@ -65,7 +65,7 @@ def test_linear():
assert np.all(np.isclose(result.data, expect))
def test_nonlinear():
def test_nonlinear(seed):
with tempdir.tempdir():
src2ref = _random_affine()
......@@ -92,7 +92,7 @@ def test_nonlinear():
assert np.all(np.isclose(result, expect))
def test_linear_altref():
def test_linear_altref(seed):
with tempdir.tempdir():
src2ref = affine.scaleOffsetXform([1, 1, 1], [5, 5, 5])
......@@ -119,7 +119,7 @@ def test_linear_altref():
assert np.all(result.data == expect)
def test_nonlinear_altref():
def test_nonlinear_altref(seed):
with tempdir.tempdir():
src2ref = affine.scaleOffsetXform([1, 1, 1], [5, 5, 5])
......@@ -149,7 +149,7 @@ def test_nonlinear_altref():
assert np.all(result.data == expect)
def test_linear_altsrc():
def test_linear_altsrc(seed):
with tempdir.tempdir():
src2ref = _random_affine()
......@@ -198,7 +198,7 @@ def test_linear_altsrc():
assert np.all(np.isclose(outoff.data, expoff))
def test_nonlinear_altsrc():
def test_nonlinear_altsrc(seed):
with tempdir.tempdir():
src2ref = _random_affine()
......@@ -252,13 +252,20 @@ def test_nonlinear_altsrc():
assert outhi .sameSpace(ref)
assert outoff.sameSpace(ref)
# We get boundary issues just at the first
# voxel, so I'm masking that voxel out
for img in (out, outlo, outhi, outoff,
exp, explo, exphi, expoff):
img[0, 0, 0] = 0
assert np.all(np.isclose(out .data, exp))
assert np.all(np.isclose(outlo .data, explo))
assert np.all(np.isclose(outhi .data, exphi))
assert np.all(np.isclose(outoff.data, expoff))
# We get boundary cropping,
# so ignore edge slices
out = out .data[1:-1, 1:-1, 1:-1]
outlo = outlo .data[1:-1, 1:-1, 1:-1]
outhi = outhi .data[1:-1, 1:-1, 1:-1]
outoff = outoff.data[1:-1, 1:-1, 1:-1]
exp = exp[ 1:-1, 1:-1, 1:-1]
explo = explo[ 1:-1, 1:-1, 1:-1]
exphi = exphi[ 1:-1, 1:-1, 1:-1]
expoff = expoff[ 1:-1, 1:-1, 1:-1]
tol = dict(atol=1e-3, rtol=1e-3)
assert np.all(np.isclose(out, exp, **tol))
assert np.all(np.isclose(outlo, explo, **tol))
assert np.all(np.isclose(outhi, exphi, **tol))
assert np.all(np.isclose(outoff, expoff, **tol))
......@@ -11,6 +11,7 @@ import tempfile
import pytest
import numpy as np
import nibabel as nib
from fsl.utils.tempdir import tempdir
......@@ -269,7 +270,7 @@ def test_imcp_script_shouldPass(move=False):
for inf in infiles:
img = nib.load(op.join(tindir, inf),
mmap=False)
imghash = hash(img.get_data().tobytes())
imghash = hash(np.asanyarray(img.dataobj).tobytes())
img = None
imageHashes.append(imghash)
......
......@@ -231,8 +231,8 @@ def test_fileOrImage():
@wutils.fileOrImage('img1', 'img2', 'output')
def func(img1, **kwargs):
img1 = nib.load(img1).get_data()
img2 = nib.load(kwargs['img2']).get_data()
img1 = np.asanyarray(nib.load(img1).dataobj)
img2 = np.asanyarray(nib.load(kwargs['img2']).dataobj)
output = nib.nifti1.Nifti1Image(img1 * img2, np.eye(4))
nib.save(output, kwargs['output'])
......@@ -247,43 +247,43 @@ def test_fileOrImage():
# file file file
func('img1.nii', img2='img2.nii', output='output.nii')
assert np.all(nib.load('output.nii').get_data() == expected)
assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
os.remove('output.nii')
# file file array
result = func('img1.nii', img2='img2.nii', output=wutils.LOAD)['output']
assert np.all(result.get_data() == expected)
assert np.all(np.asanyarray(result.dataobj) == expected)
# file array file
func('img1.nii', img2=img2, output='output.nii')
assert np.all(nib.load('output.nii').get_data() == expected)
assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
os.remove('output.nii')
# file array array
result = func('img1.nii', img2=img2, output=wutils.LOAD)['output']
assert np.all(result.get_data() == expected)
assert np.all(np.asanyarray(result.dataobj) == expected)
# array file file
func(img1, img2='img2.nii', output='output.nii')
assert np.all(nib.load('output.nii').get_data() == expected)
assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
os.remove('output.nii')
# array file array
result = func(img1, img2='img2.nii', output=wutils.LOAD)['output']
assert np.all(result.get_data() == expected)
assert np.all(np.asanyarray(result.dataobj) == expected)
# array array file
func(img1, img2=img2, output='output.nii')
assert np.all(nib.load('output.nii').get_data() == expected)
assert np.all(np.asanyarray(nib.load('output.nii').dataobj) == expected)
os.remove('output.nii')
# array array array
result = func(img1, img2=img2, output=wutils.LOAD)['output']