"""Read/Write Nifti-1 files
"""
# Copyright (c) 2013-2018 Erling Andersen, Haukeland University Hospital, Bergen, Norway
import os.path
import tempfile
import logging
import numpy as np
import imagedata.formats
import imagedata.axis
from imagedata.formats.abstractplugin import AbstractPlugin
import nibabel
import nibabel.spatialimages
# noinspection PyUnresolvedReferences
[docs]class NiftiPlugin(AbstractPlugin):
"""Read/write Nifti-1 files."""
name = "nifti"
description = "Read and write Nifti-1 files."
authors = "Erling Andersen"
version = "1.0.0"
url = "www.helse-bergen.no"
"""
data - getter and setter - NumPy array
read() method
write() method
"""
def __init__(self):
super(NiftiPlugin, self).__init__(self.name, self.description,
self.authors, self.version, self.url)
self.shape = None
self.slices = None
self.spacing = None
self.transformationMatrix = None
self.imagePositions = None
self.tags = None
self.origin = None
self.orientation = None
self.normal = None
self.output_sort = None
def _read_image(self, f, opts, hdr):
"""Read image data from given file handle
Args:
self: format plugin instance
f: file handle or filename (depending on self._need_local_file)
opts: Input options (dict)
hdr: Header dict
Returns:
Tuple of
hdr: Header dict
Return values:
- info: Internal data for the plugin
None if the given file should not be included (e.g. raw file)
si: numpy array (multi-dimensional)
"""
logging.debug("niftiplugin::read filehandle {}".format(f))
# TODO: Read nifti directly from open file object
# Should be able to do something like:
#
# with archive.open(member_name) as member:
# # Create a nibabel image using
# # the existing file handle.
# fmap = nibabel.nifti1.Nifti1Image.make_file_map()
# #nibabel.nifti1.Nifti1Header
# fmap['image'].fileobj = member
# img = nibabel.Nifti1Image.from_file_map(fmap)
#
logging.debug("niftiplugin::read load f {}".format(f))
try:
img = nibabel.load(f)
except nibabel.spatialimages.ImageFileError:
raise imagedata.formats.NotImageError(
'{} does not look like a nifti file.'.format(f))
except Exception:
raise
info = img.header
si = self._reorder_to_dicom(img.get_data(), flip=False, flipud=True)
return info, si
def _need_local_file(self):
"""Do the plugin need access to local files?
Returns:
Boolean:
- True: The plugin need access to local filenames
- False: The plugin can access files given by an open file handle
"""
return True
def _set_tags(self, image_list, hdr, si):
"""Set header tags.
Args:
self: format plugin instance
image_list: list with (info,img) tuples
hdr: Header dict
si: numpy array (multi-dimensional)
Returns:
hdr: Header dict
"""
info, si = image_list[0]
_data_shape = info.get_data_shape()
nt = nz = 1
nx, ny = _data_shape[:2]
if len(_data_shape) > 2:
nz = _data_shape[2]
if len(_data_shape) > 3:
nt = _data_shape[3]
logging.debug("_set_tags: ny {}, nx {}, nz {}, nt {}".format(ny, nx, nz, nt))
logging.debug('NiftiPlugin.read: get_qform\n{}'.format(info.get_qform()))
logging.debug('NiftiPlugin.read: info.get_zooms() {}'.format(info.get_zooms()))
_xyzt_units = info.get_xyzt_units()
_data_zooms = info.get_zooms()
_dim_info = info.get_dim_info()
logging.debug("_set_tags: get_dim_info(): {}".format(info.get_dim_info()))
logging.debug("_set_tags: get_xyzt_units(): {}".format(info.get_xyzt_units()))
dt = dz = 1
dx, dy = _data_zooms[:2]
if len(_data_zooms) > 2:
dz = _data_zooms[2]
if len(_data_zooms) > 3:
dt = _data_zooms[3]
if _xyzt_units[0] == 'meter':
dx, dy, dz = dx * 1000., dy * 1000., dz * 1000.
elif _xyzt_units[0] == 'micron':
dx, dy, dz = dx / 1000., dy / 1000., dz / 1000.
if _xyzt_units[1] == 'msec':
dt = dt / 1000.
elif _xyzt_units[1] == 'usec':
dt = dt / 1000000.
self.spacing = (float(dz), float(dy), float(dx))
hdr['spacing'] = (float(dz), float(dy), float(dx))
# Simplify shape
self._reduce_shape(si)
sform, scode = info.get_sform(coded=True)
qform, qcode = info.get_qform(coded=True)
qfac = info['pixdim'][0]
if qfac not in (-1, 1):
raise ValueError('qfac (pixdim[0]) should be 1 or -1')
# Image orientation and positions
hdr['imagePositions'] = {}
if sform is not None and scode != 0:
logging.debug("Method 3 - sform: orientation")
# Note: rz, ry, rx, cz, cy, cx
iop = np.array([
sform[2, 0] / dx,
- sform[1, 0] / dx, # NIfTI is RAS+, DICOM is LPS+
- sform[0, 0] / dx, # NIfTI is RAS+, DICOM is LPS+
sform[2, 1] / dy,
- sform[1, 1] / dy, # NIfTI is RAS+, DICOM is LPS+
- sform[0, 1] / dy # NIfTI is RAS+, DICOM is LPS+
# - sform[2,1] / dy,
# sform[1,1] / dy, # NIfTI is RAS+, DICOM is LPS+
# sform[0,1] / dy # NIfTI is RAS+, DICOM is LPS+
])
for _slice in range(nz):
_p = np.array([
- (sform[0, 2] * _slice + sform[0, 3]), # NIfTI is RAS+, DICOM is LPS+
- (sform[1, 2] * _slice + sform[1, 3]), # NIfTI is RAS+, DICOM is LPS+
(sform[2, 2] * _slice + sform[2, 3])
])
hdr['imagePositions'][_slice] = _p[::-1] # Reverse x,y,z
elif qform is not None and qcode != 0:
logging.debug("Method 2 - qform: orientation")
qoffset_x, qoffset_y, qoffset_z = qform[0:3, 3]
a, b, c, d = info.get_qform_quaternion()
rx = - (a * a + b * b - c * c - d * d)
ry = - (2 * b * c + 2 * a * d)
rz = (2 * b * d - 2 * a * c)
cx = - (2 * b * c - 2 * a * d)
cy = - (a * a + c * c - b * b - d * d)
cz = (2 * c * d + 2 * a * b)
# normal from quaternion derived once and saved for position calculation ... do not handle qfac here ... do it later
tx = - (2 * b * d + 2 * a * c) # NIfTI is RAS+, DICOM is LPS+
ty = - (2 * c * d - 2 * a * b) # NIfTI is RAS+, DICOM is LPS+
tz = (a * a + d * d - c * c - b * b)
iop = np.array([rz, ry, rx, cz, cy, cx])
for _slice in range(nz):
_p = np.array([
tx * qfac * dz * _slice - qoffset_x, # NIfTI is RAS+, DICOM is LPS+
ty * qfac * dz * _slice - qoffset_y, # NIfTI is RAS+, DICOM is LPS+
tz * qfac * dz * _slice + qoffset_z
])
hdr['imagePositions'][_slice] = _p[::-1] # Reverse x,y,z
else:
logging.debug("Method 1 - assume axial: orientation")
iop = np.array([0, 0, 1, 0, 1, 0])
for _slice in range(nz):
_p = np.array([
0, # NIfTI is RAS+, DICOM is LPS+
0, # NIfTI is RAS+, DICOM is LPS+
dz * _slice
])
hdr['imagePositions'][_slice] = _p[::-1] # Reverse x,y,z
hdr['orientation'] = iop
self.shape = si.shape
times = [0]
if nt > 1:
times = np.arange(0, nt * dt, dt)
assert len(times) == nt, "Wrong timeline calculated (times={}) (nt={})".format(len(times), nt)
logging.debug("_set_tags: times {}".format(times))
tags = {}
for z in range(nz):
tags[z] = np.array(times)
hdr['tags'] = tags
axes = list()
if si.ndim > 3:
axes.append(imagedata.axis.UniformLengthAxis(
imagedata.formats.input_order_to_dirname_str(hdr['input_order']),
0,
nt,
dt)
)
if si.ndim > 2:
axes.append(imagedata.axis.UniformLengthAxis(
'slice',
0,
nz,
dz)
)
axes.append(imagedata.axis.UniformLengthAxis(
'row',
0,
ny,
dy)
)
axes.append(imagedata.axis.UniformLengthAxis(
'column',
0,
nx,
dx)
)
hdr['axes'] = axes
hdr['photometricInterpretation'] = 'MONOCHROME2'
hdr['color'] = False
# def nifti_to_affine(self, affine, shape):
#
# if len(shape) != 4:
# raise ValueError("4D only (was: %dD)" % len(shape))
#
# q = affine.copy()
#
# logging.debug("q from nifti_to_affine():\n{}".format(q))
# # Swap row 0 (z) and 2 (x)
# q[[0, 2],:] = q[[2, 0],:]
# # Swap column 0 (z) and 2 (x)
# q[:,[0, 2]] = q[:,[2, 0]]
# logging.debug("q swap nifti_to_affine():\n{}".format(q))
#
# analyze_to_dicom = np.eye(4)
# analyze_to_dicom[0,3] = 1
# analyze_to_dicom[1,3] = 1
# analyze_to_dicom[2,3] = 1
# dicom_to_analyze = np.linalg.inv(analyze_to_dicom)
# q = np.dot(q,dicom_to_analyze)
# logging.debug("q after dicom_to_analyze:\n{}".format(q))
#
# analyze_to_dicom = np.eye(4)
# analyze_to_dicom[0,3] = -1
# analyze_to_dicom[1,1] = -1
# rows = shape[2]
# analyze_to_dicom[1,3] = rows
# analyze_to_dicom[2,3] = -1
# dicom_to_analyze = np.linalg.inv(analyze_to_dicom)
# q = np.dot(q,dicom_to_analyze)
# logging.debug("q after rows dicom_to_analyze:\n{}".format(q))
#
# patient_to_tal = np.eye(4)
# patient_to_tal[0,0] = -1
# patient_to_tal[1,1] = -1
# tal_to_patient = np.linalg.inv(patient_to_tal)
# q = np.dot(tal_to_patient,q)
# logging.debug("q after tal_to_patient:\n{}".format(q))
#
# return q
# def affine_to_nifti(self, shape):
# q = self.transformationMatrix.copy()
# logging.debug("Affine from self.transformationMatrix:\n{}".format(q))
# # Swap row 0 (z) and 2 (x)
# q[[0, 2],:] = q[[2, 0],:]
# # Swap column 0 (z) and 2 (x)
# q[:,[0, 2]] = q[:,[2, 0]]
# logging.debug("Affine swap self.transformationMatrix:\n{}".format(q))
#
# # q now equals dicom_to_patient in spm_dicom_convert
#
# # Convert space
# analyze_to_dicom = np.eye(4)
# analyze_to_dicom[0,3] = -1
# analyze_to_dicom[1,1] = -1
# #if len(shape) == 3:
# # rows = shape[1]
# #else:
# # rows = shape[2]
# rows = shape[-2]
# analyze_to_dicom[1,3] = rows
# analyze_to_dicom[2,3] = -1
# logging.debug("analyze_to_dicom:\n{}".format(analyze_to_dicom))
#
# patient_to_tal = np.eye(4)
# patient_to_tal[0,0] = -1
# patient_to_tal[1,1] = -1
# logging.debug("patient_to_tal:\n{}".format(patient_to_tal))
#
# q = np.dot(patient_to_tal,q)
# logging.debug("q with patient_to_tal:\n{}".format(q))
# q = np.dot(q,analyze_to_dicom)
# # q now equals mat in spm_dicom_convert
#
# analyze_to_dicom = np.eye(4)
# analyze_to_dicom[0,3] = 1
# analyze_to_dicom[1,3] = 1
# analyze_to_dicom[2,3] = 1
# logging.debug("analyze_to_dicom:\n{}".format(analyze_to_dicom))
# q = np.dot(q,analyze_to_dicom)
#
# logging.debug("q nifti:\n{}".format(q))
# return q
@staticmethod
def _get_geometry_from_affine(hdr, q):
"""Extract geometry attributes from Nifti header
Args:
self: NiftiPlugin instance
q: nifti Qform
hdr['spacing']
Returns:
hdr: header dict
- hdr['imagePositions'][0]
- hdr['orientation']
- hdr['transformationMatrix']
"""
# Swap back from nifti patient space, flip x and y directions
affine = np.dot(np.diag([-1, -1, 1, 1]), q)
# Set imagePositions for first slice
x, y, z = affine[0:3, 3]
hdr['imagePositions'] = {0: np.array([z, y, x])}
logging.debug("getGeometryFromAffine: hdr imagePositions={}".format(hdr['imagePositions']))
# Set slice orientation
ds, dr, dc = hdr['spacing']
logging.debug("getGeometryFromAffine: spacing ds {}, dr {}, dc {}".format(ds, dr, dc))
colr = affine[:3, 0][::-1] / dr
colc = affine[:3, 1][::-1] / dc
# T0 = affine[:3,3][::-1]
orient = []
logging.debug("getGeometryFromAffine: affine\n{}".format(affine))
for i in range(3):
orient.append(colc[i])
for i in range(3):
orient.append(colr[i])
logging.debug("getGeometryFromAffine: orient {}".format(orient))
hdr['orientation'] = orient
return
# noinspection PyPep8Naming
[docs] def create_affine_xyz(self):
"""Create affine in xyz.
"""
def normalize(v):
"""Normalize a vector
https://stackoverflow.com/questions/21030391/how-to-normalize-an-array-in-numpy
Args:
v: 3D vector
Returns:
normalized 3D vector
"""
norm = np.linalg.norm(v, ord=1)
if norm == 0:
norm = np.finfo(v.dtype).eps
return v / norm
ds, dr, dc = self.spacing
colr = normalize(np.array(self.orientation[3:6])).reshape((3,)) * [-1, -1, 1]
colc = normalize(np.array(self.orientation[0:3])).reshape((3,)) * [-1, -1, 1]
# T0 = self.imagePositions[0][::-1].reshape(3, ) # x,y,z
if self.slices > 1:
# Tn = self.imagePositions[self.slices - 1][::-1].reshape(3, ) # x,y,z
# k = Tn
k = np.cross(colc, colr, axis=0)
k = k * ds
else:
k = np.cross(colc, colr, axis=0)
k = k * ds
L = np.zeros((4, 4))
L[:3, 1] = colr * dr
L[:3, 0] = colc * dc
L[:3, 2] = k
L[:3, 3] = self.origin * [-1, -1, 1]
L[3, 3] = 1
return L
# def getQformFromTransformationMatrix(self):
# # def matrix_from_orientation(orientation, normal):
# # oT = orientation.reshape((2,3)).T
# # colr = oT[:,0].reshape((3,1))
# # colc = oT[:,1].reshape((3,1))
# # coln = normal.reshape((3,1))
# # if len(self.shape) < 3:
# # M = np.hstack((colr[:2], colc[:2])).reshape((2,2))
# # else:
# # M = np.hstack((colr, colc, coln)).reshape((3,3))
# # return M
#
# def normalize(v):
# """Normalize a vector
#
# https://stackoverflow.com/questions/21030391/how-to-normalize-an-array-in-numpy
#
# :param v: 3D vector
# :return: normalized 3D vector
# """
# norm = np.linalg.norm(v, ord=1)
# if norm == 0:
# norm = np.finfo(v.dtype).eps
# return v / norm
#
# def L_from_orientation(orientation, normal, spacing):
# """
# orientation: row, then column index direction cosines
# """
# _ds, _dr, _dc = spacing
# _colr = normalize(np.array(orientation[3:6])).reshape((3,))
# _colc = normalize(np.array(orientation[0:3])).reshape((3,))
# _t0 = self.imagePositions[0][::-1].reshape(3, ) # x,y,z
# if self.slices > 1:
# _tn = self.imagePositions[self.slices - 1][::-1].reshape(3, ) # x,y,z
# # k = _tn
# _k = np.cross(_colr, _colc, axis=0)
# _k = _k * _ds
# else:
# _k = np.cross(_colr, _colc, axis=0)
# _k = _k * _ds
#
# _L = np.zeros((4, 4))
# _L[:3, 0] = _t0[:]
# _L[3, 0] = 1
# _L[:3, 1] = _k
# _L[3, 1] = 1 if self.slices > 1 else 0
# _L[:3, 2] = _colr * [-1, -1, 1] * _dr
# _L[:3, 3] = _colc * [-1, -1, 1] * _dc
# return _L
#
# # M = self.transformationMatrix
# # M = matrix_from_orientation(self.orientation, self.normal)
# # ipp = self.origin
# # q = np.array([[M[2,2], M[2,1], M[2,0], ipp[0]],
# # [M[1,2], M[1,1], M[1,0], ipp[1]],
# # [M[0,2], M[0,1], M[0,0], ipp[2]],
# # [ 0, 0, 0, 1 ]]
# # )
#
# if self.slices > 1:
# r = np.array([[1, 1, 1, 0], [1, 1, 0, 1], [1, self.slices, 0, 0], [1, 1, 0, 0]])
# else:
# r = np.array([[1, 0, 1, 0], [1, 0, 0, 1], [1, self.slices, 0, 0], [1, 0, 0, 0]])
# l = L_from_orientation(self.orientation, self.normal, self.spacing)
#
# # Linv = np.linalg.inv(L)
# # Aspm = np.dot(r, np.linalg.inv(l))
# to_ones = np.eye(4)
# to_ones[:, 3] = 1
# # A = np.dot(Aspm, to_ones)
#
# ds, dr, dc = self.spacing
# colr = normalize(np.array(self.orientation[3:6])).reshape((3,))
# colc = normalize(np.array(self.orientation[0:3])).reshape((3,))
# coln = normalize(np.cross(colc, colr, axis=0))
# t_0 = self.imagePositions[0][::-1].reshape(3, ) # x,y,z
# if self.slices > 1:
# t_n = self.imagePositions[self.slices - 1][::-1].reshape((3,)) # x,y,z
# abcd = np.array([1, 1, self.slices, 1]).reshape((4,))
# one = np.ones((1,))
# efgh = np.concatenate((t_n, one))
# else:
# abcd = np.array([0, 0, 1, 0]).reshape((4,))
# # zero = np.zeros((1,))
# efgh = np.concatenate((n * ds, zeros))
#
# # From derivations/spm_dicom_orient.py
#
# # premultiplication matrix to go from 0 to 1 based indexing
# one_based = np.eye(4)
# one_based[:3, 3] = (1, 1, 1)
# # premult for swapping row and column indices
# row_col_swap = np.eye(4)
# row_col_swap[:, 0] = np.eye(4)[:, 1]
# row_col_swap[:, 1] = np.eye(4)[:, 0]
#
# # various worming matrices
# orient_pat = np.hstack([colr.reshape(3, 1), colc.reshape(3, 1)])
# orient_cross = coln
# pos_pat_0 = t_0
# if self.slices > 1:
# missing_r_col = (t_0 - t_n) / (1 - self.slices)
# pos_pat_N = t_n
# pixel_spacing = [dr, dc]
# NZ = self.slices
# slice_thickness = ds
#
# R3 = np.dot(orient_pat, np.diag(pixel_spacing))
# # R3 = orient_pat * np.diag(pixel_spacing)
# r = np.zeros((4, 2))
# r[:3, :] = R3
#
# # The following is specific to the SPM algorithm.
# x1 = np.ones(4)
# y1 = np.ones(4)
# y1[:3] = pos_pat_0
#
# to_inv = np.zeros((4, 4))
# to_inv[:, 0] = x1
# to_inv[:, 1] = abcd
# to_inv[0, 2] = 1
# to_inv[1, 3] = 1
# inv_lhs = np.zeros((4, 4))
# inv_lhs[:, 0] = y1
# inv_lhs[:, 1] = efgh
# inv_lhs[:, 2:] = r
#
# def spm_full_matrix(x2, y2):
# rhs = to_inv[:, :]
# rhs[:, 1] = x2
# lhs = inv_lhs[:, :]
# lhs[:, 1] = y2
# return np.dot(lhs, np.linalg.inv(rhs))
#
# if self.slices > 1:
# x2_ms = np.array([1, 1, NZ, 1])
# y2_ms = np.ones((4,))
# y2_ms[:3] = pos_pat_N
# A_ms = spm_full_matrix(x2_ms, y2_ms)
# A = A_ms
# else:
# orient = np.zeros((3, 3))
# orient[:3, :2] = orient_pat
# orient[:, 2] = orient_cross
# x2_ss = np.array([0, 0, 1, 0])
# y2_ss = np.zeros((4,))
# # y2_ss[:3] = orient * np.array([0, 0, slice_thickness])
# y2_ss[:3] = np.dot(orient, np.array([0, 0, slice_thickness]))
# A_ss = spm_full_matrix(x2_ss, y2_ss)
# A = A_ss
#
# A = np.dot(A, row_col_swap)
#
# multi_aff = np.eye(4)
# multi_aff[:3, :2] = R3
# trans_z_N = np.array([0, 0, self.slices - 1, 1])
# multi_aff[:3, 2] = missing_r_col
# multi_aff[:3, 3] = pos_pat_0
# # est_pos_pat_N = np.dot(multi_aff, trans_z_N)
#
# # Flip voxels in y
# analyze_to_dicom = np.eye(4)
# analyze_to_dicom[1, 1] = -1
# # analyze_to_dicom[1,3] = shape[1]+1
# analyze_to_dicom[1, 3] = self.slices
# logging.debug("getQformFromTransformationMatrix: analyze_to_dicom\n{}".format(analyze_to_dicom))
# # dicom_to_analyze = np.linalg.inv(analyze_to_dicom)
# # q = np.dot(q,dicom_to_analyze)
# q = np.dot(A, analyze_to_dicom)
# # ## 2019.07.03 # q = np.dot(q,analyze_to_dicom)
# # ## 2019.07.03 # logging.debug("q after rows dicom_to_analyze:\n{}".format(q))
# # Flip mm coords in x and y directions
# patient_to_tal = np.diag([1, -1, -1, 1])
# # patient_to_tal = np.eye(4)
# # patient_to_tal[0,0] = -1
# # patient_to_tal[1,1] = -1
# # tal_to_patient = np.linalg.inv(patient_to_tal)
# # q = np.dot(tal_to_patient,q)
# logging.debug("getQformFromTransformationMatrix: patient_to_tal\n{}".format(patient_to_tal))
# q = np.dot(patient_to_tal, q)
# logging.debug("getQformFromTransformationMatrix: q after\n{}".format(q))
#
# return q
# def create_affine(self, sorted_dicoms):
# """
# Function to generate the affine matrix for a dicom series
# From dicom2nifti:common.py: https://github.com/icometrix/dicom2nifti/blob/master/dicom2nifti/common.py
# This method was based on (http://nipy.org/nibabel/dicom/dicom_orientation.html)
# :param sorted_dicoms: list with sorted dicom files
# """
#
# # Create affine matrix (http://nipy.sourceforge.net/nibabel/dicom/dicom_orientation.html#dicom-slice-affine)
# image_orient1 = np.array(sorted_dicoms[0].ImageOrientationPatient)[0:3]
# image_orient2 = np.array(sorted_dicoms[0].ImageOrientationPatient)[3:6]
#
# delta_r = float(sorted_dicoms[0].PixelSpacing[0])
# delta_c = float(sorted_dicoms[0].PixelSpacing[1])
#
# image_pos = np.array(sorted_dicoms[0].ImagePositionPatient)
#
# last_image_pos = np.array(sorted_dicoms[-1].ImagePositionPatient)
#
# if len(sorted_dicoms) == 1:
# # Single slice
# step = [0, 0, -1]
# else:
# step = (image_pos - last_image_pos) / (1 - len(sorted_dicoms))
#
# # check if this is actually a volume and not all slices on the same location
# if np.linalg.norm(step) == 0.0:
# raise imagedata.formats.NotImageError("Not a volume")
#
# affine = np.array(
# [[-image_orient1[0] * delta_c, -image_orient2[0] * delta_r, -step[0], -image_pos[0]],
# [-image_orient1[1] * delta_c, -image_orient2[1] * delta_r, -step[1], -image_pos[1]],
# [image_orient1[2] * delta_c, image_orient2[2] * delta_r, step[2], image_pos[2]],
# [0, 0, 0, 1]]
# )
# return affine, np.linalg.norm(step)
[docs] def write_3d_numpy(self, si, destination, opts):
"""Write 3D numpy image as Nifti file
Args:
self: NiftiPlugin instance
si: Series array (3D or 4D), including these attributes:
- slices,
- spacing,
- imagePositions,
- transformationMatrix,
- orientation,
- tags
destination: dict of archive and filenames
opts: Output options (dict)
"""
if si.color:
raise imagedata.formats.WriteNotImplemented(
"Writing color Nifti images not implemented.")
logging.debug('NiftiPlugin.write_3d_numpy: destination {}'.format(destination))
archive = destination['archive']
filename_template = 'Image.nii.gz'
if len(destination['files']) > 0 and len(destination['files'][0]) > 0:
filename_template = destination['files'][0]
self.shape = si.shape
self.slices = si.slices
self.spacing = si.spacing
self.transformationMatrix = si.transformationMatrix
self.imagePositions = si.imagePositions
self.tags = si.tags
self.origin, self.orientation, self.normal = si.get_transformation_components_xyz()
logging.info("Data shape write: {}".format(imagedata.formats.shape_to_str(si.shape)))
assert si.ndim == 2 or si.ndim == 3, "write_3d_series: input dimension %d is not 3D." % si.ndim
fsi = self._reorder_from_dicom(si, flip=False, flipud=True)
shape = fsi.shape
affine_xyz = self.create_affine_xyz()
nifti_header = nibabel.Nifti1Header()
nifti_header.set_dim_info(freq=0, phase=1, slice=2)
nifti_header.set_data_shape(shape)
dz, dy, dx = self.spacing
if si.ndim < 3:
nifti_header.set_zooms((dx, dy))
else:
nifti_header.set_zooms((dx, dy, dz))
nifti_header.set_data_dtype(fsi.dtype)
nifti_header.set_sform(affine_xyz, code=1)
nifti_header.set_xyzt_units(xyz='mm')
img = nibabel.Nifti1Image(fsi, None, nifti_header)
try:
filename = filename_template % 0
except TypeError:
filename = filename_template
self.write_numpy_nifti(img, archive, filename)
[docs] def write_4d_numpy(self, si, destination, opts):
"""Write 4D numpy image as Nifti file
Args:
self: NiftiPlugin instance
si[tag,slice,rows,columns]: Series array, including these attributes:
- slices,
- spacing,
- imagePositions,
- transformationMatrix,
- orientation,
- tags
destination: dict of archive and filenames
opts: Output options (dict)
"""
if si.color:
raise imagedata.formats.WriteNotImplemented(
"Writing color Nifti images not implemented.")
logging.debug('ITKPlugin.write_4d_numpy: destination {}'.format(destination))
archive = destination['archive']
filename_template = 'Image.nii.gz'
if len(destination['files']) > 0 and len(destination['files'][0]) > 0:
filename_template = destination['files'][0]
self.shape = si.shape
self.slices = si.slices
self.spacing = si.spacing
self.transformationMatrix = si.transformationMatrix
self.imagePositions = si.imagePositions
self.tags = si.tags
self.origin, self.orientation, self.normal = si.get_transformation_components_xyz()
# Defaults
self.output_sort = imagedata.formats.SORT_ON_SLICE
if 'output_sort' in opts:
self.output_sort = opts['output_sort']
# Should we allow to write 3D volume?
if si.ndim == 2:
si.shape = (1, 1,) + si.shape
elif si.ndim == 3:
si.shape = (1,) + si.shape
if si.ndim != 4:
raise ValueError("write_4d_numpy: input dimension {} is not 4D.".format(si.ndim))
logging.debug("write_4d_numpy: si dtype {}, shape {}, sort {}".format(
si.dtype, si.shape,
imagedata.formats.sort_on_to_str(self.output_sort)))
steps = si.shape[0]
slices = si.shape[1]
if steps != len(si.tags[0]):
raise ValueError(
"write_4d_series: tags of dicom template ({}) differ from input array ({}).".format(len(si.tags[0]),
steps))
if slices != si.slices:
raise ValueError(
"write_4d_series: slices of dicom template ({}) differ from input array ({}).".format(si.slices,
slices))
fsi = self._reorder_from_dicom(si, flip=False, flipud=True)
shape = fsi.shape
affine_xyz = self.create_affine_xyz()
nifti_header = nibabel.Nifti1Header()
nifti_header.set_dim_info(freq=0, phase=1, slice=2)
nifti_header.set_data_shape(shape)
dz, dy, dx = self.spacing
nifti_header.set_zooms((dx, dy, dz, 1))
nifti_header.set_data_dtype(fsi.dtype)
nifti_header.set_sform(affine_xyz, code=1)
# NiftiHeader.set_slice_duration()
# NiftiHeader.set_slice_times(times)
nifti_header.set_xyzt_units(xyz='mm', t='sec')
img = nibabel.Nifti1Image(fsi, None, nifti_header)
try:
filename = filename_template % 0
except TypeError:
filename = filename_template
self.write_numpy_nifti(img, archive, filename)
[docs] @staticmethod
def write_numpy_nifti(img, archive, filename):
"""Write nifti data to file
Args:
self: ITKPlugin instance, including these attributes:
- slices (not used)
- spacing
- imagePositions
- transformationMatrix
- orientation (not used)
- tags (not used)
img: Nifti1Image
archive: archive object
filename: file name, possibly without extentsion
"""
if len(os.path.splitext(filename)[1]) == 0:
filename = filename + '.nii.gz'
ext = os.path.splitext(filename)[1]
if filename.endswith('.nii.gz'):
ext = '.nii.gz'
logging.debug('write_numpy_nifti: ext %s' % ext)
f = tempfile.NamedTemporaryFile(
suffix=ext, delete=False)
logging.debug('write_numpy_nifti: write local file %s' % f.name)
img.to_filename(f.name)
f.close()
logging.debug('write_numpy_nifti: copy to file %s' % filename)
_ = archive.add_localfile(f.name, filename)
os.unlink(f.name)