def load_ROI_mask(self):
proxy = nib.load(self.FLAIR_FILE)
image_array = np.asarray(proxy.dataobj)
mask = np.ones_like(image_array)
mask[np.where(image_array < 90)] = 0
# img = nib.Nifti1Image(mask, proxy.affine)
# nib.save(img, join(modalities_path,'mask.nii.gz'))
struct_element_size = (20, 20, 20)
mask_augmented = np.pad(mask, [(21, 21), (21, 21), (21, 21)], 'constant', constant_values=(0, 0))
mask_augmented = binary_closing(mask_augmented, structure=np.ones(struct_element_size, dtype=bool)).astype(
np.int)
return mask_augmented[21:-21, 21:-21, 21:-21].astype('bool')
python类Nifti1Image()的实例源码
def export_data(prediction_dir, nii_image_dir, tfrecords_dir, export_dir, transformation=None):
for file_path in os.listdir(prediction_dir):
name, prediction, probability = read_prediction_file(os.path.join(prediction_dir, file_path))
if transformation:
image, ground_truth = get_original_image(os.path.join(tfrecords_dir, name + '.tfrecord'), False)
prediction = transformation.transform_image(prediction, probability, image)
# build cv_predictions .nii image
img = nib.Nifti1Image(prediction, np.eye(4))
img.set_data_dtype(dtype=np.uint8)
path = os.path.join(nii_image_dir, name)
adc_name = next(l for l in os.listdir(path) if 'MR_ADC' in l)
export_image = nib.load(os.path.join(nii_image_dir, name, adc_name, adc_name + '.nii'))
i = export_image.get_data()
i[:] = img.get_data()
# set name to specification and export
_id = next(l for l in os.listdir(path) if 'MR_MTT' in l).split('.')[-1]
export_path = os.path.join(export_dir, 'SMIR.' + name + '.' + _id + '.nii')
nib.save(export_image, os.path.join(export_path))
def batch_works(k):
if k == n_processes - 1:
paths = all_paths[k * int(len(all_paths) / n_processes) : ]
else:
paths = all_paths[k * int(len(all_paths) / n_processes) : (k + 1) * int(len(all_paths) / n_processes)]
for path in paths:
probs = np.load(os.path.join(input_path, path))
pred = np.argmax(probs, axis=3)
fg_prob = 1 - probs[..., 0]
pred = clean_contour(fg_prob, pred)
seg = np.zeros(pred.shape, dtype=np.int16)
seg[pred == 1] = 1
seg[pred == 2] = 2
seg[pred == 3] = 4
img = nib.Nifti1Image(seg, np.eye(4))
nib.save(img, os.path.join(output_path, path.replace('_probs.npy', '.nii.gz')))
def _make_test_data(n_subjects=8, noisy=False):
rng = np.random.RandomState(0)
shape = (20, 20, 1)
components = _make_components(shape)
if noisy: # Creating noisy non positive data
components[rng.randn(*components.shape) > .8] *= -5.
for component in components:
assert(component.max() <= -component.min()) # Goal met ?
# Create a "multi-subject" dataset
data = _make_data_from_components(components, n_subjects=n_subjects)
affine = np.eye(4)
mask_img = nibabel.Nifti1Image(np.ones(shape, dtype=np.int8),
affine)
masker = MultiNiftiMasker(mask_img).fit()
init = components + 1 * rng.randn(*components.shape)
components = masker.inverse_transform(components)
init = masker.inverse_transform(init)
data = masker.inverse_transform(data)
return data, mask_img, components, init
def _binarize(in_file, threshold=0.5, out_file=None):
import os.path as op
import numpy as np
import nibabel as nb
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_bin{}'.format(fname, ext))
nii = nb.load(in_file)
data = nii.get_data()
data[data <= threshold] = 0
data[data > 0] = 1
hdr = nii.header.copy()
hdr.set_data_dtype(np.uint8)
nb.Nifti1Image(data.astype(np.uint8), nii.affine, hdr).to_filename(
out_file)
return out_file
def image_gradient(in_file, snr, out_file=None):
"""Computes the magnitude gradient of an image using numpy"""
import os.path as op
import numpy as np
import nibabel as nb
from scipy.ndimage import gaussian_gradient_magnitude as gradient
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_grad{}'.format(fname, ext))
imnii = nb.load(in_file)
data = imnii.get_data().astype(np.float32) # pylint: disable=no-member
datamax = np.percentile(data.reshape(-1), 99.5)
data *= 100 / datamax
grad = gradient(data, 3.0)
gradmax = np.percentile(grad.reshape(-1), 99.5)
grad *= 100.
grad /= gradmax
nb.Nifti1Image(grad, imnii.get_affine(), imnii.get_header()).to_filename(out_file)
return out_file
def thresh_image(in_file, thres=0.5, out_file=None):
"""Thresholds an image"""
import os.path as op
import nibabel as nb
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_thresh{}'.format(fname, ext))
im = nb.load(in_file)
data = im.get_data()
data[data < thres] = 0
data[data > 0] = 1
nb.Nifti1Image(
data, im.get_affine(), im.get_header()).to_filename(out_file)
return out_file
def test_set_new_data_simple_modifications():
aff = np.eye(4); aff[2, 1] = 42.0
im_0 = nib.Nifti1Image(np.zeros([3,3,3]), affine=aff)
im_0_header = im_0.header
# default intent_code
assert_equals(im_0_header['intent_code'], 0)
# change intento code
im_0_header['intent_code'] = 5
# generate new nib from the old with new data
im_1 = set_new_data(im_0, np.ones([3,3,3]))
im_1_header = im_1.header
# see if the infos are the same as in the modified header
assert_array_equal(im_1.get_data()[:], np.ones([3,3,3]))
assert_equals(im_1_header['intent_code'], 5)
assert_array_equal(im_1.get_affine(), aff)
def test_stack_images_cascade():
d = 2
im1 = nib.Nifti1Image(np.zeros([d,d]), affine=np.eye(4))
assert_array_equal(im1.shape, (d, d))
list_images1 = [im1] * d
im2 = stack_images(list_images1)
assert_array_equal(im2.shape, (d,d,d))
list_images2 = [im2] * d
im3 = stack_images(list_images2)
assert_array_equal(im3.shape, (d,d,d,d))
list_images3 = [im3] * d
im4 = stack_images(list_images3)
assert_array_equal(im4.shape, (d, d, d, d, d))
def test_remove_nan():
data_ts = np.array([[[0, 1, 2, 3],
[4, 5, np.nan, 7],
[8, 9, 10, np.nan]],
[[12, 13, 14, 15],
[16, np.nan, 18, 19],
[20, 21, 22, 23]]])
im = nib.Nifti1Image(data_ts, np.eye(4))
im_no_nan = remove_nan(im)
data_no = np.array([[[0, 1, 2, 3],
[4, 5, 0, 7],
[8, 9, 10, 0]],
[[12, 13, 14, 15],
[16, 0, 18, 19],
[20, 21, 22, 23]]])
assert_array_equal(im_no_nan.get_data(), data_no)
def tensor2fa(tensors, tensor_name, dti, derivdir, qcdir):
'''
outdir: location of output directory.
fname: name of output fa map file. default is none (name created based on
input file)
'''
dti_data = nb.load(dti)
affine = dti_data.get_affine()
dti_data = dti_data.get_data()
# create FA map
FA = fractional_anisotropy(tensors.evals)
FA[np.isnan(FA)] = 0
# generate the RGB FA map
FA = np.clip(FA, 0, 1)
RGB = color_fa(FA, tensors.evecs)
fname = os.path.split(tensor_name)[1].split(".")[0] + '_fa_rgb.nii.gz'
fa = nb.Nifti1Image(np.array(255 * RGB, 'uint8'), affine)
nb.save(fa, derivdir + fname)
fa_pngs(fa, fname, qcdir)
def test_ROIsPlot(oasis_dir):
""" the BET report capable test """
import nibabel as nb
import numpy as np
labels = nb.load(os.path.join(oasis_dir, 'T_template0_glm_4labelsJointFusion.nii.gz'))
data = labels.get_data()
out_files = []
ldata = np.zeros_like(data)
for i, l in enumerate([1, 3, 4, 2]):
ldata[data == l] = 1
out_files.append(os.path.abspath('label%d.nii.gz' % i))
lfile = nb.Nifti1Image(ldata, labels.affine, labels.header)
lfile.to_filename(out_files[-1])
roi_rpt = ROIsPlot(
generate_report=True,
in_file=os.path.join(oasis_dir, 'T_template0.nii.gz'),
in_mask=out_files[-1],
in_rois=out_files[:-1],
colors=['g', 'y']
)
_smoke_test_report(roi_rpt, 'testROIsPlot.svg')
def test_register_dwi():
fdata, fbval, fbvec = dpd.get_data('small_64D')
with nbtmp.InTemporaryDirectory() as tmpdir:
# Use an abbreviated data-set:
img = nib.load(fdata)
data = img.get_data()[..., :10]
nib.save(nib.Nifti1Image(data, img.affine),
op.join(tmpdir, 'data.nii.gz'))
# Convert from npy to txt:
bvals = np.load(fbval)
bvecs = np.load(fbvec)
np.savetxt(op.join(tmpdir, 'bvals.txt'), bvals[:10])
np.savetxt(op.join(tmpdir, 'bvecs.txt'), bvecs[:10])
reg_file = register_dwi(op.join(tmpdir, 'data.nii.gz'),
op.join(tmpdir, 'bvals.txt'),
op.join(tmpdir, 'bvecs.txt'))
npt.assert_(op.exists(reg_file))
def test_set_new_data_basic():
data_1 = np.random.normal(5, 10, [10, 10, 5])
affine_1 = np.diag([1, 2, 3, 1])
data_2 = np.random.normal(5, 10, [3, 2, 4]).astype(np.float32)
im_data_1 = nib.Nifti1Image(data_1, affine_1)
im_data_1.set_data_dtype(np.uint8)
im_data_1.header['descrip'] = 'Spam'
im_data_2 = set_new_data(im_data_1, data_2)
assert_array_equal(im_data_2.get_data(), data_2)
assert_array_equal(im_data_2.get_affine(), affine_1)
assert_equal(im_data_2.header['descrip'], b'Spam')
assert_equal(im_data_1.get_data_dtype(), np.uint8)
assert_equal(im_data_2.get_data_dtype(), np.float32)
def preprocess(inputfile, outputfile, order=0, df=None, input_key=None, output_key=None):
img = nib.load(inputfile)
data = img.get_data()
affine = img.affine
zoom = img.header.get_zooms()[:3]
data, affine = reslice(data, affine, zoom, (1., 1., 1.), order)
data = np.squeeze(data)
data = np.pad(data, [(0, 256 - len_) for len_ in data.shape], "constant")
if order == 0:
if df is not None:
tmp = np.zeros_like(data)
for target, source in zip(df[output_key], df[input_key]):
tmp[np.where(data == source)] = target
data = tmp
data = np.int32(data)
assert data.ndim == 3, data.ndim
else:
data_sub = data - gaussian_filter(data, sigma=1)
img = sitk.GetImageFromArray(np.copy(data_sub))
img = sitk.AdaptiveHistogramEqualization(img)
data_clahe = sitk.GetArrayFromImage(img)[:, :, :, None]
data = np.concatenate((data_clahe, data[:, :, :, None]), 3)
data = (data - np.mean(data, (0, 1, 2))) / np.std(data, (0, 1, 2))
assert data.ndim == 4, data.ndim
assert np.allclose(np.mean(data, (0, 1, 2)), 0.), np.mean(data, (0, 1, 2))
assert np.allclose(np.std(data, (0, 1, 2)), 1.), np.std(data, (0, 1, 2))
data = np.float32(data)
img = nib.Nifti1Image(data, affine)
nib.save(img, outputfile)
def roi_from_bbox(
input_file,
bbox,
output_file):
""" Create a ROI image from a bounding box.
Parameters
----------
input_file: str
the reference image where the bbox is defined.
bbox: 6-uplet
the corner of the bbox in voxel coordinates: xmin, xmax, ymin, ymax,
zmin, zmax.
output_file: str
the desired ROI image file.
"""
# Load the reference image and generate a grid
im = nibabel.load(input_file)
xv, yv, zv = numpy.meshgrid(
numpy.linspace(0, im.shape[0] - 1, im.shape[0]),
numpy.linspace(0, im.shape[1] - 1, im.shape[1]),
numpy.linspace(0, im.shape[2] - 1, im.shape[2]))
xv = xv.astype(int)
yv = yv.astype(int)
zv = zv.astype(int)
# Intersect the grid with the bbox
xa = numpy.bitwise_and(xv >= bbox[0], xv <= bbox[1])
ya = numpy.bitwise_and(yv >= bbox[2], yv <= bbox[3])
za = numpy.bitwise_and(zv >= bbox[4], zv <= bbox[5])
# Generate bbox indices
indices = numpy.bitwise_and(numpy.bitwise_and(xa, ya), za)
# Generate/save ROI
roi = numpy.zeros(im.shape, dtype=int)
roi[xv[indices].tolist(), yv[indices].tolist(), zv[indices].tolist()] = 1
roi_im = nibabel.Nifti1Image(roi, affine=im.get_affine())
nibabel.save(roi_im, output_file)
def extract_image(in_file, index, out_file=None):
""" Extract the image at 'index' position.
Parameters
----------
in_file: str (mandatory)
the input image.
index: int (mandatory)
the index of last image dimention to extract.
out_file: str (optional, default None)
the name of the extracted image file.
Returns
-------
out_file: str
the name of the extracted image file.
"""
# Set default output if necessary
dirname = os.path.dirname(in_file)
basename = os.path.basename(in_file).split(".")[0]
if out_file is None:
out_file = os.path.join(
dirname, "extract{0}_{1}.nii.gz".format(index, basename))
# Extract the image of interest
image = nibabel.load(in_file)
affine = image.get_affine()
extracted_array = image.get_data()[..., index]
extracted_image = nibabel.Nifti1Image(extracted_array, affine)
nibabel.save(extracted_image, out_file)
return out_file
def save_nii(img_path, data, affine, header):
'''
Shortcut to save a nifty file
'''
nimg = nib.Nifti1Image(data, affine=affine, header=header)
nimg.to_filename(img_path)
def saveImageAsNifti(imageToSave,
imageName,
imageOriginalName,
imageType):
printFileNames = False
if printFileNames == True:
print(" ... Saving image in {}".format(imageName))
[imageData,img_proxy] = load_nii(imageOriginalName, printFileNames)
# Generate the nii file
niiToSave = nib.Nifti1Image(imageToSave, img_proxy.affine)
niiToSave.set_data_dtype(imageType)
dim = len(imageToSave.shape)
zooms = list(img_proxy.header.get_zooms()[:dim])
if len(zooms) < dim :
zooms = zooms + [1.0]*(dim-len(zooms))
niiToSave.header.set_zooms(zooms)
nib.save(niiToSave, imageName)
print "... Image succesfully saved..."
## ========= For Matlab files ========== ##
def extract_data(nifti_file, mask_file, out_file, zscore, detrend, smoothing_fwmw):
if mask_file is None:
#whole brain, get coordinate info from nifti_file itself
mask = nib.load(nifti_file.name)
else:
mask = nib.load(mask_file.name)
affine = mask.get_affine()
if mask_file is None:
mask_data = mask.get_data()
if mask_data.ndim == 4:
#get mask in 3D
img_data_type = mask.header.get_data_dtype()
n_tr = mask_data.shape[3]
mask_data = mask_data[:,:,:,n_tr//2].astype(bool)
mask = nib.Nifti1Image(mask_data.astype(img_data_type), affine)
else:
mask_data = mask_data.astype(bool)
else:
mask_data = mask.get_data().astype(bool)
#get voxel coordinates
R = np.float64(np.argwhere(mask_data))
#get scanner RAS coordinates based on voxel coordinates
if affine is not []:
R = (np.dot(affine[:3,:3], R.T) + affine[:3,3:4]).T
#get ROI data, and run preprocessing
nifti_masker = NiftiMasker(mask_img=mask, standardize=zscore, detrend=detrend, smoothing_fwhm=smoothing_fwmw)
img = nib.load(nifti_file.name)
all_images = np.float64(nifti_masker.fit_transform(img))
data = all_images.T.copy()
#save data
subj_data = {'data': data, 'R': R}
scipy.io.savemat(out_file.name, subj_data)
def save(self, image, outpath):
""" A method that save the image data and associated metadata.
Parameters
----------
image: Image
the image to be saved.
outpath: str
the path where the the image will be saved.
"""
diag = (1. / image.spacing).tolist() + [1]
_image = nibabel.Nifti1Image(image.data, numpy.diag(diag))
nibabel.save(_image, outpath)
def write_nifti(data, output_fname, header=None, affine=None, use_data_dtype=True, **kwargs):
"""Write data to a nifti file.
This will write the output directory if it does not exist yet.
Args:
data (ndarray): the data to write to that nifti file
output_fname (str): the name of the resulting nifti file, this function will append .nii.gz if no
suitable extension is given.
header (nibabel header): the nibabel header to use as header for the nifti file. If None we will use
a default header.
affine (ndarray): the affine transformation matrix
use_data_dtype (boolean): if we want to use the dtype from the data instead of that from the header
when saving the nifti.
**kwargs: other arguments to Nifti2Image from NiBabel
"""
if header is None:
header = nib.nifti2.Nifti2Header()
if use_data_dtype:
header = copy.deepcopy(header)
dtype = data.dtype
if data.dtype == np.bool:
dtype = np.char
try:
header.set_data_dtype(dtype)
except nib.spatialimages.HeaderDataError:
pass
if not (output_fname.endswith('.nii.gz') or output_fname.endswith('.nii')):
output_fname += '.nii.gz'
if not os.path.exists(os.path.dirname(output_fname)):
os.makedirs(os.path.dirname(output_fname))
if isinstance(header, nib.nifti2.Nifti2Header):
format = nib.Nifti2Image
else:
format = nib.Nifti1Image
format(data, affine, header=header, **kwargs).to_filename(output_fname)
def remove_minor_cc(vol_data, rej_ratio, rename_map):
"""Remove small connected components refer to rejection ratio"""
"""Usage
# rename_map = [0, 205, 420, 500, 550, 600, 820, 850]
# nii_path = '/home/xinyang/project_xy/mmwhs2017/dataset/ct_output/test/test_4.nii'
# vol_file = nib.load(nii_path)
# vol_data = vol_file.get_data().copy()
# ref_affine = vol_file.affine
# rem_vol = remove_minor_cc(vol_data, rej_ratio=0.2, class_n=8, rename_map=rename_map)
# # save
# rem_path = 'rem_cc.nii'
# rem_vol_file = nib.Nifti1Image(rem_vol, ref_affine)
# nib.save(rem_vol_file, rem_path)
#===# possible be parallel in future
"""
rem_vol = copy.deepcopy(vol_data)
class_n = len(rename_map)
# retrieve all classes
for c in range(1, class_n):
print 'processing class %d...' % c
class_idx = (vol_data==rename_map[c])*1
class_vol = np.sum(class_idx)
labeled_cc, num_cc = measurements.label(class_idx)
# retrieve all connected components in this class
for cc in range(1, num_cc+1):
single_cc = ((labeled_cc==cc)*1)
single_vol = np.sum(single_cc)
# remove if too small
if single_vol / (class_vol*1.0) < rej_ratio:
rem_vol[labeled_cc==cc] = 0
return rem_vol
def read_series(dicoms,return_nifti=True):
'''read_series will read in a series of dicoms belonging to a group
:param dicoms: a list of dicom files to parse, assumed in same series and
equal size
:param return_nifti: If True (default) will turn image as Nifti File
'''
# Sort the dicoms
dicoms.sort()
# Get the size of the image
params = sniff_header(dicoms[0])
xdim = params['xdim']
ydim = params['ydim']
windom_center = params['window_center']
bot.debug("First dicom found with dimension %s by %s, using as standard." %(xdim,ydim))
# Let's get ordering of images based on InstanceNumber
ordered = dict()
for d in range(len(dicoms)):
ds = dicom.read_file(dicoms[d])
if ds.Rows == xdim and ds.Columns == ydim:
ordered[int(ds.InstanceNumber)] = ds.pixel_array
# Sort by order
zdim = len(ordered)
data = numpy.ndarray((xdim,ydim,zdim))
# Start at window center, then go back to zero
index = 0
for key in sorted(ordered.keys()):
data[:,:,index] = ordered[key]
index +=1
if return_nifti == True:
affine = numpy.diag((1,1,1,0))
data = nibabel.Nifti1Image(data,affine=affine)
return data
def _enhance(in_file, out_file=None):
import os.path as op
import numpy as np
import nibabel as nb
if out_file is None:
fname, ext = op.splitext(op.basename(in_file))
if ext == '.gz':
fname, ext2 = op.splitext(fname)
ext = ext2 + ext
out_file = op.abspath('{}_enhanced{}'.format(fname, ext))
imnii = nb.load(in_file)
data = imnii.get_data().astype(np.float32) # pylint: disable=no-member
range_max = np.percentile(data[data > 0], 99.98)
range_min = np.median(data[data > 0])
# Resample signal excess pixels
excess = np.where(data > range_max)
data[excess] = 0
data[excess] = np.random.choice(data[data > range_min], size=len(excess[0]))
nb.Nifti1Image(data, imnii.get_affine(), imnii.get_header()).to_filename(
out_file)
return out_file
def reorient_and_discard_non_steady(in_file, float32=False):
import nibabel as nb
import os
import numpy as np
import nibabel as nb
from statsmodels.robust.scale import mad
_, outfile = os.path.split(in_file)
nii = nb.as_closest_canonical(nb.load(in_file))
in_data = nii.get_data()
# downcast to reduce space consumption and improve performance
if float32 and np.dtype(in_data.dtype).itemsize > 4:
in_data = in_data.astype(np.float32)
data = in_data[:, :, :, :50]
timeseries = data.max(axis=0).max(axis=0).max(axis=0)
outlier_timecourse = (timeseries - np.median(timeseries)) / mad(
timeseries)
exclude_index = 0
for i in range(10):
if outlier_timecourse[i] > 10:
exclude_index += 1
else:
break
nb.Nifti1Image(in_data[:, :, :, exclude_index:], nii.affine, nii.header).to_filename(outfile)
nii.uncache()
return exclude_index, os.path.abspath(outfile)
def test_compare_two_nib_equals():
im_0 = nib.Nifti1Image(np.zeros([3, 3, 3]), affine=np.eye(4))
im_1 = nib.Nifti1Image(np.zeros([3, 3, 3]), affine=np.eye(4))
assert_equals(compare_two_nib(im_0, im_1), True)
def test_compare_two_nib_different_nifti_version():
im_0 = nib.Nifti1Image(np.zeros([3, 3, 3]), affine=np.eye(4))
im_1 = nib.Nifti2Image(np.zeros([3, 3, 3]), affine=np.eye(4))
assert_equals(compare_two_nib(im_0, im_1), False)
def test_compare_two_nib_different_affine():
aff_1 = np.eye(4)
aff_1[3,3] = 5
im_0 = nib.Nifti1Image(np.zeros([3, 3, 3]), affine=np.eye(4))
im_1 = nib.Nifti1Image(np.zeros([3, 3, 3]), affine=aff_1)
assert_equals(compare_two_nib(im_0, im_1), False)
# def test_eliminates_consecutive_duplicates():
#
# l_in = [0,0,0,1,1,2,3,4,5,5,5,6,7,8,9]
# l_out = range(10)
# assert_array_equal(eliminates_consecutive_duplicates(l_in), l_out)
def test_covariance_matrices():
arr_1 = np.array([[[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 2, 2, 0],
[0, 1, 1, 1, 1, 1, 2, 2, 0],
[0, 1, 1, 1, 1, 1, 0, 3, 3],
[0, 0, 0, 0, 0, 0, 0, 3, 3],
[0, 0, 0, 0, 0, 3, 3, 3, 3]],
[[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 0, 0, 0],
[0, 1, 1, 1, 1, 1, 2, 2, 0],
[0, 1, 1, 1, 1, 1, 2, 2, 0],
[0, 1, 1, 1, 1, 1, 0, 0, 3],
[0, 0, 0, 0, 0, 0, 0, 3, 3],
[0, 0, 0, 0, 0, 0, 3, 3, 3]]
])
im1 = nib.Nifti1Image(arr_1, np.eye(4))
cov = covariance_matrices(im1, [1, 2, 3])
assert len(cov) == 3
for i in cov:
assert_array_equal(i.shape, [3, 3])
if np.count_nonzero(i - np.diag(np.diagonal(i))) == 0:
assert_array_equal(np.diag(np.diag(i)), i)
cov1 = covariance_matrices(im1, [1, 2, 3, 4])
assert_array_equal(cov1[-1], np.nan * np.ones([3, 3]))
