From ece0bd63e9b5e43d8610d5050ecf6a7f9b7e9038 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 7 Dec 2017 16:25:36 +1030
Subject: [PATCH] Including a copy of dcmstack, which will be removed when it
is brought up to date.
---
fsl/data/dcmstack/COPYING | 35 +
fsl/data/dcmstack/__init__.py | 7 +
fsl/data/dcmstack/dcmmeta.py | 1804 +++++++++++++++++++++++++++++
fsl/data/dcmstack/dcmstack.py | 1162 +++++++++++++++++++
fsl/data/dcmstack/dcmstack_cli.py | 357 ++++++
fsl/data/dcmstack/extract.py | 520 +++++++++
fsl/data/dcmstack/info.py | 54 +
fsl/data/dcmstack/nitool_cli.py | 251 ++++
8 files changed, 4190 insertions(+)
create mode 100644 fsl/data/dcmstack/COPYING
create mode 100644 fsl/data/dcmstack/__init__.py
create mode 100644 fsl/data/dcmstack/dcmmeta.py
create mode 100644 fsl/data/dcmstack/dcmstack.py
create mode 100755 fsl/data/dcmstack/dcmstack_cli.py
create mode 100644 fsl/data/dcmstack/extract.py
create mode 100644 fsl/data/dcmstack/info.py
create mode 100644 fsl/data/dcmstack/nitool_cli.py
diff --git a/fsl/data/dcmstack/COPYING b/fsl/data/dcmstack/COPYING
new file mode 100644
index 000000000..7374d5d05
--- /dev/null
+++ b/fsl/data/dcmstack/COPYING
@@ -0,0 +1,35 @@
+**********************
+Copyright and Licenses
+**********************
+
+dcmstack
+=======
+
+The dcmstack package, including all examples, code snippets and attached
+documentation is covered by the MIT license.
+
+::
+
+ The MIT License
+
+ Copyright (c) 2011-2012 Brendan Moloney <moloney@ohsu.edu>
+
+ Permission is hereby granted, free of charge, to any person obtaining a copy
+ of this software and associated documentation files (the "Software"), to deal
+ in the Software without restriction, including without limitation the rights
+ to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ copies of the Software, and to permit persons to whom the Software is
+ furnished to do so, subject to the following conditions:
+
+ The above copyright notice and this permission notice shall be included in
+ all copies or substantial portions of the Software.
+
+ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ THE SOFTWARE.
+
+
diff --git a/fsl/data/dcmstack/__init__.py b/fsl/data/dcmstack/__init__.py
new file mode 100644
index 000000000..2663c2259
--- /dev/null
+++ b/fsl/data/dcmstack/__init__.py
@@ -0,0 +1,7 @@
+"""
+Package for stacking DICOM images into multi dimensional volumes, extracting
+the DICOM meta data, converting the result to Nifti files with the meta
+data stored in a header extension, and work with these extended Nifti files.
+"""
+from .info import __version__
+from .dcmstack import *
diff --git a/fsl/data/dcmstack/dcmmeta.py b/fsl/data/dcmstack/dcmmeta.py
new file mode 100644
index 000000000..9049d7d7c
--- /dev/null
+++ b/fsl/data/dcmstack/dcmmeta.py
@@ -0,0 +1,1804 @@
+"""
+DcmMeta header extension and NiftiWrapper for working with extended Niftis.
+"""
+from __future__ import print_function
+
+import sys
+import json, warnings
+from copy import deepcopy
+import numpy as np
+import nibabel as nb
+from nibabel.nifti1 import Nifti1Extension
+from nibabel.spatialimages import HeaderDataError
+
+try:
+ from collections import OrderedDict
+except ImportError:
+ from ordereddict import OrderedDict
+
+with warnings.catch_warnings():
+ warnings.simplefilter('ignore')
+ from nibabel.nicom.dicomwrappers import wrapper_from_data
+
+dcm_meta_ecode = 0
+
+_meta_version = 0.6
+
+_req_base_keys_map= {0.5 : set(('dcmmeta_affine',
+ 'dcmmeta_slice_dim',
+ 'dcmmeta_shape',
+ 'dcmmeta_version',
+ 'global',
+ )
+ ),
+ 0.6 : set(('dcmmeta_affine',
+ 'dcmmeta_reorient_transform',
+ 'dcmmeta_slice_dim',
+ 'dcmmeta_shape',
+ 'dcmmeta_version',
+ 'global',
+ )
+ ),
+ }
+'''Minimum required keys in the base dictionaty to be considered valid'''
+
+def is_constant(sequence, period=None):
+ '''Returns true if all elements in (each period of) the sequence are equal.
+
+ Parameters
+ ----------
+ sequence : sequence
+ The sequence of elements to check.
+
+ period : int
+ If not None then each subsequence of that length is checked.
+ '''
+ if period is None:
+ return all(val == sequence[0] for val in sequence)
+ else:
+ if period <= 1:
+ raise ValueError('The period must be greater than one')
+ seq_len = len(sequence)
+ if seq_len % period != 0:
+ raise ValueError('The sequence length is not evenly divisible by '
+ 'the period length.')
+
+ for period_idx in range(seq_len / period):
+ start_idx = period_idx * period
+ end_idx = start_idx + period
+ if not all(val == sequence[start_idx]
+ for val in sequence[start_idx:end_idx]):
+ return False
+
+ return True
+
+def is_repeating(sequence, period):
+ '''Returns true if the elements in the sequence repeat with the given
+ period.
+
+ Parameters
+ ----------
+ sequence : sequence
+ The sequence of elements to check.
+
+ period : int
+ The period over which the elements should repeat.
+ '''
+ seq_len = len(sequence)
+ if period <= 1 or period >= seq_len:
+ raise ValueError('The period must be greater than one and less than '
+ 'the length of the sequence')
+ if seq_len % period != 0:
+ raise ValueError('The sequence length is not evenly divisible by the '
+ 'period length.')
+
+ for period_idx in range(1, seq_len / period):
+ start_idx = period_idx * period
+ end_idx = start_idx + period
+ if sequence[start_idx:end_idx] != sequence[:period]:
+ return False
+
+ return True
+
+class InvalidExtensionError(Exception):
+ def __init__(self, msg):
+ '''Exception denoting than a DcmMetaExtension is invalid.'''
+ self.msg = msg
+
+ def __str__(self):
+ return 'The extension is not valid: %s' % self.msg
+
+
+class DcmMetaExtension(Nifti1Extension):
+ '''Nifti extension for storing a summary of the meta data from the source
+ DICOM files.
+ '''
+
+ @property
+ def reorient_transform(self):
+ '''The transformation due to reorientation of the data array. Can be
+ used to update directional DICOM meta data (after converting to RAS if
+ needed) into the same space as the affine.'''
+ if self.version < 0.6:
+ return None
+ if self._content['dcmmeta_reorient_transform'] is None:
+ return None
+ return np.array(self._content['dcmmeta_reorient_transform'])
+
+ @reorient_transform.setter
+ def reorient_transform(self, value):
+ if not value is None and value.shape != (4, 4):
+ raise ValueError("The reorient_transform must be none or (4,4) "
+ "array")
+ if value is None:
+ self._content['dcmmeta_reorient_transform'] = None
+ else:
+ self._content['dcmmeta_reorient_transform'] = value.tolist()
+
+ @property
+ def affine(self):
+ '''The affine associated with the meta data. If this differs from the
+ image affine, the per-slice meta data will not be used. '''
+ return np.array(self._content['dcmmeta_affine'])
+
+ @affine.setter
+ def affine(self, value):
+ if value.shape != (4, 4):
+ raise ValueError("Invalid shape for affine")
+ self._content['dcmmeta_affine'] = value.tolist()
+
+ @property
+ def slice_dim(self):
+ '''The index of the slice dimension associated with the per-slice
+ meta data.'''
+ return self._content['dcmmeta_slice_dim']
+
+ @slice_dim.setter
+ def slice_dim(self, value):
+ if not value is None and not (0 <= value < 3):
+ raise ValueError("The slice dimension must be between zero and "
+ "two")
+ self._content['dcmmeta_slice_dim'] = value
+
+ @property
+ def shape(self):
+ '''The shape of the data associated with the meta data. Defines the
+ number of values for the meta data classifications.'''
+ return tuple(self._content['dcmmeta_shape'])
+
+ @shape.setter
+ def shape(self, value):
+ if not (3 <= len(value) < 6):
+ raise ValueError("The shape must have a length between three and "
+ "six")
+ self._content['dcmmeta_shape'][:] = value
+
+ @property
+ def version(self):
+ '''The version of the meta data extension.'''
+ return self._content['dcmmeta_version']
+
+ @version.setter
+ def version(self, value):
+ '''Set the version of the meta data extension.'''
+ self._content['dcmmeta_version'] = value
+
+ @property
+ def slice_normal(self):
+ '''The slice normal associated with the per-slice meta data.'''
+ slice_dim = self.slice_dim
+ if slice_dim is None:
+ return None
+ return np.array(self.affine[slice_dim][:3])
+
+ @property
+ def n_slices(self):
+ '''The number of slices associated with the per-slice meta data.'''
+ slice_dim = self.slice_dim
+ if slice_dim is None:
+ return None
+ return self.shape[slice_dim]
+
+ classifications = (('global', 'const'),
+ ('global', 'slices'),
+ ('time', 'samples'),
+ ('time', 'slices'),
+ ('vector', 'samples'),
+ ('vector', 'slices'),
+ )
+ '''The classifications used to separate meta data based on if and how the
+ values repeat. Each class is a tuple with a base class and a sub class.'''
+
+ def get_valid_classes(self):
+ '''Return the meta data classifications that are valid for this
+ extension.
+
+ Returns
+ -------
+ valid_classes : tuple
+ The classifications that are valid for this extension (based on its
+ shape).
+
+ '''
+ shape = self.shape
+ n_dims = len(shape)
+ if n_dims == 3:
+ return self.classifications[:2]
+ elif n_dims == 4:
+ return self.classifications[:4]
+ elif n_dims == 5:
+ if shape[3] != 1:
+ return self.classifications
+ else:
+ return self.classifications[:2] + self.classifications[-2:]
+ else:
+ raise ValueError("There must be 3 to 5 dimensions.")
+
+ def get_multiplicity(self, classification):
+ '''Get the number of meta data values for all meta data of the provided
+ classification.
+
+ Parameters
+ ----------
+ classification : tuple
+ The meta data classification.
+
+ Returns
+ -------
+ multiplicity : int
+ The number of values for any meta data of the provided
+ `classification`.
+ '''
+ if not classification in self.get_valid_classes():
+ raise ValueError("Invalid classification: %s" % classification)
+
+ base, sub = classification
+ shape = self.shape
+ n_vals = 1
+ if sub == 'slices':
+ n_vals = self.n_slices
+ if n_vals is None:
+ return 0
+ if base == 'vector':
+ n_vals *= shape[3]
+ elif base == 'global':
+ for dim_size in shape[3:]:
+ n_vals *= dim_size
+ elif sub == 'samples':
+ if base == 'time':
+ n_vals = shape[3]
+ if len(shape) == 5:
+ n_vals *= shape[4]
+ elif base == 'vector':
+ n_vals = shape[4]
+
+ return n_vals
+
+ def check_valid(self):
+ '''Check if the extension is valid.
+
+ Raises
+ ------
+ InvalidExtensionError
+ The extension is missing required meta data or classifications, or
+ some element(s) have the wrong number of values for their
+ classification.
+ '''
+ #Check for the required base keys in the json data
+ if not _req_base_keys_map[self.version] <= set(self._content):
+ raise InvalidExtensionError('Missing one or more required keys')
+
+ #Check the orientation/shape/version
+ if self.affine.shape != (4, 4):
+ raise InvalidExtensionError('Affine has incorrect shape')
+ slice_dim = self.slice_dim
+ if slice_dim is not None:
+ if not (0 <= slice_dim < 3):
+ raise InvalidExtensionError('Slice dimension is not valid')
+ if not (3 <= len(self.shape) < 6):
+ raise InvalidExtensionError('Shape is not valid')
+
+ #Check all required meta dictionaries, make sure values have correct
+ #multiplicity
+ valid_classes = self.get_valid_classes()
+ for classes in valid_classes:
+ if not classes[0] in self._content:
+ raise InvalidExtensionError('Missing required base '
+ 'classification %s' % classes[0])
+ if not classes[1] in self._content[classes[0]]:
+ raise InvalidExtensionError(('Missing required sub '
+ 'classification %s in base '
+ 'classification %s') % classes)
+ cls_meta = self.get_class_dict(classes)
+ cls_mult = self.get_multiplicity(classes)
+ if cls_mult == 0 and len(cls_meta) != 0:
+ raise InvalidExtensionError('Slice dim is None but per-slice '
+ 'meta data is present')
+ elif cls_mult > 1:
+ for key, vals in cls_meta.iteritems():
+ n_vals = len(vals)
+ if n_vals != cls_mult:
+ msg = (('Incorrect number of values for key %s with '
+ 'classification %s, expected %d found %d') %
+ (key, classes, cls_mult, n_vals)
+ )
+ raise InvalidExtensionError(msg)
+
+ #Check that all keys are uniquely classified
+ for classes in valid_classes:
+ for other_classes in valid_classes:
+ if classes == other_classes:
+ continue
+ intersect = (set(self.get_class_dict(classes)) &
+ set(self.get_class_dict(other_classes))
+ )
+ if len(intersect) != 0:
+ raise InvalidExtensionError("One or more keys have "
+ "multiple classifications")
+
+ def get_keys(self):
+ '''Get a list of all the meta data keys that are available.'''
+ keys = []
+ for base_class, sub_class in self.get_valid_classes():
+ keys += self._content[base_class][sub_class].keys()
+ return keys
+
+ def get_classification(self, key):
+ '''Get the classification for the given `key`.
+
+ Parameters
+ ----------
+ key : str
+ The meta data key.
+
+ Returns
+ -------
+ classification : tuple or None
+ The classification tuple for the provided key or None if the key is
+ not found.
+
+ '''
+ for base_class, sub_class in self.get_valid_classes():
+ if key in self._content[base_class][sub_class]:
+ return (base_class, sub_class)
+
+ return None
+
+ def get_class_dict(self, classification):
+ '''Get the dictionary for the given classification.
+
+ Parameters
+ ----------
+ classification : tuple
+ The meta data classification.
+
+ Returns
+ -------
+ meta_dict : dict
+ The dictionary for the provided classification.
+ '''
+ base, sub = classification
+ return self._content[base][sub]
+
+ def get_values(self, key):
+ '''Get all values for the provided key.
+
+ Parameters
+ ----------
+ key : str
+ The meta data key.
+
+ Returns
+ -------
+ values
+ The value or values for the given key. The number of values
+ returned depends on the classification (see 'get_multiplicity').
+ '''
+ classification = self.get_classification(key)
+ if classification is None:
+ return None
+ return self.get_class_dict(classification)[key]
+
+ def get_values_and_class(self, key):
+ '''Get the values and the classification for the provided key.
+
+ Parameters
+ ----------
+ key : str
+ The meta data key.
+
+ Returns
+ -------
+ vals_and_class : tuple
+ None for both the value and classification if the key is not found.
+
+ '''
+ classification = self.get_classification(key)
+ if classification is None:
+ return (None, None)
+ return (self.get_class_dict(classification)[key], classification)
+
+ def filter_meta(self, filter_func):
+ '''Filter the meta data.
+
+ Parameters
+ ----------
+ filter_func : callable
+ Must take a key and values as parameters and return True if they
+ should be filtered out.
+
+ '''
+ for classes in self.get_valid_classes():
+ filtered = []
+ curr_dict = self.get_class_dict(classes)
+ for key, values in curr_dict.iteritems():
+ if filter_func(key, values):
+ filtered.append(key)
+ for key in filtered:
+ del curr_dict[key]
+
+ def clear_slice_meta(self):
+ '''Clear all meta data that is per slice.'''
+ for base_class, sub_class in self.get_valid_classes():
+ if sub_class == 'slices':
+ self.get_class_dict((base_class, sub_class)).clear()
+
+ def get_subset(self, dim, idx):
+ '''Get a DcmMetaExtension containing a subset of the meta data.
+
+ Parameters
+ ----------
+ dim : int
+ The dimension we are taking the subset along.
+
+ idx : int
+ The position on the dimension `dim` for the subset.
+
+ Returns
+ -------
+ result : DcmMetaExtension
+ A new DcmMetaExtension corresponding to the subset.
+
+ '''
+ if not 0 <= dim < 5:
+ raise ValueError("The argument 'dim' must be in the range [0, 5).")
+
+ shape = self.shape
+ valid_classes = self.get_valid_classes()
+
+ #Make an empty extension for the result
+ result_shape = list(shape)
+ result_shape[dim] = 1
+ while result_shape[-1] == 1 and len(result_shape) > 3:
+ result_shape = result_shape[:-1]
+ result = self.make_empty(result_shape,
+ self.affine,
+ self.reorient_transform,
+ self.slice_dim
+ )
+
+ for src_class in valid_classes:
+ #Constants remain constant
+ if src_class == ('global', 'const'):
+ for key, val in self.get_class_dict(src_class).iteritems():
+ result.get_class_dict(src_class)[key] = deepcopy(val)
+ continue
+
+ if dim == self.slice_dim:
+ if src_class[1] != 'slices':
+ for key, vals in self.get_class_dict(src_class).iteritems():
+ result.get_class_dict(src_class)[key] = deepcopy(vals)
+ else:
+ result._copy_slice(self, src_class, idx)
+ elif dim < 3:
+ for key, vals in self.get_class_dict(src_class).iteritems():
+ result.get_class_dict(src_class)[key] = deepcopy(vals)
+ elif dim == 3:
+ result._copy_sample(self, src_class, 'time', idx)
+ else:
+ result._copy_sample(self, src_class, 'vector', idx)
+
+ return result
+
+ def to_json(self):
+ '''Return the extension encoded as a JSON string.'''
+ self.check_valid()
+ return self._mangle(self._content)
+
+ @classmethod
+ def from_json(klass, json_str):
+ '''Create an extension from the JSON string representation.'''
+ result = klass(dcm_meta_ecode, json_str)
+ result.check_valid()
+ return result
+
+ @classmethod
+ def make_empty(klass, shape, affine, reorient_transform=None,
+ slice_dim=None):
+ '''Make an empty DcmMetaExtension.
+
+ Parameters
+ ----------
+ shape : tuple
+ The shape of the data associated with this extension.
+
+ affine : array
+ The RAS affine for the data associated with this extension.
+
+ reorient_transform : array
+ The transformation matrix representing any reorientation of the
+ data array.
+
+ slice_dim : int
+ The index of the slice dimension for the data associated with this
+ extension
+
+ Returns
+ -------
+ result : DcmMetaExtension
+ An empty DcmMetaExtension with the required values set to the
+ given arguments.
+
+ '''
+ result = klass(dcm_meta_ecode, '{}')
+ result._content['global'] = OrderedDict()
+ result._content['global']['const'] = OrderedDict()
+ result._content['global']['slices'] = OrderedDict()
+
+ if len(shape) > 3 and shape[3] != 1:
+ result._content['time'] = OrderedDict()
+ result._content['time']['samples'] = OrderedDict()
+ result._content['time']['slices'] = OrderedDict()
+
+ if len(shape) > 4:
+ result._content['vector'] = OrderedDict()
+ result._content['vector']['samples'] = OrderedDict()
+ result._content['vector']['slices'] = OrderedDict()
+
+ result._content['dcmmeta_shape'] = []
+ result.shape = shape
+ result.affine = affine
+ result.reorient_transform = reorient_transform
+ result.slice_dim = slice_dim
+ result.version = _meta_version
+
+ return result
+
+ @classmethod
+ def from_runtime_repr(klass, runtime_repr):
+ '''Create an extension from the Python runtime representation (nested
+ dictionaries).
+ '''
+ result = klass(dcm_meta_ecode, '{}')
+ result._content = runtime_repr
+ result.check_valid()
+ return result
+
+ @classmethod
+ def from_sequence(klass, seq, dim, affine=None, slice_dim=None):
+ '''Create an extension from a sequence of extensions.
+
+ Parameters
+ ----------
+ seq : sequence
+ The sequence of DcmMetaExtension objects.
+
+ dim : int
+ The dimension to merge the extensions along.
+
+ affine : array
+ The affine to use in the resulting extension. If None, the affine
+ from the first extension in `seq` will be used.
+
+ slice_dim : int
+ The slice dimension to use in the resulting extension. If None, the
+ slice dimension from the first extension in `seq` will be used.
+
+ Returns
+ -------
+ result : DcmMetaExtension
+ The result of merging the extensions in `seq` along the dimension
+ `dim`.
+ '''
+ if not 0 <= dim < 5:
+ raise ValueError("The argument 'dim' must be in the range [0, 5).")
+
+ n_inputs = len(seq)
+ first_input = seq[0]
+ input_shape = first_input.shape
+
+ if len(input_shape) > dim and input_shape[dim] != 1:
+ raise ValueError("The dim must be singular or not exist for the "
+ "inputs.")
+
+ output_shape = list(input_shape)
+ while len(output_shape) <= dim:
+ output_shape.append(1)
+ output_shape[dim] = n_inputs
+
+ if affine is None:
+ affine = first_input.affine
+ if slice_dim is None:
+ slice_dim = first_input.slice_dim
+
+ result = klass.make_empty(output_shape,
+ affine,
+ None,
+ slice_dim)
+
+ #Need to initialize the result with the first extension in 'seq'
+ result_slc_norm = result.slice_normal
+ first_slc_norm = first_input.slice_normal
+ use_slices = (not result_slc_norm is None and
+ not first_slc_norm is None and
+ np.allclose(result_slc_norm, first_slc_norm))
+ for classes in first_input.get_valid_classes():
+ if classes[1] == 'slices' and not use_slices:
+ continue
+ result._content[classes[0]][classes[1]] = \
+ deepcopy(first_input.get_class_dict(classes))
+
+ #Adjust the shape to what the extension actually contains
+ shape = list(result.shape)
+ shape[dim] = 1
+ result.shape = shape
+
+ #Initialize reorient transform
+ reorient_transform = first_input.reorient_transform
+
+ #Add the other extensions, updating the shape as we go
+ for input_ext in seq[1:]:
+ #If the affines or reorient_transforms don't match, we set the
+ #reorient_transform to None as we can not reliably use it to update
+ #directional meta data
+ if ((reorient_transform is None or
+ input_ext.reorient_transform is None) or
+ not (np.allclose(input_ext.affine, affine) or
+ np.allclose(input_ext.reorient_transform,
+ reorient_transform)
+ )
+ ):
+ reorient_transform = None
+ result._insert(dim, input_ext)
+ shape[dim] += 1
+ result.shape = shape
+
+ #Set the reorient transform
+ result.reorient_transform = reorient_transform
+
+ #Try simplifying any keys in global slices
+ for key in result.get_class_dict(('global', 'slices')).keys():
+ result._simplify(key)
+
+ return result
+
+ def __str__(self):
+ return self._mangle(self._content)
+
+ def __eq__(self, other):
+ if not np.allclose(self.affine, other.affine):
+ return False
+ if self.shape != other.shape:
+ return False
+ if self.slice_dim != other.slice_dim:
+ return False
+ if self.version != other.version:
+ return False
+ for classes in self.get_valid_classes():
+ if (dict(self.get_class_dict(classes)) !=
+ dict(other.get_class_dict(classes))):
+ return False
+
+ return True
+
+ def _unmangle(self, value):
+ '''Go from extension data to runtime representation.'''
+ #Its not possible to preserve order while loading with python 2.6
+ kwargs = {}
+ if sys.version_info >= (2, 7):
+ kwargs['object_pairs_hook'] = OrderedDict
+ return json.loads(value, **kwargs)
+
+ def _mangle(self, value):
+ '''Go from runtime representation to extension data.'''
+ return json.dumps(value, indent=4)
+
+ _const_tests = {('global', 'slices') : (('global', 'const'),
+ ('vector', 'samples'),
+ ('time', 'samples')
+ ),
+ ('vector', 'slices') : (('global', 'const'),
+ ('time', 'samples')
+ ),
+ ('time', 'slices') : (('global', 'const'),
+ ),
+ ('time', 'samples') : (('global', 'const'),
+ ('vector', 'samples'),
+ ),
+ ('vector', 'samples') : (('global', 'const'),)
+ }
+ '''Classification mapping showing possible reductions in multiplicity for
+ values that are constant with some period.'''
+
+ def _get_const_period(self, src_cls, dest_cls):
+ '''Get the period over which we test for const-ness with for the
+ given classification change.'''
+ if dest_cls == ('global', 'const'):
+ return None
+ elif src_cls == ('global', 'slices'):
+ return self.get_multiplicity(src_cls) / self.get_multiplicity(dest_cls)
+ elif src_cls == ('vector', 'slices'): #implies dest_cls == ('time', 'samples'):
+ return self.n_slices
+ elif src_cls == ('time', 'samples'): #implies dest_cls == ('vector', 'samples')
+ return self.shape[3]
+ assert False #Should take one of the above branches
+
+ _repeat_tests = {('global', 'slices') : (('time', 'slices'),
+ ('vector', 'slices')
+ ),
+ ('vector', 'slices') : (('time', 'slices'),),
+ }
+ '''Classification mapping showing possible reductions in multiplicity for
+ values that are repeating with some period.'''
+
+ def _simplify(self, key):
+ '''Try to simplify (reduce the multiplicity) of a single meta data
+ element by changing its classification. Return True if the
+ classification is changed, otherwise False.
+
+ Looks for values that are constant or repeating with some pattern.
+ Constant elements with a value of None will be deleted.
+ '''
+ values, curr_class = self.get_values_and_class(key)
+
+ #If the class is global const then just delete it if the value is None
+ if curr_class == ('global', 'const'):
+ if values is None:
+ del self.get_class_dict(curr_class)[key]
+ return True
+ return False
+
+ #Test if the values are constant with some period
+ dests = self._const_tests[curr_class]
+ for dest_cls in dests:
+ if dest_cls[0] in self._content:
+ period = self._get_const_period(curr_class, dest_cls)
+ #If the period is one, the two classifications have the
+ #same multiplicity so we are dealing with a degenerate
+ #case (i.e. single slice data). Just change the
+ #classification to the "simpler" one in this case
+ if period == 1 or is_constant(values, period):
+ if period is None:
+ self.get_class_dict(dest_cls)[key] = \
+ values[0]
+ else:
+ self.get_class_dict(dest_cls)[key] = \
+ values[::period]
+ break
+ else: #Otherwise test if values are repeating with some period
+ if curr_class in self._repeat_tests:
+ for dest_cls in self._repeat_tests[curr_class]:
+ if dest_cls[0] in self._content:
+ dest_mult = self.get_multiplicity(dest_cls)
+ if is_repeating(values, dest_mult):
+ self.get_class_dict(dest_cls)[key] = \
+ values[:dest_mult]
+ break
+ else: #Can't simplify
+ return False
+ else:
+ return False
+
+ del self.get_class_dict(curr_class)[key]
+ return True
+
+ _preserving_changes = {None : (('global', 'const'),
+ ('vector', 'samples'),
+ ('time', 'samples'),
+ ('time', 'slices'),
+ ('vector', 'slices'),
+ ('global', 'slices'),
+ ),
+ ('global', 'const') : (('vector', 'samples'),
+ ('time', 'samples'),
+ ('time', 'slices'),
+ ('vector', 'slices'),
+ ('global', 'slices'),
+ ),
+ ('vector', 'samples') : (('time', 'samples'),
+ ('global', 'slices'),
+ ),
+ ('time', 'samples') : (('global', 'slices'),
+ ),
+ ('time', 'slices') : (('vector', 'slices'),
+ ('global', 'slices'),
+ ),
+ ('vector', 'slices') : (('global', 'slices'),
+ ),
+ ('global', 'slices') : tuple(),
+ }
+ '''Classification mapping showing allowed changes when increasing the
+ multiplicity.'''
+
+ def _get_changed_class(self, key, new_class, slice_dim=None):
+ '''Get an array of values corresponding to a single meta data
+ element with its classification changed by increasing its
+ multiplicity. This will preserve all the meta data and allow easier
+ merging of values with different classifications.'''
+ values, curr_class = self.get_values_and_class(key)
+ if curr_class == new_class:
+ return values
+
+ if not new_class in self._preserving_changes[curr_class]:
+ raise ValueError("Classification change would lose data.")
+
+ if curr_class is None:
+ curr_mult = 1
+ per_slice = False
+ else:
+ curr_mult = self.get_multiplicity(curr_class)
+ per_slice = curr_class[1] == 'slices'
+ if new_class in self.get_valid_classes():
+ new_mult = self.get_multiplicity(new_class)
+ #Only way we get 0 for mult is if slice dim is undefined
+ if new_mult == 0:
+ new_mult = self.shape[slice_dim]
+ else:
+ new_mult = 1
+ mult_fact = new_mult / curr_mult
+ if curr_mult == 1:
+ values = [values]
+
+
+ if per_slice:
+ result = values * mult_fact
+ else:
+ result = []
+ for value in values:
+ result.extend([deepcopy(value)] * mult_fact)
+
+ if new_class == ('global', 'const'):
+ result = result[0]
+
+ return result
+
+
+ def _change_class(self, key, new_class):
+ '''Change the classification of the meta data element in place. See
+ _get_changed_class.'''
+ values, curr_class = self.get_values_and_class(key)
+ if curr_class == new_class:
+ return
+
+ self.get_class_dict(new_class)[key] = self._get_changed_class(key,
+ new_class)
+
+ if not curr_class is None:
+ del self.get_class_dict(curr_class)[key]
+
+
+
+ def _copy_slice(self, other, src_class, idx):
+ '''Get a copy of the meta data from the 'other' instance with
+ classification 'src_class', corresponding to the slice with index
+ 'idx'.'''
+ if src_class[0] == 'global':
+ for classes in (('time', 'samples'),
+ ('vector', 'samples'),
+ ('global', 'const')):
+ if classes in self.get_valid_classes():
+ dest_class = classes
+ break
+ elif src_class[0] == 'vector':
+ for classes in (('time', 'samples'),
+ ('global', 'const')):
+ if classes in self.get_valid_classes():
+ dest_class = classes
+ break
+ else:
+ dest_class = ('global', 'const')
+
+ src_dict = other.get_class_dict(src_class)
+ dest_dict = self.get_class_dict(dest_class)
+ dest_mult = self.get_multiplicity(dest_class)
+ stride = other.n_slices
+ for key, vals in src_dict.iteritems():
+ subset_vals = vals[idx::stride]
+
+ if len(subset_vals) < dest_mult:
+ full_vals = []
+ for val_idx in xrange(dest_mult / len(subset_vals)):
+ full_vals += deepcopy(subset_vals)
+ subset_vals = full_vals
+ if len(subset_vals) == 1:
+ subset_vals = subset_vals[0]
+ dest_dict[key] = deepcopy(subset_vals)
+ self._simplify(key)
+
+ def _global_slice_subset(self, key, sample_base, idx):
+ '''Get a subset of the meta data values with the classificaion
+ ('global', 'slices') corresponding to a single sample along the
+ time or vector dimension (as specified by 'sample_base' and 'idx').
+ '''
+ n_slices = self.n_slices
+ shape = self.shape
+ src_dict = self.get_class_dict(('global', 'slices'))
+ if sample_base == 'vector':
+ slices_per_vec = n_slices * shape[3]
+ start_idx = idx * slices_per_vec
+ end_idx = start_idx + slices_per_vec
+ return src_dict[key][start_idx:end_idx]
+ else:
+ if not ('vector', 'samples') in self.get_valid_classes():
+ start_idx = idx * n_slices
+ end_idx = start_idx + n_slices
+ return src_dict[key][start_idx:end_idx]
+ else:
+ result = []
+ slices_per_vec = n_slices * shape[3]
+ for vec_idx in xrange(shape[4]):
+ start_idx = (vec_idx * slices_per_vec) + (idx * n_slices)
+ end_idx = start_idx + n_slices
+ result.extend(src_dict[key][start_idx:end_idx])
+ return result
+
+ def _copy_sample(self, other, src_class, sample_base, idx):
+ '''Get a copy of meta data from 'other' instance with classification
+ 'src_class', corresponding to one sample along the time or vector
+ dimension.'''
+ assert src_class != ('global', 'const')
+ src_dict = other.get_class_dict(src_class)
+ if src_class[1] == 'samples':
+ #If we are indexing on the same dim as the src_class we need to
+ #change the classification
+ if src_class[0] == sample_base:
+ #Time samples may become vector samples, otherwise const
+ best_dest = None
+ for dest_cls in (('vector', 'samples'),
+ ('global', 'const')):
+ if (dest_cls != src_class and
+ dest_cls in self.get_valid_classes()
+ ):
+ best_dest = dest_cls
+ break
+
+ dest_mult = self.get_multiplicity(dest_cls)
+ if dest_mult == 1:
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(dest_cls)[key] = \
+ deepcopy(vals[idx])
+ else: #We must be doing time samples -> vector samples
+ stride = other.shape[3]
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(dest_cls)[key] = \
+ deepcopy(vals[idx::stride])
+ for key in src_dict.keys():
+ self._simplify(key)
+
+ else: #Otherwise classification does not change
+ #The multiplicity will change for time samples if splitting
+ #vector dimension
+ if src_class == ('time', 'samples'):
+ dest_mult = self.get_multiplicity(src_class)
+ start_idx = idx * dest_mult
+ end_idx = start_idx + dest_mult
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(src_class)[key] = \
+ deepcopy(vals[start_idx:end_idx])
+ self._simplify(key)
+ else: #Otherwise multiplicity is unchanged
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(src_class)[key] = deepcopy(vals)
+ else: #The src_class is per slice
+ if src_class[0] == sample_base:
+ best_dest = None
+ for dest_class in self._preserving_changes[src_class]:
+ if dest_class in self.get_valid_classes():
+ best_dest = dest_class
+ break
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(best_dest)[key] = deepcopy(vals)
+ elif src_class[0] != 'global':
+ if sample_base == 'time':
+ #Take a subset of vector slices
+ n_slices = self.n_slices
+ start_idx = idx * n_slices
+ end_idx = start_idx + n_slices
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(src_class)[key] = \
+ deepcopy(vals[start_idx:end_idx])
+ self._simplify(key)
+ else:
+ #Time slices are unchanged
+ for key, vals in src_dict.iteritems():
+ self.get_class_dict(src_class)[key] = deepcopy(vals)
+ else:
+ #Take a subset of global slices
+ for key, vals in src_dict.iteritems():
+ subset_vals = \
+ other._global_slice_subset(key, sample_base, idx)
+ self.get_class_dict(src_class)[key] = deepcopy(subset_vals)
+ self._simplify(key)
+
+ def _insert(self, dim, other):
+ self_slc_norm = self.slice_normal
+ other_slc_norm = other.slice_normal
+
+ #If we are not using slice meta data, temporarily remove it from the
+ #other dcmmeta object
+ use_slices = (not self_slc_norm is None and
+ not other_slc_norm is None and
+ np.allclose(self_slc_norm, other_slc_norm))
+ other_slc_meta = {}
+ if not use_slices:
+ for classes in other.get_valid_classes():
+ if classes[1] == 'slices':
+ other_slc_meta[classes] = other.get_class_dict(classes)
+ other._content[classes[0]][classes[1]] = {}
+ missing_keys = list(set(self.get_keys()) - set(other.get_keys()))
+ for other_classes in other.get_valid_classes():
+ other_keys = other.get_class_dict(other_classes).keys()
+
+ #Treat missing keys as if they were in global const and have a value
+ #of None
+ if other_classes == ('global', 'const'):
+ other_keys += missing_keys
+
+ #When possible, reclassify our meta data so it matches the other
+ #classification
+ for key in other_keys:
+ local_classes = self.get_classification(key)
+ if local_classes != other_classes:
+ local_allow = self._preserving_changes[local_classes]
+ other_allow = self._preserving_changes[other_classes]
+
+ if other_classes in local_allow:
+ self._change_class(key, other_classes)
+ elif not local_classes in other_allow:
+ best_dest = None
+ for dest_class in local_allow:
+ if (dest_class[0] in self._content and
+ dest_class in other_allow):
+ best_dest = dest_class
+ break
+ self._change_class(key, best_dest)
+
+ #Insert new meta data and further reclassify as necessary
+ for key in other_keys:
+ if dim == self.slice_dim:
+ self._insert_slice(key, other)
+ elif dim < 3:
+ self._insert_non_slice(key, other)
+ elif dim == 3:
+ self._insert_sample(key, other, 'time')
+ elif dim == 4:
+ self._insert_sample(key, other, 'vector')
+
+ #Restore per slice meta if needed
+ if not use_slices:
+ for classes in other.get_valid_classes():
+ if classes[1] == 'slices':
+ other._content[classes[0]][classes[1]] = \
+ other_slc_meta[classes]
+
+ def _insert_slice(self, key, other):
+ local_vals, classes = self.get_values_and_class(key)
+ other_vals = other._get_changed_class(key, classes, self.slice_dim)
+
+
+ #Handle some common / simple insertions with special cases
+ if classes == ('global', 'const'):
+ if local_vals != other_vals:
+ for dest_base in ('time', 'vector', 'global'):
+ if dest_base in self._content:
+ self._change_class(key, (dest_base, 'slices'))
+ other_vals = other._get_changed_class(key,
+ (dest_base,
+ 'slices'),
+ self.slice_dim
+ )
+ self.get_values(key).extend(other_vals)
+ break
+ elif classes == ('time', 'slices'):
+ local_vals.extend(other_vals)
+ else:
+ #Default to putting in global slices and simplifying later
+ if classes != ('global', 'slices'):
+ self._change_class(key, ('global', 'slices'))
+ local_vals = self.get_class_dict(('global', 'slices'))[key]
+ other_vals = other._get_changed_class(key,
+ ('global', 'slices'),
+ self.slice_dim)
+
+ #Need to interleave slices from different volumes
+ n_slices = self.n_slices
+ other_n_slices = other.n_slices
+ shape = self.shape
+ n_vols = 1
+ for dim_size in shape[3:]:
+ n_vols *= dim_size
+
+ intlv = []
+ loc_start = 0
+ oth_start = 0
+ for vol_idx in xrange(n_vols):
+ intlv += local_vals[loc_start:loc_start + n_slices]
+ intlv += other_vals[oth_start:oth_start + other_n_slices]
+ loc_start += n_slices
+ oth_start += other_n_slices
+
+ self.get_class_dict(('global', 'slices'))[key] = intlv
+
+ def _insert_non_slice(self, key, other):
+ local_vals, classes = self.get_values_and_class(key)
+ other_vals = other._get_changed_class(key, classes, self.slice_dim)
+
+ if local_vals != other_vals:
+ del self.get_class_dict(classes)[key]
+
+ def _insert_sample(self, key, other, sample_base):
+ local_vals, classes = self.get_values_and_class(key)
+ other_vals = other._get_changed_class(key, classes, self.slice_dim)
+
+ if classes == ('global', 'const'):
+ if local_vals != other_vals:
+ self._change_class(key, (sample_base, 'samples'))
+ local_vals = self.get_values(key)
+ other_vals = other._get_changed_class(key,
+ (sample_base, 'samples'),
+ self.slice_dim
+ )
+ local_vals.extend(other_vals)
+ elif classes == (sample_base, 'samples'):
+ local_vals.extend(other_vals)
+ else:
+ if classes != ('global', 'slices'):
+ self._change_class(key, ('global', 'slices'))
+ local_vals = self.get_values(key)
+ other_vals = other._get_changed_class(key,
+ ('global', 'slices'),
+ self.slice_dim)
+
+ shape = self.shape
+ n_dims = len(shape)
+ if sample_base == 'time' and n_dims == 5:
+ #Need to interleave values from the time points in each vector
+ #component
+ n_slices = self.n_slices
+ slices_per_vec = n_slices * shape[3]
+ oth_slc_per_vec = n_slices * other.shape[3]
+
+ intlv = []
+ loc_start = 0
+ oth_start = 0
+ for vec_idx in xrange(shape[4]):
+ intlv += local_vals[loc_start:loc_start+slices_per_vec]
+ intlv += other_vals[oth_start:oth_start+oth_slc_per_vec]
+ loc_start += slices_per_vec
+ oth_start += oth_slc_per_vec
+
+ self.get_class_dict(('global', 'slices'))[key] = intlv
+ else:
+ local_vals.extend(other_vals)
+
+#Add our extension to nibabel
+nb.nifti1.extension_codes.add_codes(((dcm_meta_ecode,
+ "dcmmeta",
+ DcmMetaExtension),)
+ )
+
+class MissingExtensionError(Exception):
+ '''Exception denoting that there is no DcmMetaExtension in the Nifti header.
+ '''
+ def __str__(self):
+ return 'No dcmmeta extension found.'
+
+def patch_dcm_ds_is(dcm):
+ '''Convert all elements in `dcm` with VR of 'DS' or 'IS' to floats and ints.
+ This is a hackish work around for the backwards incompatibility of pydicom
+ 0.9.7 and should not be needed once nibabel is updated.
+ '''
+ for elem in dcm:
+ if elem.VM == 1:
+ if elem.VR in ('DS', 'IS'):
+ if elem.value == '':
+ continue
+ if elem.VR == 'DS':
+ elem.VR = 'FD'
+ elem.value = float(elem.value)
+ else:
+ elem.VR = 'SL'
+ elem.value = int(elem.value)
+ else:
+ if elem.VR in ('DS', 'IS'):
+ if elem.value == '':
+ continue
+ if elem.VR == 'DS':
+ elem.VR = 'FD'
+ elem.value = [float(val) for val in elem.value]
+ else:
+ elem.VR = 'SL'
+ elem.value = [int(val) for val in elem.value]
+
+
+class NiftiWrapper(object):
+ '''Wraps a Nifti1Image object containing a DcmMeta header extension.
+ Provides access to the meta data and the ability to split or merge the
+ data array while updating the meta data.
+
+ Parameters
+ ----------
+ nii_img : nibabel.nifti1.Nifti1Image
+ The Nifti1Image to wrap.
+
+ make_empty : bool
+ If True an empty DcmMetaExtension will be created if none is found.
+
+ Raises
+ ------
+ MissingExtensionError
+ No valid DcmMetaExtension was found.
+
+ ValueError
+ More than one valid DcmMetaExtension was found.
+ '''
+
+ def __init__(self, nii_img, make_empty=False):
+ self.nii_img = nii_img
+ hdr = nii_img.get_header()
+ self.meta_ext = None
+ for extension in hdr.extensions:
+ if extension.get_code() == dcm_meta_ecode:
+ try:
+ extension.check_valid()
+ except InvalidExtensionError as e:
+ print("Found candidate extension, but invalid: %s" % e)
+ else:
+ if not self.meta_ext is None:
+ raise ValueError('More than one valid DcmMeta '
+ 'extension found.')
+ self.meta_ext = extension
+ if not self.meta_ext:
+ if make_empty:
+ slice_dim = hdr.get_dim_info()[2]
+ self.meta_ext = \
+ DcmMetaExtension.make_empty(self.nii_img.shape,
+ hdr.get_best_affine(),
+ None,
+ slice_dim)
+ hdr.extensions.append(self.meta_ext)
+ else:
+ raise MissingExtensionError
+ self.meta_ext.check_valid()
+
+ def __getitem__(self, key):
+ '''Get the value for the given meta data key. Only considers meta data
+ that is globally constant. To access varying meta data you must use the
+ method 'get_meta'.'''
+ return self.meta_ext.get_class_dict(('global', 'const'))[key]
+
+ def meta_valid(self, classification):
+ '''Return true if the meta data with the given classification appears
+ to be valid for the wrapped Nifti image. Considers the shape and
+ orientation of the image and the meta data extension.'''
+ if classification == ('global', 'const'):
+ return True
+
+ img_shape = self.nii_img.get_shape()
+ meta_shape = self.meta_ext.shape
+ if classification == ('vector', 'samples'):
+ return meta_shape[4:] == img_shape[4:]
+ if classification == ('time', 'samples'):
+ return meta_shape[3:] == img_shape[3:]
+
+ hdr = self.nii_img.get_header()
+ if self.meta_ext.n_slices != hdr.get_n_slices():
+ return False
+
+ slice_dim = hdr.get_dim_info()[2]
+ slice_dir = self.nii_img.get_affine()[slice_dim, :3]
+ slices_aligned = np.allclose(slice_dir,
+ self.meta_ext.slice_normal,
+ atol=1e-6)
+
+ if classification == ('time', 'slices'):
+ return slices_aligned
+ if classification == ('vector', 'slices'):
+ return meta_shape[3] == img_shape[3] and slices_aligned
+ if classification == ('global', 'slices'):
+ return meta_shape[3:] == img_shape[3:] and slices_aligned
+
+ def get_meta(self, key, index=None, default=None):
+ '''Return the meta data value for the provided `key`.
+
+ Parameters
+ ----------
+ key : str
+ The meta data key.
+
+ index : tuple
+ The voxel index we are interested in.
+
+ default
+ This will be returned if the meta data for `key` is not found.
+
+ Returns
+ -------
+ value
+ The meta data value for the given `key` (and optionally `index`)
+
+ Notes
+ -----
+ The per-sample and per-slice meta data will only be considered if the
+ `samples_valid` and `slices_valid` methods return True (respectively),
+ and an `index` is specified.
+ '''
+ #Get the value(s) and classification for the key
+ values, classes = self.meta_ext.get_values_and_class(key)
+ if classes is None:
+ return default
+
+ #Check if the value is constant
+ if classes == ('global', 'const'):
+ return values
+
+ #Check if the classification is valid
+ if not self.meta_valid(classes):
+ return default
+
+ #If an index is provided check the varying values
+ if not index is None:
+ #Test if the index is valid
+ shape = self.nii_img.get_shape()
+ if len(index) != len(shape):
+ raise IndexError('Incorrect number of indices.')
+ for dim, ind_val in enumerate(index):
+ if not 0 <= ind_val < shape[dim]:
+ raise IndexError('Index is out of bounds.')
+
+ #First try per time/vector sample values
+ if classes == ('time', 'samples'):
+ return values[index[3]]
+ if classes == ('vector', 'samples'):
+ return values[index[4]]
+
+ #Finally, if aligned, try per-slice values
+ slice_dim = self.nii_img.get_header().get_dim_info()[2]
+ n_slices = shape[slice_dim]
+ if classes == ('global', 'slices'):
+ val_idx = index[slice_dim]
+ for count, idx_val in enumerate(index[3:]):
+ val_idx += idx_val * n_slices
+ n_slices *= shape[count+3]
+ return values[val_idx]
+ elif classes == ('time', 'slices'):
+ val_idx = index[slice_dim]
+ return values[val_idx]
+ elif classes == ('vector', 'slices'):
+ val_idx = index[slice_dim]
+ val_idx += index[3]*n_slices
+ return values[val_idx]
+
+ return default
+
+ def remove_extension(self):
+ '''Remove the DcmMetaExtension from the header of nii_img. The
+ attribute `meta_ext` will still point to the extension.'''
+ hdr = self.nii_img.get_header()
+ target_idx = None
+ for idx, ext in enumerate(hdr.extensions):
+ if id(ext) == id(self.meta_ext):
+ target_idx = idx
+ break
+ else:
+ raise IndexError('Extension not found in header')
+ del hdr.extensions[target_idx]
+ # Nifti1Image.update_header will increase this if necessary
+ hdr['vox_offset'] = 0
+
+ def replace_extension(self, dcmmeta_ext):
+ '''Replace the DcmMetaExtension.
+
+ Parameters
+ ----------
+ dcmmeta_ext : DcmMetaExtension
+ The new DcmMetaExtension.
+
+ '''
+ self.remove_extension()
+ self.nii_img.get_header().extensions.append(dcmmeta_ext)
+ self.meta_ext = dcmmeta_ext
+
+ def split(self, dim=None):
+ '''Generate splits of the array and meta data along the specified
+ dimension.
+
+ Parameters
+ ----------
+ dim : int
+ The dimension to split the voxel array along. If None it will
+ prefer the vector, then time, then slice dimensions.
+
+ Returns
+ -------
+ result
+ Generator which yields a NiftiWrapper result for each index
+ along `dim`.
+
+ '''
+ shape = self.nii_img.get_shape()
+ data = self.nii_img.get_data()
+ header = self.nii_img.get_header()
+ slice_dim = header.get_dim_info()[2]
+
+ #If dim is None, choose the vector/time/slice dim in that order
+ if dim is None:
+ dim = len(shape) - 1
+ if dim == 2:
+ if slice_dim is None:
+ raise ValueError("Slice dimension is not known")
+ dim = slice_dim
+
+ #If we are splitting on a spatial dimension, we need to update the
+ #translation
+ trans_update = None
+ if dim < 3:
+ trans_update = header.get_best_affine()[:3, dim]
+
+ split_hdr = header.copy()
+ slices = [slice(None)] * len(shape)
+ for idx in xrange(shape[dim]):
+ #Grab the split data, get rid of trailing singular dimensions
+ if dim >= 3 and dim == len(shape) - 1:
+ slices[dim] = idx
+ else:
+ slices[dim] = slice(idx, idx+1)
+
+ split_data = data[slices].copy()
+
+ #Update the translation in any affines if needed
+ if not trans_update is None and idx != 0:
+ qform = split_hdr.get_qform()
+ if not qform is None:
+ qform[:3, 3] += trans_update
+ split_hdr.set_qform(qform)
+ sform = split_hdr.get_sform()
+ if not sform is None:
+ sform[:3, 3] += trans_update
+ split_hdr.set_sform(sform)
+
+ #Create the initial Nifti1Image object
+ split_nii = nb.Nifti1Image(split_data,
+ split_hdr.get_best_affine(),
+ header=split_hdr)
+
+ #Replace the meta data with the appropriate subset
+ meta_dim = dim
+ if dim == slice_dim:
+ meta_dim = self.meta_ext.slice_dim
+ split_meta = self.meta_ext.get_subset(meta_dim, idx)
+ result = NiftiWrapper(split_nii)
+ result.replace_extension(split_meta)
+
+ yield result
+
+ def to_filename(self, out_path):
+ '''Write out the wrapped Nifti to a file
+
+ Parameters
+ ----------
+ out_path : str
+ The path to write out the file to
+
+ Notes
+ -----
+ Will check that the DcmMetaExtension is valid before writing the file.
+ '''
+ self.meta_ext.check_valid()
+ self.nii_img.to_filename(out_path)
+
+ @classmethod
+ def from_filename(klass, path):
+ '''Create a NiftiWrapper from a file.
+
+ Parameters
+ ----------
+ path : str
+ The path to the Nifti file to load.
+ '''
+ return klass(nb.load(path))
+
+ @classmethod
+ def from_dicom_wrapper(klass, dcm_wrp, meta_dict=None):
+ '''Create a NiftiWrapper from a nibabel DicomWrapper.
+
+ Parameters
+ ----------
+ dcm_wrap : nicom.dicomwrappers.DicomWrapper
+ The dataset to convert into a NiftiWrapper.
+
+ meta_dict : dict
+ An optional dictionary of meta data extracted from `dcm_data`. See
+ the `extract` module for generating this dict.
+
+ '''
+ data = dcm_wrp.get_data()
+
+ #The Nifti patient space flips the x and y directions
+ affine = np.dot(np.diag([-1., -1., 1., 1.]), dcm_wrp.get_affine())
+
+ #Make 2D data 3D
+ if len(data.shape) == 2:
+ data = data.reshape(data.shape + (1,))
+
+ #Create the nifti image and set header data
+ nii_img = nb.nifti1.Nifti1Image(data, affine)
+ hdr = nii_img.get_header()
+ hdr.set_xyzt_units('mm', 'sec')
+ dim_info = {'freq' : None,
+ 'phase' : None,
+ 'slice' : 2
+ }
+ if hasattr(dcm_wrp.dcm_data, 'InplanePhaseEncodingDirection'):
+ if dcm_wrp['InplanePhaseEncodingDirection'] == 'ROW':
+ dim_info['phase'] = 1
+ dim_info['freq'] = 0
+ else:
+ dim_info['phase'] = 0
+ dim_info['freq'] = 1
+ hdr.set_dim_info(**dim_info)
+
+ #Embed the meta data extension
+ result = klass(nii_img, make_empty=True)
+
+ result.meta_ext.reorient_transform = np.eye(4)
+ if meta_dict:
+ result.meta_ext.get_class_dict(('global', 'const')).update(meta_dict)
+
+ return result
+
+ @classmethod
+ def from_dicom(klass, dcm_data, meta_dict=None):
+ '''Create a NiftiWrapper from a single DICOM dataset.
+
+ Parameters
+ ----------
+ dcm_data : dicom.dataset.Dataset
+ The DICOM dataset to convert into a NiftiWrapper.
+
+ meta_dict : dict
+ An optional dictionary of meta data extracted from `dcm_data`. See
+ the `extract` module for generating this dict.
+
+ '''
+ dcm_wrp = wrapper_from_data(dcm_data)
+ return klass.from_dicom_wrapper(dcm_wrp, meta_dict)
+
+ @classmethod
+ def from_sequence(klass, seq, dim=None):
+ '''Create a NiftiWrapper by joining a sequence of NiftiWrapper objects
+ along the given dimension.
+
+ Parameters
+ ----------
+ seq : sequence
+ The sequence of NiftiWrapper objects.
+
+ dim : int
+ The dimension to join the NiftiWrapper objects along. If None,
+ 2D inputs will become 3D, 3D inputs will become 4D, and 4D inputs
+ will become 5D.
+
+ Returns
+ -------
+ result : NiftiWrapper
+ The merged NiftiWrapper with updated meta data.
+ '''
+ n_inputs = len(seq)
+ first_input = seq[0]
+ first_nii = first_input.nii_img
+ first_hdr = first_nii.get_header()
+ shape = first_nii.shape
+ affine = first_nii.get_affine().copy()
+
+ #If dim is None, choose a sane default
+ if dim is None:
+ if len(shape) == 3:
+ singular_dim = None
+ for dim_idx, dim_size in enumerate(shape):
+ if dim_size == 1:
+ singular_dim = dim_idx
+ if singular_dim is None:
+ dim = 3
+ else:
+ dim = singular_dim
+ if len(shape) == 4:
+ dim = 4
+ else:
+ if not 0 <= dim < 5:
+ raise ValueError("The argument 'dim' must be in the range "
+ "[0, 5).")
+ if dim < len(shape) and shape[dim] != 1:
+ raise ValueError('The dimension must be singular or not exist')
+
+ #Pull out the three axes vectors for validation of other input affines
+ axes = []
+ for axis_idx in xrange(3):
+ axis_vec = affine[:3, axis_idx]
+ if axis_idx == dim:
+ axis_vec = axis_vec.copy()
+ axis_vec /= np.sqrt(np.dot(axis_vec, axis_vec))
+ axes.append(axis_vec)
+ #Pull out the translation
+ trans = affine[:3, 3]
+
+ #Determine the shape of the result data array and create it
+ result_shape = list(shape)
+ while dim >= len(result_shape):
+ result_shape.append(1)
+ result_shape[dim] = n_inputs
+
+ result_dtype = max(input_wrp.nii_img.get_data().dtype
+ for input_wrp in seq)
+ result_data = np.empty(result_shape, dtype=result_dtype)
+
+ #Start with the header info from the first input
+ hdr_info = {'qform' : first_hdr.get_qform(),
+ 'qform_code' : first_hdr['qform_code'],
+ 'sform' : first_hdr.get_sform(),
+ 'sform_code' : first_hdr['sform_code'],
+ 'dim_info' : list(first_hdr.get_dim_info()),
+ 'xyzt_units' : list(first_hdr.get_xyzt_units()),
+ }
+
+ try:
+ hdr_info['slice_duration'] = first_hdr.get_slice_duration()
+ except HeaderDataError:
+ hdr_info['slice_duration'] = None
+ try:
+ hdr_info['intent'] = first_hdr.get_intent()
+ except HeaderDataError:
+ hdr_info['intent'] = None
+ try:
+ hdr_info['slice_times'] = first_hdr.get_slice_times()
+ except HeaderDataError:
+ hdr_info['slice_times'] = None
+
+ #Fill the data array, check header consistency
+ data_slices = [slice(None)] * len(result_shape)
+ for dim_idx, dim_size in enumerate(result_shape):
+ if dim_size == 1:
+ data_slices[dim_idx] = 0
+ last_trans = None #Keep track of the translation from last input
+ for input_idx in range(n_inputs):
+
+ input_wrp = seq[input_idx]
+ input_nii = input_wrp.nii_img
+ input_aff = input_nii.get_affine()
+ input_hdr = input_nii.get_header()
+
+ #Check that the affines match appropriately
+ for axis_idx, axis_vec in enumerate(axes):
+ in_vec = input_aff[:3, axis_idx]
+
+ #If we are joining on this dimension
+ if axis_idx == dim:
+ #Allow scaling difference as it will be updated later
+ in_vec = in_vec.copy()
+ in_vec /= np.sqrt(np.dot(in_vec, in_vec))
+
+ in_trans = input_aff[:3, 3]
+ if not last_trans is None:
+ #Must be translated along the axis
+ trans_diff = in_trans - last_trans
+ if not np.allclose(trans_diff, 0.0):
+ trans_diff /= np.sqrt(np.dot(trans_diff, trans_diff))
+
+ if (np.allclose(trans_diff, 0.0) or
+ not np.allclose(np.dot(trans_diff, in_vec),
+ 1.0,
+ atol=1e-6)
+ ):
+ raise ValueError("Slices must be translated along the "
+ "normal direction")
+
+ #Update reference to last translation
+ last_trans = in_trans
+
+ #Check that axis vectors match
+ if not np.allclose(in_vec, axis_vec, atol=5e-4):
+ raise ValueError("Cannot join images with different "
+ "orientations.")
+
+
+ data_slices[dim] = input_idx
+ result_data[data_slices] = input_nii.get_data().squeeze()
+
+ if input_idx != 0:
+ if (hdr_info['qform'] is None or
+ input_hdr.get_qform() is None or
+ not np.allclose(input_hdr.get_qform(), hdr_info['qform'])
+ ):
+ hdr_info['qform'] = None
+ if input_hdr['qform_code'] != hdr_info['qform_code']:
+ hdr_info['qform_code'] = None
+ if (hdr_info['sform'] is None or
+ input_hdr.get_sform() is None or
+ not np.allclose(input_hdr.get_sform(), hdr_info['sform'])
+ ):
+ hdr_info['sform'] = None
+ if input_hdr['sform_code'] != hdr_info['sform_code']:
+ hdr_info['sform_code'] = None
+ in_dim_info = list(input_hdr.get_dim_info())
+ if in_dim_info != hdr_info['dim_info']:
+ for idx in xrange(3):
+ if in_dim_info[idx] != hdr_info['dim_info'][idx]:
+ hdr_info['dim_info'][idx] = None
+ in_xyzt_units = list(input_hdr.get_xyzt_units())
+ if in_xyzt_units != hdr_info['xyzt_units']:
+ for idx in xrange(2):
+ if in_xyzt_units[idx] != hdr_info['xyzt_units'][idx]:
+ hdr_info['xyzt_units'][idx] = None
+
+ try:
+ if input_hdr.get_slice_duration() != hdr_info['slice_duration']:
+ hdr_info['slice_duration'] = None
+ except HeaderDataError:
+ hdr_info['slice_duration'] = None
+ try:
+ if input_hdr.get_intent() != hdr_info['intent']:
+ hdr_info['intent'] = None
+ except HeaderDataError:
+ hdr_info['intent'] = None
+ try:
+ if input_hdr.get_slice_times() != hdr_info['slice_times']:
+ hdr_info['slice_times'] = None
+ except HeaderDataError:
+ hdr_info['slice_times'] = None
+
+ #If we joined along a spatial dim, rescale the appropriate axis
+ scaled_dim_dir = None
+ if dim < 3:
+ scaled_dim_dir = seq[1].nii_img.get_affine()[:3, 3] - trans
+ affine[:3, dim] = scaled_dim_dir
+
+ #Create the resulting Nifti and wrapper
+ result_nii = nb.Nifti1Image(result_data, affine)
+ result_hdr = result_nii.get_header()
+
+ #Update the header with any info that is consistent across inputs
+ if hdr_info['qform'] is not None and hdr_info['qform_code'] is not None:
+ if not scaled_dim_dir is None:
+ hdr_info['qform'][:3, dim] = scaled_dim_dir
+ result_nii.set_qform(hdr_info['qform'],
+ int(hdr_info['qform_code']),
+ update_affine=True)
+ if hdr_info['sform'] is not None and hdr_info['sform_code'] is not None:
+ if not scaled_dim_dir is None:
+ hdr_info['sform'][:3, dim] = scaled_dim_dir
+ result_nii.set_sform(hdr_info['sform'],
+ int(hdr_info['sform_code']),
+ update_affine=True)
+ if hdr_info['dim_info'] is not None:
+ result_hdr.set_dim_info(*hdr_info['dim_info'])
+ slice_dim = hdr_info['dim_info'][2]
+ else:
+ slice_dim = None
+ if hdr_info['intent'] is not None:
+ result_hdr.set_intent(*hdr_info['intent'])
+ if hdr_info['xyzt_units'] is not None:
+ result_hdr.set_xyzt_units(*hdr_info['xyzt_units'])
+ if hdr_info['slice_duration'] is not None:
+ result_hdr.set_slice_duration(hdr_info['slice_duration'])
+ if hdr_info['slice_times'] is not None:
+ result_hdr.set_slice_times(hdr_info['slice_times'])
+
+ #Create the meta data extension and insert it
+ seq_exts = [elem.meta_ext for elem in seq]
+ result_ext = DcmMetaExtension.from_sequence(seq_exts,
+ dim,
+ affine,
+ slice_dim)
+ result_hdr.extensions.append(result_ext)
+
+ return NiftiWrapper(result_nii)
diff --git a/fsl/data/dcmstack/dcmstack.py b/fsl/data/dcmstack/dcmstack.py
new file mode 100644
index 000000000..ee7d93f7c
--- /dev/null
+++ b/fsl/data/dcmstack/dcmstack.py
@@ -0,0 +1,1162 @@
+"""
+Stack DICOM datasets into volumes. The contents of this module are imported
+into the package namespace.
+"""
+import warnings, re, dicom
+from copy import deepcopy
+import nibabel as nb
+from nibabel.nifti1 import Nifti1Extensions
+from nibabel.spatialimages import HeaderDataError
+from nibabel.orientations import (io_orientation,
+ apply_orientation,
+ inv_ornt_aff)
+import numpy as np
+from .dcmmeta import DcmMetaExtension, NiftiWrapper
+
+with warnings.catch_warnings():
+ warnings.simplefilter('ignore')
+ from nibabel.nicom.dicomwrappers import wrapper_from_data
+
+def make_key_regex_filter(exclude_res, force_include_res=None):
+ '''Make a meta data filter using regular expressions.
+
+ Parameters
+ ----------
+ exclude_res : sequence
+ Sequence of regular expression strings. Any meta data where the key
+ matches one of these expressions will be excluded, unless it matches
+ one of the `force_include_res`.
+ force_include_res : sequence
+ Sequence of regular expression strings. Any meta data where the key
+ matches one of these expressions will be included.
+
+ Returns
+ -------
+ A callable which can be passed to `DicomStack` as the `meta_filter`.
+ '''
+ exclude_re = re.compile('|'.join(['(?:' + regex + ')'
+ for regex in exclude_res])
+ )
+ include_re = None
+ if force_include_res:
+ include_re = re.compile('|'.join(['(?:' + regex + ')'
+ for regex in force_include_res])
+ )
+
+ def key_regex_filter(key, value):
+ return (exclude_re.search(key) and
+ not (include_re and include_re.search(key)))
+ return key_regex_filter
+
+default_key_excl_res = ['Patient',
+ 'Physician',
+ 'Operator',
+ 'Date',
+ 'Birth',
+ 'Address',
+ 'Institution',
+ 'Station',
+ 'SiteName',
+ 'Age',
+ 'Comment',
+ 'Phone',
+ 'Telephone',
+ 'Insurance',
+ 'Religious',
+ 'Language',
+ 'Military',
+ 'MedicalRecord',
+ 'Ethnic',
+ 'Occupation',
+ 'Unknown',
+ 'PrivateTagData',
+ 'UID',
+ 'StudyDescription',
+ 'DeviceSerialNumber',
+ 'ReferencedImageSequence',
+ 'RequestedProcedureDescription',
+ 'PerformedProcedureStepDescription',
+ 'PerformedProcedureStepID',
+ ]
+'''A list of regexes passed to `make_key_regex_filter` as `exclude_res` to
+create the `default_meta_filter`.'''
+
+default_key_incl_res = ['ImageOrientationPatient',
+ 'ImagePositionPatient',
+ ]
+'''A list of regexes passed to `make_key_regex_filter` as `force_include_res`
+to create the `default_meta_filter`.'''
+
+default_meta_filter = make_key_regex_filter(default_key_excl_res,
+ default_key_incl_res)
+'''Default meta_filter for `DicomStack`.'''
+
+def ornt_transform(start_ornt, end_ornt):
+ '''Return the orientation that transforms from `start_ornt` to `end_ornt`.
+
+ Parameters
+ ----------
+ start_ornt : (n,2) orientation array
+ Initial orientation.
+
+ end_ornt : (n,2) orientation array
+ Final orientation.
+
+ Returns
+ -------
+ orientations : (p, 2) ndarray
+ The orientation that will transform the `start_ornt` to the `end_ornt`.
+ '''
+ start_ornt = np.asarray(start_ornt)
+ end_ornt = np.asarray(end_ornt)
+ if start_ornt.shape != end_ornt.shape:
+ raise ValueError("The orientations must have the same shape")
+ if start_ornt.shape[1] != 2:
+ raise ValueError("Invalid shape for an orientation: %s" %
+ start_ornt.shape)
+ result = np.empty_like(start_ornt)
+ for end_in_idx, (end_out_idx, end_flip) in enumerate(end_ornt):
+ for start_in_idx, (start_out_idx, start_flip) in enumerate(start_ornt):
+ if end_out_idx == start_out_idx:
+ if start_flip == end_flip:
+ flip = 1
+ else:
+ flip = -1
+ result[start_in_idx, :] = [end_in_idx, flip]
+ break
+ else:
+ raise ValueError("Unable to find out axis %d in start_ornt" %
+ end_out_idx)
+ return result
+
+def axcodes2ornt(axcodes, labels=None):
+ """ Convert axis codes `axcodes` to an orientation
+
+ Parameters
+ ----------
+ axcodes : (N,) tuple
+ axis codes - see ornt2axcodes docstring
+ labels : optional, None or sequence of (2,) sequences
+ (2,) sequences are labels for (beginning, end) of output axis. That
+ is, if the first element in `axcodes` is ``front``, and the second
+ (2,) sequence in `labels` is ('back', 'front') then the first
+ row of `ornt` will be ``[1, 1]``. If None, equivalent to
+ ``(('L','R'),('P','A'),('I','S'))`` - that is - RAS axes.
+
+ Returns
+ -------
+ ornt : (N,2) array-like
+ oritation array - see io_orientation docstring
+
+ Examples
+ --------
+ >>> axcodes2ornt(('F', 'L', 'U'), (('L','R'),('B','F'),('D','U')))
+ [[1, 1],[0,-1],[2,1]]
+ """
+
+ if labels is None:
+ labels = zip('LPI', 'RAS')
+
+ n_axes = len(axcodes)
+ ornt = np.ones((n_axes, 2), dtype=np.int8) * np.nan
+ for code_idx, code in enumerate(axcodes):
+ for label_idx, codes in enumerate(labels):
+ if code is None:
+ continue
+ if code in codes:
+ if code == codes[0]:
+ ornt[code_idx, :] = [label_idx, -1]
+ else:
+ ornt[code_idx, :] = [label_idx, 1]
+ break
+ return ornt
+
+def reorder_voxels(vox_array, affine, voxel_order):
+ '''Reorder the given voxel array and corresponding affine.
+
+ Parameters
+ ----------
+ vox_array : array
+ The array of voxel data
+
+ affine : array
+ The affine for mapping voxel indices to Nifti patient space
+
+ voxel_order : str
+ A three character code specifing the desired ending point for rows,
+ columns, and slices in terms of the orthogonal axes of patient space:
+ (l)eft, (r)ight, (a)nterior, (p)osterior, (s)uperior, and (i)nferior.
+
+ Returns
+ -------
+ out_vox : array
+ An updated view of vox_array.
+
+ out_aff : array
+ A new array with the updated affine
+
+ reorient_transform : array
+ The transform used to update the affine.
+
+ ornt_trans : tuple
+ The orientation transform used to update the orientation.
+
+ '''
+ #Check if voxel_order is valid
+ voxel_order = voxel_order.upper()
+ if len(voxel_order) != 3:
+ raise ValueError('The voxel_order must contain three characters')
+ dcm_axes = ['LR', 'AP', 'SI']
+ for char in voxel_order:
+ if not char in 'LRAPSI':
+ raise ValueError('The characters in voxel_order must be one '
+ 'of: L,R,A,P,I,S')
+ for idx, axis in enumerate(dcm_axes):
+ if char in axis:
+ del dcm_axes[idx]
+ if len(dcm_axes) != 0:
+ raise ValueError('No character in voxel_order corresponding to '
+ 'axes: %s' % dcm_axes)
+
+ #Check the vox_array and affine have correct shape/size
+ if len(vox_array.shape) < 3:
+ raise ValueError('The vox_array must be at least three dimensional')
+ if affine.shape != (4, 4):
+ raise ValueError('The affine must be 4x4')
+
+ #Pull the current index directions from the affine
+ orig_ornt = io_orientation(affine)
+ new_ornt = axcodes2ornt(voxel_order)
+ ornt_trans = ornt_transform(orig_ornt, new_ornt)
+ orig_shape = vox_array.shape
+ vox_array = apply_orientation(vox_array, ornt_trans)
+ aff_trans = inv_ornt_aff(ornt_trans, orig_shape)
+ affine = np.dot(affine, aff_trans)
+
+ return (vox_array, affine, aff_trans, ornt_trans)
+
+def dcm_time_to_sec(time_str):
+ '''Convert a DICOM time value (value representation of 'TM') to the number
+ of seconds past midnight.
+
+ Parameters
+ ----------
+ time_str : str
+ The DICOM time value string
+
+ Returns
+ -------
+ A floating point representing the number of seconds past midnight
+ '''
+ #Allow ACR/NEMA style format by removing any colon chars
+ time_str = time_str.replace(':', '')
+
+ #Only the hours portion is required
+ result = int(time_str[:2]) * 3600
+
+ str_len = len(time_str)
+ if str_len > 2:
+ result += int(time_str[2:4]) * 60
+ if str_len > 4:
+ result += float(time_str[4:])
+
+ return float(result)
+
+class IncongruentImageError(Exception):
+ def __init__(self, msg):
+ '''An exception denoting that a DICOM with incorrect size or orientation
+ was passed to `DicomStack.add_dcm`.'''
+ self.msg = msg
+
+ def __str__(self):
+ return 'The image is not congruent to the existing stack: %s' % self.msg
+
+class ImageCollisionError(Exception):
+ '''An exception denoting that a DICOM which collides with one already in
+ the stack was passed to a `DicomStack.add_dcm`.'''
+ def __str__(self):
+ return 'The image collides with one already in the stack'
+
+class InvalidStackError(Exception):
+ def __init__(self, msg):
+ '''An exception denoting that a `DicomStack` is not currently valid'''
+ self.msg = msg
+
+ def __str__(self):
+ return 'The DICOM stack is not valid: %s' % self.msg
+
+class DicomOrdering(object):
+ '''Object defining an ordering for a set of dicom datasets. Create a
+ DicomOrdering with the given DICOM element keyword.
+
+ Parameters
+ ----------
+ key : str
+ The DICOM keyword to use for ordering the datasets
+
+ abs_ordering : sequence
+ A sequence specifying the absolute order for values corresponding
+ to the `key`. Instead of ordering by the value associated with the
+ `key`, the index of the value in this sequence will be used.
+
+ abs_as_str : bool
+ If true, the values will be converted to strings before looking up
+ the index in `abs_ordering`.
+ '''
+
+ def __init__(self, key, abs_ordering=None, abs_as_str=False):
+ self.key = key
+ self.abs_ordering = abs_ordering
+ self.abs_as_str = abs_as_str
+
+ def get_ordinate(self, ds):
+ '''Get the ordinate for the given DICOM data set.
+
+ Parameters
+ ----------
+ ds : dict like
+ The DICOM data set we want the ordinate of. Should allow
+ dict like access where DICOM keywords return the corresponing
+ value.
+
+ Returns
+ -------
+ An ordinate for the data set. If `abs_ordering` is None then this will
+ just be the value for the keyword `key`. Otherwise it will be an
+ integer.
+ '''
+ try:
+ val = ds[self.key]
+ except KeyError:
+ return None
+
+ if self.abs_ordering:
+ if self.abs_as_str:
+ val = str(val)
+ return self.abs_ordering.index(val)
+
+ return val
+
+def _make_dummy(reference, meta, iop):
+ '''Make a "dummy" NiftiWrapper (no valid pixel data).'''
+ #Create the dummy data array filled with largest representable value
+ data = np.empty_like(reference.nii_img.get_data())
+ data[...] = np.iinfo(np.int16).max
+
+ #Create the nifti image and set header data
+ aff = reference.nii_img.get_affine().copy()
+ aff[:3, 3] = [iop[1], iop[0], iop[2]]
+ nii_img = nb.nifti1.Nifti1Image(data, aff)
+ hdr = nii_img.get_header()
+ hdr.set_xyzt_units('mm', 'sec')
+ dim_info = {'freq' : None,
+ 'phase' : None,
+ 'slice' : 2
+ }
+ if 'InplanePhaseEncodingDirection' in meta:
+ if meta['InplanePhaseEncodingDirection'] == 'ROW':
+ dim_info['phase'] = 1
+ dim_info['freq'] = 0
+ else:
+ dim_info['phase'] = 0
+ dim_info['freq'] = 1
+ hdr.set_dim_info(**dim_info)
+
+ #Embed the meta data extension
+ result = NiftiWrapper(nii_img, make_empty=True)
+ result.meta_ext.reorient_transform = np.diag([-1., -1., 1., 1.])
+ result.meta_ext.get_class_dict(('global', 'const')).update(meta)
+
+ return result
+
+default_group_keys = ('SeriesInstanceUID',
+ 'SeriesNumber',
+ 'ProtocolName',
+ 'ImageOrientationPatient')
+'''Default keys for grouping DICOM files that belong in the same
+multi-dimensional array together.'''
+
+class DicomStack(object):
+ '''Defines a method for stacking together DICOM data sets into a multi
+ dimensional volume.
+
+ Tailored towards creating NiftiImage output, but can also just create numpy
+ arrays. Can summarize all of the meta data from the input DICOM data sets
+ into a Nifti header extension (see `dcmmeta.DcmMetaExtension`).
+
+ Parameters
+ ----------
+ time_order : str or DicomOrdering
+ The DICOM keyword or DicomOrdering object specifying how to order
+ the DICOM data sets along the time dimension.
+
+ vector_order : str or DicomOrdering
+ The DICOM keyword or DicomOrdering object specifying how to order
+ the DICOM data sets along the vector dimension.
+
+ allow_dummies : bool
+ If True then data sets without pixel data can be added to the stack.
+ The "dummy" voxels will have the maximum representable value for
+ the datatype.
+
+ meta_filter : callable
+ A callable that takes a meta data key and value, and returns True if
+ that meta data element should be excluded from the DcmMeta extension.
+
+ Notes
+ -----
+ If both time_order and vector_order are None, the time_order will be
+ guessed based off the data sets.
+ '''
+
+ sort_guesses = ['EchoTime',
+ 'InversionTime',
+ 'RepetitionTime',
+ 'FlipAngle',
+ 'TriggerTime',
+ 'AcquisitionTime',
+ 'ContentTime',
+ 'AcquisitionNumber',
+ 'InstanceNumber',
+ ]
+ '''The meta data keywords used when trying to guess the sorting order.
+ Keys that come earlier in the list are given higher priority.'''
+
+ minimal_keys = set(sort_guesses +
+ ['Rows',
+ 'Columns',
+ 'PixelSpacing',
+ 'ImageOrientationPatient',
+ 'InPlanePhaseEncodingDirection',
+ 'RepetitionTime',
+ 'AcquisitionTime'
+ ] +
+ list(default_group_keys)
+ )
+ '''Set of minimal meta data keys that should be provided if they exist in
+ the source DICOM files.'''
+
+ def __init__(self, time_order=None, vector_order=None,
+ allow_dummies=False, meta_filter=None):
+ if isinstance(time_order, str):
+ self._time_order = DicomOrdering(time_order)
+ else:
+ self._time_order = time_order
+ if isinstance(vector_order, str):
+ self._vector_order = DicomOrdering(vector_order)
+ else:
+ self._vector_order = vector_order
+
+ if meta_filter is None:
+ self._meta_filter = default_meta_filter
+ else:
+ self._meta_filter = meta_filter
+
+ self._allow_dummies = allow_dummies
+
+ #Sets all the state variables to their defaults
+ self.clear()
+
+ def _chk_equal(self, keys, meta1, meta2):
+ for key in keys:
+ if meta1[key] != meta2[key]:
+ raise IncongruentImageError("%s does not match" % key)
+
+ def _chk_close(self, keys, meta1, meta2):
+ for key in keys:
+ if not np.allclose(meta1[key], meta2[key], atol=5e-5):
+ raise IncongruentImageError("%s is not close to matching" %
+ key)
+
+ def _chk_congruent(self, meta):
+ is_dummy = not 'Rows' in meta or not 'Columns' in meta
+ if is_dummy and not self._allow_dummies:
+ raise IncongruentImageError('Missing Rows/Columns')
+
+ if not self._ref_input is None:
+ self._chk_close(('PixelSpacing',
+ 'ImageOrientationPatient'),
+ meta,
+ self._ref_input
+ )
+ if not is_dummy:
+ self._chk_equal(('Rows', 'Columns'), meta, self._ref_input)
+ elif len(self._dummies) != 0:
+ self._chk_close(('PixelSpacing',
+ 'ImageOrientationPatient'),
+ meta,
+ self._dummies[0][0]
+ )
+ return is_dummy
+
+
+ def add_dcm(self, dcm, meta=None):
+ '''Add a pydicom dataset to the stack.
+
+ Parameters
+ ----------
+ dcm : dicom.dataset.Dataset
+ The data set being added to the stack
+
+ meta : dict
+ The extracted meta data for the DICOM data set `dcm`. If None
+ extract.default_extractor will be used.
+
+ Raises
+ ------
+ IncongruentImageError
+ The provided `dcm` does not match the orientation or dimensions of
+ those already in the stack.
+
+ ImageCollisionError
+ The provided `dcm` has the same slice location and time/vector
+ values.
+
+ '''
+ if meta is None:
+ from .extract import default_extractor
+ meta = default_extractor(dcm)
+
+ dw = wrapper_from_data(dcm)
+
+ is_dummy = self._chk_congruent(meta)
+
+ self._phase_enc_dirs.add(meta.get('InPlanePhaseEncodingDirection'))
+ self._repetition_times.add(meta.get('RepetitionTime'))
+
+ #Pull the info used for sorting
+ slice_pos = dw.slice_indicator
+ self._slice_pos_vals.add(slice_pos)
+ time_val = None
+ if self._time_order:
+ time_val = self._time_order.get_ordinate(meta)
+ self._time_vals.add(time_val)
+ vector_val = None
+ if self._vector_order:
+ vector_val = self._vector_order.get_ordinate(meta)
+ self._vector_vals.add(vector_val)
+
+ #Create a tuple with the sorting values
+ sorting_tuple = (vector_val, time_val, slice_pos)
+
+ #If a explicit order was specified, raise an exception if image
+ #collides with another already in the stack
+ if ((not self._time_order is None or
+ not self._vector_order is None) and
+ sorting_tuple in self._sorting_tuples
+ ):
+ raise ImageCollisionError()
+ self._sorting_tuples.add(sorting_tuple)
+
+ #Create a NiftiWrapper for this input if possible
+ nii_wrp = None
+ if not is_dummy:
+ nii_wrp = NiftiWrapper.from_dicom_wrapper(dw, meta)
+ if self._ref_input is None:
+ #We don't have a reference input yet, use this one
+ self._ref_input = nii_wrp
+ #Convert any dummies that we have stashed previously
+ for dummy_meta, dummy_tuple, iop in self._dummies:
+ dummy_wrp = _make_dummy(self._ref_input, dummy_meta, iop)
+ self._files_info.append((dummy_wrp, dummy_tuple))
+ else:
+ if self._ref_input is None:
+ #We don't have a reference input, so stash the dummy for now
+ self._dummies.append((meta, sorting_tuple, dcm.ImagePositionPatient))
+ else:
+ #Convert dummy using the reference input
+ nii_wrp = _make_dummy(self._ref_input, meta, dcm.ImagePositionPatient)
+
+ #If we made a NiftiWrapper add it to the stack
+ if not nii_wrp is None:
+ self._files_info.append((nii_wrp, sorting_tuple))
+
+ #Set the dirty flags
+ self._shape_dirty = True
+ self._meta_dirty = True
+
+ def clear(self):
+ '''Remove any DICOM datasets from the stack.'''
+ self._slice_pos_vals = set()
+ self._time_vals = set()
+ self._vector_vals = set()
+ self._sorting_tuples = set()
+
+ self._phase_enc_dirs = set()
+ self._repetition_times = set()
+
+ self._dummies = []
+ self._ref_input = None
+
+ self._shape_dirty = True
+ self._shape = None
+
+ self._meta_dirty = True
+ self._meta = None
+
+ self._files_info = []
+
+ def _chk_order(self, slice_positions, files_per_vol, num_volumes,
+ num_time_points, num_vec_comps):
+ #Sort the files
+ self._files_info.sort(key=lambda x: x[1])
+ if files_per_vol > 1:
+ for vol_idx in range(num_volumes):
+ start_slice = vol_idx * files_per_vol
+ end_slice = start_slice + files_per_vol
+ self._files_info[start_slice:end_slice] = \
+ sorted(self._files_info[start_slice:end_slice],
+ key=lambda x: x[1][-1])
+
+ #Do a thorough check for correctness
+ for vec_idx in xrange(num_vec_comps):
+ file_idx = vec_idx*num_time_points*files_per_vol
+ curr_vec_val = self._files_info[file_idx][1][0]
+ for time_idx in xrange(num_time_points):
+ for slice_idx in xrange(files_per_vol):
+ file_idx = (vec_idx*num_time_points*files_per_vol +
+ time_idx*files_per_vol + slice_idx)
+ file_info = self._files_info[file_idx]
+ if file_info[1][0] != curr_vec_val:
+ raise InvalidStackError("Not enough images with the " +
+ "vector value of " +
+ str(curr_vec_val))
+ if (file_info[1][2] != slice_positions[slice_idx]):
+ if (file_info[1][2] == slice_positions[slice_idx-1]):
+ error_msg = ["Duplicate slice position"]
+ else:
+ error_msg = ["Missing slice position"]
+ error_msg.append(" at slice index %d" % slice_idx)
+ if num_time_points > 1:
+ error_msg.append(' in time point %d' % time_idx)
+ if num_vec_comps > 1:
+ error_msg.append(' for vector component %s' %
+ str(curr_vec_val))
+ raise InvalidStackError(''.join(error_msg))
+
+
+ def get_shape(self):
+ '''Get the shape of the stack.
+
+ Returns
+ -------
+ A tuple of integers giving the size of the dimensions of the stack.
+
+ Raises
+ ------
+ InvalidStackError
+ The stack is incomplete or invalid.
+ '''
+ #If the dirty flag is not set, return the cached value
+ if not self._shape_dirty:
+ return self._shape
+
+ #We need at least one non-dummy file in the stack
+ if len(self._files_info) == 0:
+ raise InvalidStackError("No (non-dummy) files in the stack")
+
+ #Figure out number of files and slices per volume
+ files_per_vol = len(self._slice_pos_vals)
+ slice_positions = sorted(list(self._slice_pos_vals))
+
+ #If more than one file per volume, check that slice spacing is equal
+ if files_per_vol > 1:
+ spacings = []
+ for idx in xrange(files_per_vol - 1):
+ spacings.append(slice_positions[idx+1] - slice_positions[idx])
+ spacings = np.array(spacings)
+ avg_spacing = np.mean(spacings)
+ if not np.allclose(avg_spacing, spacings, rtol=4e-2):
+ raise InvalidStackError("Slice spacings are not consistent")
+
+ #Simple check for an incomplete stack
+ if len(self._files_info) % files_per_vol != 0:
+ raise InvalidStackError("Number of files is not an even multiple "
+ "of the number of unique slice positions.")
+ num_volumes = len(self._files_info) / files_per_vol
+
+ #Figure out the number of vector components and time points
+ num_vec_comps = len(self._vector_vals)
+ if num_vec_comps > num_volumes:
+ raise InvalidStackError("Vector variable varies within volumes")
+ if num_volumes % num_vec_comps != 0:
+ raise InvalidStackError("Number of volumes not an even multiple "
+ "of the number of vector components.")
+ num_time_points = num_volumes / num_vec_comps
+
+ #If both sort keys are None try to guess
+ if (num_volumes > 1 and self._time_order is None and
+ self._vector_order is None):
+ #Get a list of possible sort orders
+ possible_orders = []
+ for key in self.sort_guesses:
+ vals = set([file_info[0].get_meta(key)
+ for file_info in self._files_info]
+ )
+ if len(vals) == num_volumes or len(vals) == len(self._files_info):
+ possible_orders.append(key)
+ if len(possible_orders) == 0:
+ raise InvalidStackError("Unable to guess key for sorting the "
+ "fourth dimension")
+
+ #Try out each possible sort order
+ for time_order in possible_orders:
+ #Update sorting tuples
+ for idx in xrange(len(self._files_info)):
+ nii_wrp, curr_tuple = self._files_info[idx]
+ self._files_info[idx] = (nii_wrp,
+ (curr_tuple[0],
+ nii_wrp[time_order],
+ curr_tuple[2]
+ )
+ )
+
+ #Check the order
+ try:
+ self._chk_order(slice_positions,
+ files_per_vol,
+ num_volumes,
+ num_time_points,
+ num_vec_comps)
+ except InvalidStackError:
+ pass
+ else:
+ break
+ else:
+ raise InvalidStackError("Unable to guess key for sorting the "
+ "fourth dimension")
+ else: #If at least on sort key was specified, just check the order
+ self._chk_order(slice_positions,
+ files_per_vol,
+ num_volumes,
+ num_time_points,
+ num_vec_comps)
+
+ #Stack appears to be valid, build the shape tuple
+ file_shape = self._files_info[0][0].nii_img.get_shape()
+ vol_shape = list(file_shape)
+ if files_per_vol > 1:
+ vol_shape[2] = files_per_vol
+ shape = vol_shape+ [num_time_points, num_vec_comps]
+ if shape[4] == 1:
+ shape = shape[:-1]
+ if shape[3] == 1:
+ shape = shape[:-1]
+ self._shape = tuple(shape)
+
+ self._shape_dirty = False
+ return self._shape
+
+ def get_data(self):
+ '''Get an array of the voxel values.
+
+ Returns
+ -------
+ A numpy array filled with values from the DICOM data sets' pixels.
+
+ Raises
+ ------
+ InvalidStackError
+ The stack is incomplete or invalid.
+ '''
+ #Create a numpy array for storing the voxel data
+ stack_shape = self.get_shape()
+ stack_shape = tuple(list(stack_shape) + ((5 - len(stack_shape)) * [1]))
+ stack_dtype = self._files_info[0][0].nii_img.get_data_dtype()
+ #This is a hack to keep fslview happy, Shouldn't cause issues as the
+ #original data should be 12-bit and any scaling will result in float
+ #data
+ if stack_dtype == np.uint16:
+ stack_dtype = np.int16
+ vox_array = np.empty(stack_shape, dtype=stack_dtype)
+
+ #Fill the array with data
+ n_vols = 1
+ if len(stack_shape) > 3:
+ n_vols *= stack_shape[3]
+ if len(stack_shape) > 4:
+ n_vols *= stack_shape[4]
+ files_per_vol = len(self._files_info) / n_vols
+ file_shape = self._files_info[0][0].nii_img.get_shape()
+ for vec_idx in range(stack_shape[4]):
+ for time_idx in range(stack_shape[3]):
+ if files_per_vol == 1 and file_shape[2] != 1:
+ file_idx = vec_idx*(stack_shape[3]) + time_idx
+ vox_array[:, :, :, time_idx, vec_idx] = \
+ self._files_info[file_idx][0].nii_img.get_data()
+ else:
+ for slice_idx in range(files_per_vol):
+ file_idx = (vec_idx*(stack_shape[3]*stack_shape[2]) +
+ time_idx*(stack_shape[2]) + slice_idx)
+ vox_array[:, :, slice_idx, time_idx, vec_idx] = \
+ self._files_info[file_idx][0].nii_img.get_data()[:, :, 0]
+
+ #Trim unused time/vector dimensions
+ if stack_shape[4] == 1:
+ vox_array = vox_array[...,0]
+ if stack_shape[3] == 1:
+ vox_array = vox_array[...,0]
+
+ return vox_array
+
+ def get_affine(self):
+ '''Get the affine transform for mapping row/column/slice indices
+ to Nifti (RAS) patient space.
+
+ Returns
+ -------
+ A 4x4 numpy array containing the affine transform.
+
+ Raises
+ ------
+ InvalidStackError
+ The stack is incomplete or invalid.
+ '''
+ #Figure out the number of three (or two) dimensional volumes
+ shape = self.get_shape()
+ n_vols = 1
+ if len(shape) > 3:
+ n_vols *= shape[3]
+ if len(shape) > 4:
+ n_vols *= shape[4]
+
+ #Figure out the number of files in each volume
+ files_per_vol = len(self._files_info) / n_vols
+
+ #Pull the DICOM Patient Space affine from the first input
+ aff = self._files_info[0][0].nii_img.get_affine()
+
+ #If there is more than one file per volume, we need to fix slice scaling
+ if files_per_vol > 1:
+ first_offset = aff[:3, 3]
+ second_offset = self._files_info[1][0].nii_img.get_affine()[:3, 3]
+ scaled_slc_dir = second_offset - first_offset
+ aff[:3, 2] = scaled_slc_dir
+
+ return aff
+
+ def to_nifti(self, voxel_order='LAS', embed_meta=False):
+ '''Returns a NiftiImage with the data and affine from the stack.
+
+ Parameters
+ ----------
+ voxel_order : str
+ A three character string repsenting the voxel order in patient
+ space (see the function `reorder_voxels`). Can be None or an empty
+ string to disable reorientation.
+
+ embed_meta : bool
+ If true a dcmmeta.DcmMetaExtension will be embedded in the Nifti
+ header.
+
+ Returns
+ -------
+ A nibabel.nifti1.Nifti1Image created with the stack's data and affine.
+ '''
+ #Get the voxel data and affine
+ data = self.get_data()
+ affine = self.get_affine()
+
+ #Figure out the number of three (or two) dimensional volumes
+ n_vols = 1
+ if len(data.shape) > 3:
+ n_vols *= data.shape[3]
+ if len(data.shape) > 4:
+ n_vols *= data.shape[4]
+
+ files_per_vol = len(self._files_info) / n_vols
+
+ #Reorder the voxel data if requested
+ permutation = [0, 1, 2]
+ slice_dim = 2
+ reorient_transform = np.eye(4)
+ if voxel_order:
+ (data,
+ affine,
+ reorient_transform,
+ ornt_trans) = reorder_voxels(data, affine, voxel_order)
+ permutation, flips = zip(*ornt_trans)
+ permutation = [int(val) for val in permutation]
+
+ #Reverse file order in each volume's files if we flipped slice order
+ #This will keep the slice times and meta data order correct
+ if files_per_vol > 1 and flips[slice_dim] == -1:
+ self._shape_dirty = True
+ for vol_idx in xrange(n_vols):
+ start = vol_idx * files_per_vol
+ stop = start + files_per_vol
+ self._files_info[start:stop] = [self._files_info[idx]
+ for idx in xrange(stop - 1,
+ start - 1,
+ -1)
+ ]
+
+ #Update the slice dim
+ slice_dim = permutation[2]
+
+ #Create the nifti image using the data array
+ nifti_image = nb.Nifti1Image(data, affine)
+ nifti_header = nifti_image.get_header()
+
+ #Set the units and dimension info
+ nifti_header.set_xyzt_units('mm', 'msec')
+ if len(self._repetition_times) == 1 and not None in self._repetition_times:
+ nifti_header['pixdim'][4] = list(self._repetition_times)[0]
+ dim_info = {'freq' : None, 'phase' : None, 'slice' : slice_dim}
+ if len(self._phase_enc_dirs) == 1 and not None in self._phase_enc_dirs:
+ phase_dir = list(self._phase_enc_dirs)[0]
+ if phase_dir == 'ROW':
+ dim_info['phase'] = permutation[1]
+ dim_info['freq'] = permutation[0]
+ else:
+ dim_info['phase'] = permutation[0]
+ dim_info['freq'] = permutation[1]
+ nifti_header.set_dim_info(**dim_info)
+ n_slices = data.shape[slice_dim]
+
+ #Set the slice timing header info
+ has_acq_time = (self._files_info[0][0].get_meta('AcquisitionTime') !=
+ None)
+ if files_per_vol > 1 and has_acq_time:
+ #Pull out the relative slice times for the first volume
+ slice_times = np.array([dcm_time_to_sec(file_info[0]['AcquisitionTime'])
+ for file_info in self._files_info[:n_slices]]
+ )
+ slice_times -= np.min(slice_times)
+
+ #If there is more than one volume, check if times are consistent
+ is_consistent = True
+ for vol_idx in xrange(1, n_vols):
+ start_slice = vol_idx * n_slices
+ end_slice = start_slice + n_slices
+ slices_info = self._files_info[start_slice:end_slice]
+ vol_slc_times = \
+ np.array([dcm_time_to_sec(file_info[0]['AcquisitionTime'])
+ for file_info in slices_info]
+ )
+ vol_slc_times -= np.min(vol_slc_times)
+ if not np.allclose(slice_times, vol_slc_times):
+ is_consistent = False
+ break
+
+ #If the times are consistent and not all zero, try setting the slice
+ #times (sets the slice duration and code if possible).
+ if is_consistent and not np.allclose(slice_times, 0.0):
+ try:
+ nifti_header.set_slice_times(slice_times)
+ except HeaderDataError:
+ pass
+
+ #Embed the meta data extension if requested
+ if embed_meta:
+ #Build meta data for each volume if needed
+ vol_meta = []
+ if files_per_vol > 1:
+ for vol_idx in xrange(n_vols):
+ start_slice = vol_idx * n_slices
+ end_slice = start_slice + n_slices
+ exts = [file_info[0].meta_ext
+ for file_info in self._files_info[start_slice:end_slice]]
+ meta = DcmMetaExtension.from_sequence(exts, 2)
+ vol_meta.append(meta)
+ else:
+ vol_meta = [file_info[0].meta_ext
+ for file_info in self._files_info]
+
+ #Build meta data for each time point / vector component
+ if len(data.shape) == 5:
+ if data.shape[3] != 1:
+ vec_meta = []
+ for vec_idx in xrange(data.shape[4]):
+ start_idx = vec_idx * data.shape[3]
+ end_idx = start_idx + data.shape[3]
+ meta = DcmMetaExtension.from_sequence(\
+ vol_meta[start_idx:end_idx], 3)
+ vec_meta.append(meta)
+ else:
+ vec_meta = vol_meta
+
+ meta_ext = DcmMetaExtension.from_sequence(vec_meta, 4)
+ elif len(data.shape) == 4:
+ meta_ext = DcmMetaExtension.from_sequence(vol_meta, 3)
+ else:
+ meta_ext = vol_meta[0]
+ if meta_ext is file_info[0].meta_ext:
+ meta_ext = deepcopy(meta_ext)
+
+ meta_ext.shape = data.shape
+ meta_ext.slice_dim = slice_dim
+ meta_ext.affine = nifti_header.get_best_affine()
+ meta_ext.reorient_transform = reorient_transform
+
+ #Filter and embed the meta data
+ meta_ext.filter_meta(self._meta_filter)
+ nifti_header.extensions = Nifti1Extensions([meta_ext])
+
+ nifti_image.update_header()
+ return nifti_image
+
+ def to_nifti_wrapper(self, voxel_order=''):
+ '''Convienance method. Calls `to_nifti` and returns a `NiftiWrapper`
+ generated from the result.
+ '''
+ return NiftiWrapper(self.to_nifti(voxel_order, True))
+
+def parse_and_group(src_paths, group_by=default_group_keys, extractor=None,
+ force=False, warn_on_except=False,
+ close_tests=('ImageOrientationPatient',)):
+ '''Parse the given dicom files and group them together. Each group is
+ stored as a (list) value in a dict where the key is a tuple of values
+ corresponding to the keys in 'group_by'
+
+ Parameters
+ ----------
+ src_paths : sequence
+ A list of paths to the source DICOM files.
+
+ group_by : tuple
+ Meta data keys to group data sets with. Any data set with the same
+ values for these keys will be grouped together. This tuple of values
+ will also be the key in the result dictionary.
+
+ extractor : callable
+ Should take a dicom.dataset.Dataset and return a dictionary of the
+ extracted meta data.
+
+ force : bool
+ Force reading source files even if they do not appear to be DICOM.
+
+ warn_on_except : bool
+ Convert exceptions into warnings, possibly allowing some results to be
+ returned.
+
+ close_tests : sequence
+ Any `group_by` key listed here is tested with `numpy.allclose` instead
+ of straight equality when determining group membership.
+
+ Returns
+ -------
+ groups : dict
+ A dict mapping tuples of values (corresponding to 'group_by') to groups
+ of data sets. Each element in the list is a tuple containing the dicom
+ object, the parsed meta data, and the filename.
+ '''
+ if extractor is None:
+ from .extract import default_extractor
+ extractor = default_extractor
+
+ results = {}
+ close_elems = {}
+ for dcm_path in src_paths:
+ #Read the DICOM file
+ try:
+ dcm = dicom.read_file(dcm_path, force=force)
+ except Exception as e:
+ if warn_on_except:
+ warnings.warn('Error reading file %s: %s' % (dcm_path, str(e)))
+ continue
+ else:
+ raise
+
+ #Extract the meta data and group
+ meta = extractor(dcm)
+ key_list = [] # Values from group_by elems with equality testing
+ close_list = [] # Values from group_by elems with np.allclose testing
+ for grp_key in group_by:
+ key_elem = meta.get(grp_key)
+ if isinstance(key_elem, list):
+ key_elem = tuple(key_elem)
+ if grp_key in close_tests:
+ close_list.append(key_elem)
+ else:
+ key_list.append(key_elem)
+
+ # Initially each key has multiple sub_results (corresponding to
+ # different values of the "close" keys)
+ key = tuple(key_list)
+ if not key in results:
+ results[key] = [(close_list, [(dcm, meta, dcm_path)])]
+ else:
+ # Look for a matching sub_result
+ for c_list, sub_res in results[key]:
+ for c_idx, c_val in enumerate(c_list):
+ if not np.allclose(c_val, close_list[c_idx], atol=5e-5):
+ break
+ else:
+ sub_res.append((dcm, meta, dcm_path))
+ break
+ else:
+ # No match found, append another sub result
+ results[key].append((close_list, [(dcm, meta, dcm_path)]))
+
+ # Unpack sub results, using the canonical value for the close keys
+ full_results = {}
+ for eq_key, sub_res_list in results.iteritems():
+ for close_key, sub_res in sub_res_list:
+ full_key = []
+ eq_idx = 0
+ close_idx = 0
+ for grp_key in group_by:
+ if grp_key in close_tests:
+ full_key.append(close_key[close_idx])
+ close_idx += 1
+ else:
+ full_key.append(eq_key[eq_idx])
+ eq_idx += 1
+ full_key = tuple(full_key)
+ full_results[full_key] = sub_res
+
+ return full_results
+
+def stack_group(group, warn_on_except=False, **stack_args):
+ result = DicomStack(**stack_args)
+ for dcm, meta, fn in group:
+ try:
+ result.add_dcm(dcm, meta)
+ except Exception as e:
+ if warn_on_except:
+ warnings.warn('Error adding file %s to stack: %s' %
+ (fn, str(e)))
+ else:
+ raise
+ return result
+
+def parse_and_stack(src_paths, group_by=default_group_keys, extractor=None,
+ force=False, warn_on_except=False, **stack_args):
+ '''Parse the given dicom files into a dictionary containing one or more
+ DicomStack objects.
+
+ Parameters
+ ----------
+ src_paths : sequence
+ A list of paths to the source DICOM files.
+
+ group_by : tuple
+ Meta data keys to group data sets with. Any data set with the same
+ values for these keys will be grouped together. This tuple of values
+ will also be the key in the result dictionary.
+
+ extractor : callable
+ Should take a dicom.dataset.Dataset and return a dictionary of the
+ extracted meta data.
+
+ force : bool
+ Force reading source files even if they do not appear to be DICOM.
+
+ warn_on_except : bool
+ Convert exceptions into warnings, possibly allowing some results to be
+ returned.
+
+ stack_args : kwargs
+ Keyword arguments to pass to the DicomStack constructor.
+ '''
+ results = parse_and_group(src_paths,
+ group_by,
+ extractor,
+ force,
+ warn_on_except)
+
+ for key, group in results.iteritems():
+ results[key] = stack_group(group, warn_on_except, **stack_args)
+
+ return results
diff --git a/fsl/data/dcmstack/dcmstack_cli.py b/fsl/data/dcmstack/dcmstack_cli.py
new file mode 100755
index 000000000..1c041a9cb
--- /dev/null
+++ b/fsl/data/dcmstack/dcmstack_cli.py
@@ -0,0 +1,357 @@
+"""
+Command line interface to dcmstack.
+
+@author: moloney
+"""
+from __future__ import print_function
+
+import os, sys, argparse, string
+from glob import glob
+import dicom
+from . import dcmstack
+from .dcmstack import (parse_and_group, stack_group, DicomOrdering,
+ default_group_keys)
+from .dcmmeta import NiftiWrapper
+from . import extract
+from .info import __version__
+
+
+prog_descrip = """Stack DICOM files from each source directory into 2D to 5D
+volumes, optionally extracting meta data.
+"""
+
+prog_epilog = """IT IS YOUR RESPONSIBILITY TO KNOW IF THERE IS PRIVATE HEALTH
+INFORMATION IN THE METADATA EXTRACTED BY THIS PROGRAM."""
+
+def parse_tags(opt_str):
+ tag_strs = opt_str.split(',')
+ tags = []
+ for tag_str in tag_strs:
+ tokens = tag_str.split('_')
+ if len(tokens) != 2:
+ raise ValueError('Invalid str format for tags')
+ tags.append(dicom.tag.Tag(int(tokens[0].strip(), 16),
+ int(tokens[1].strip(), 16))
+ )
+ return tags
+
+def sanitize_path_comp(path_comp):
+ result = []
+ for char in path_comp:
+ if not char in string.letters + string.digits + '-_.':
+ result.append('_')
+ else:
+ result.append(char)
+ return ''.join(result)
+
+def main(argv=sys.argv):
+ #Handle command line options
+ arg_parser = argparse.ArgumentParser(description=prog_descrip,
+ epilog=prog_epilog)
+ arg_parser.add_argument('src_dirs', nargs='*', help=('The source '
+ 'directories containing DICOM files.'))
+
+ input_opt = arg_parser.add_argument_group('Input options')
+ input_opt.add_argument('--force-read', action='store_true', default=False,
+ help=('Try reading all files as DICOM, even if they '
+ 'are missing the preamble.'))
+ input_opt.add_argument('--file-ext', default='.dcm', help=('Only try reading '
+ 'files with the given extension. Default: '
+ '%(default)s'))
+ input_opt.add_argument('--allow-dummies', action='store_true', default=False,
+ help=('Allow DICOM files that are missing pixel '
+ 'data, filling that slice of the output nifti with '
+ 'the maximum representable value.'))
+
+ output_opt = arg_parser.add_argument_group('Output options')
+ output_opt.add_argument('--dest-dir', default=None,
+ help=('Destination directory, defaults to the '
+ 'source directory.'))
+ output_opt.add_argument('-o', '--output-name', default=None,
+ help=('Python format string determining the output '
+ 'filenames based on DICOM tags.'))
+ output_opt.add_argument('--output-ext', default='.nii.gz',
+ help=('The extension for the output file type. '
+ 'Default: %(default)s'))
+ output_opt.add_argument('-d', '--dump-meta', default=False,
+ action='store_true', help=('Dump the extracted '
+ 'meta data into a JSON file with the same base '
+ 'name as the generated Nifti'))
+ output_opt.add_argument('--embed-meta', default=False, action='store_true',
+ help=('Embed the extracted meta data into a Nifti '
+ 'header extension (in JSON format).'))
+
+ stack_opt = arg_parser.add_argument_group('Stacking Options')
+ stack_opt.add_argument('-g', '--group-by', default=None,
+ help=("Comma separated list of meta data keys to "
+ "group input files into stacks with."))
+ stack_opt.add_argument('--voxel-order', default='LAS',
+ help=('Order the voxels so the spatial indices '
+ 'start from these directions in patient space. '
+ 'The directions in patient space should be given '
+ 'as a three character code: (l)eft, (r)ight, '
+ '(a)nterior, (p)osterior, (s)uperior, (i)nferior. '
+ 'Passing an empty string will disable '
+ 'reorientation. '
+ 'Default: %(default)s'))
+ stack_opt.add_argument('-t', '--time-var', default=None,
+ help=('The DICOM element keyword to use for '
+ 'ordering the stack along the time dimension.'))
+ stack_opt.add_argument('--vector-var', default=None,
+ help=('The DICOM element keyword to use for '
+ 'ordering the stack along the vector dimension.'))
+ stack_opt.add_argument('--time-order', default=None,
+ help=('Provide a text file with the desired order '
+ 'for the values (one per line) of the attribute '
+ 'used as the time variable. This option is rarely '
+ 'needed.'))
+ stack_opt.add_argument('--vector-order', default=None,
+ help=('Provide a text file with the desired order '
+ 'for the values (one per line) of the attribute '
+ 'used as the vector variable. This option is rarely '
+ 'needed.'))
+
+ meta_opt = arg_parser.add_argument_group('Meta Extraction and Filtering '
+ 'Options')
+ meta_opt.add_argument('-l', '--list-translators', default=False,
+ action='store_true', help=('List enabled translators '
+ 'and exit'))
+ meta_opt.add_argument('--disable-translator', default=None,
+ help=('Disable the translators for the provided '
+ 'tags. Tags should be given in the format '
+ '"0x0_0x0". More than one can be given in a comma '
+ 'separated list. If the word "all" is provided, all '
+ 'translators will be disabled.'))
+ meta_opt.add_argument('--extract-private', default=False,
+ action='store_true',
+ help=('Extract meta data from private elements, even '
+ 'if there is no translator. If the value for the '
+ 'element contains non-ascii bytes it will still be '
+ 'ignored. The extracted meta data may still be '
+ 'filtered out by the regular expressions.'))
+ meta_opt.add_argument('-i', '--include-regex', action='append',
+ help=('Include any meta data where the key matches '
+ 'the provided regular expression. This will override '
+ 'any exclude expressions. Applies to all meta data.'))
+ meta_opt.add_argument('-e', '--exclude-regex', action='append',
+ help=('Exclude any meta data where the key matches '
+ 'the provided regular expression. This will '
+ 'supplement the default exclude expressions. Applies '
+ 'to all meta data.'))
+ meta_opt.add_argument('--default-regexes', default=False,
+ action='store_true',
+ help=('Print the list of default include and exclude '
+ 'regular expressions and exit.'))
+
+ gen_opt = arg_parser.add_argument_group('General Options')
+ gen_opt.add_argument('-v', '--verbose', default=False, action='store_true',
+ help=('Print additional information.'))
+ gen_opt.add_argument('--strict', default=False, action='store_true',
+ help=('Fail on the first exception instead of '
+ 'showing a warning.'))
+ gen_opt.add_argument('--version', default=False, action='store_true',
+ help=('Show the version and exit.'))
+
+ args = arg_parser.parse_args(argv[1:])
+
+ if args.version:
+ print(__version__)
+ return 0
+
+ #Check if we are just listing the translators
+ if args.list_translators:
+ for translator in extract.default_translators:
+ print('%s -> %s' % (translator.tag, translator.name))
+ return 0
+
+ #Check if we are just listing the default exclude regular expressions
+ if args.default_regexes:
+ print('Default exclude regular expressions:')
+ for regex in dcmstack.default_key_excl_res:
+ print('\t' + regex)
+ print('Default include regular expressions:')
+ for regex in dcmstack.default_key_incl_res:
+ print('\t' + regex)
+ return 0
+
+ #Check if we are generating meta data
+ gen_meta = args.embed_meta or args.dump_meta
+
+ if gen_meta:
+ #Start with the module defaults
+ ignore_rules = extract.default_ignore_rules
+ translators = extract.default_translators
+
+ #Disable translators if requested
+ if args.disable_translator:
+ if args.disable_translator.lower() == 'all':
+ translators = tuple()
+ else:
+ try:
+ disable_tags = parse_tags(args.disable_translator)
+ except:
+ arg_parser.error('Invalid tag format to --disable-translator.')
+ new_translators = []
+ for translator in translators:
+ if not translator.tag in disable_tags:
+ new_translators.append(translator)
+ translators = new_translators
+
+ #Include non-translated private elements if requested
+ if args.extract_private:
+ ignore_rules = (extract.ignore_pixel_data,
+ extract.ignore_overlay_data,
+ extract.ignore_color_lut_data)
+
+ extractor = extract.MetaExtractor(ignore_rules, translators)
+
+ else:
+ extractor = extract.minimal_extractor
+
+ #Add include/exclude regexes to meta filter
+ include_regexes = dcmstack.default_key_incl_res
+ if args.include_regex:
+ include_regexes += args.include_regex
+ exclude_regexes = dcmstack.default_key_excl_res
+ if args.exclude_regex:
+ exclude_regexes += args.exclude_regex
+ meta_filter = dcmstack.make_key_regex_filter(exclude_regexes,
+ include_regexes)
+
+ #Figure out time and vector ordering
+ if args.time_var:
+ if args.time_order:
+ order_file = open(args.time_order)
+ abs_order = [line.strip() for line in order_file.readlines()]
+ order_file.close()
+ time_order = DicomOrdering(args.time_var, abs_order, True)
+ else:
+ time_order = DicomOrdering(args.time_var)
+ else:
+ time_order = None
+ if args.vector_var:
+ if args.vector_order:
+ order_file = open(args.vector_order)
+ abs_order = [line.strip() for line in order_file.readlines()]
+ order_file.close()
+ vector_order = DicomOrdering(args.vector_var, abs_order, True)
+ else:
+ vector_order = DicomOrdering(args.vector_var)
+ else:
+ vector_order = None
+
+ if len(args.src_dirs) == 0:
+ arg_parser.error('No source directories were provided.')
+
+ #Handle group-by option
+ if not args.group_by is None:
+ group_by = args.group_by.split(',')
+ else:
+ group_by = default_group_keys
+
+ #Handle each source directory individually
+ for src_dir in args.src_dirs:
+ if not os.path.isdir(src_dir):
+ print('%s is not a directory, skipping' % src_dir, file=sys.stderr)
+
+ if args.verbose:
+ print("Processing source directory %s" % src_dir)
+
+ #Build a list of paths to source files
+ glob_str = os.path.join(src_dir, '*')
+ if args.file_ext:
+ glob_str += args.file_ext
+ src_paths = glob(glob_str)
+
+ if args.verbose:
+ print("Found %d source files in the directory" % len(src_paths))
+
+ #Group the files in this directory
+ groups = parse_and_group(src_paths,
+ group_by,
+ extractor,
+ args.force_read,
+ not args.strict,
+ )
+
+ if args.verbose:
+ print("Found %d groups of DICOM images" % len(groups))
+
+ if len(groups) == 0:
+ print("No DICOM files found in %s" % src_dir)
+
+ out_idx = 0
+ generated_outs = set()
+ for key, group in groups.iteritems():
+ stack = stack_group(group,
+ warn_on_except=not args.strict,
+ time_order=time_order,
+ vector_order=vector_order,
+ allow_dummies=args.allow_dummies,
+ meta_filter=meta_filter)
+ meta = group[0][1]
+
+ #Build an appropriate output format string if none was specified
+ if args.output_name is None:
+ out_fmt = []
+ if 'SeriesNumber' in meta:
+ out_fmt.append('%(SeriesNumber)03d')
+ if 'ProtocolName' in meta:
+ out_fmt.append('%(ProtocolName)s')
+ elif 'SeriesDescription' in meta:
+ out_fmt.append('%(SeriesDescription)s')
+ else:
+ out_fmt.append('series')
+ out_fmt = '-'.join(out_fmt)
+ else:
+ out_fmt = args.output_name
+
+ #Get the output filename from the format string, make sure the
+ #result is unique for this source directory
+ out_fn = sanitize_path_comp(out_fmt % meta)
+ if out_fn in generated_outs:
+ out_fn += '-%03d' % out_idx
+ generated_outs.add(out_fn)
+ out_idx += 1
+ out_fn = out_fn + args.output_ext
+
+ if args.dest_dir:
+ out_path = os.path.join(args.dest_dir, out_fn)
+ else:
+ out_path = os.path.join(src_dir, out_fn)
+
+ if args.verbose:
+ print("Writing out stack to path %s" % out_path)
+
+ nii = stack.to_nifti(args.voxel_order, gen_meta)
+
+ if args.dump_meta:
+ nii_wrp = NiftiWrapper(nii)
+ path_tokens = out_path.split('.')
+ if path_tokens[-1] == 'gz':
+ path_tokens = path_tokens[:-1]
+ if path_tokens[-1] == 'nii':
+ path_tokens = path_tokens[:-1]
+ meta_path = '.'.join(path_tokens + ['json'])
+ out_file = open(meta_path, 'w')
+ out_file.write(nii_wrp.meta_ext.to_json())
+ out_file.close()
+
+ if not args.embed_meta:
+ nii_wrp.remove_extension()
+
+ del nii_wrp
+
+ nii.to_filename(out_path)
+
+ del key
+ del group
+ del stack
+ del meta
+ del nii
+ del groups
+
+ return 0
+
+if __name__ == '__main__':
+ sys.exit(main())
diff --git a/fsl/data/dcmstack/extract.py b/fsl/data/dcmstack/extract.py
new file mode 100644
index 000000000..40fd91157
--- /dev/null
+++ b/fsl/data/dcmstack/extract.py
@@ -0,0 +1,520 @@
+"""
+Extract meta data from a DICOM data set.
+"""
+import struct, warnings
+from collections import namedtuple, defaultdict
+import dicom
+from dicom.datadict import keyword_for_tag
+from nibabel.nicom import csareader
+from .dcmstack import DicomStack
+try:
+ from collections import OrderedDict
+except ImportError:
+ from ordereddict import OrderedDict
+try:
+ import chardet
+ have_chardet = True
+except ImportError:
+ have_chardet = False
+ pass
+
+#This is needed to allow extraction on files with invalid values (e.g. too
+#long of a decimal string)
+dicom.config.enforce_valid_values = False
+
+# Python 2 / 3 compatibility
+unicode_str = unicode if sys.version_info[0] < 3 else str
+
+def is_ascii(in_str):
+ '''Return true if the given string is valid ASCII.'''
+ if all(' ' <= c <= '~' for c in in_str):
+ return True
+ return False
+
+def ignore_private(elem):
+ '''Ignore rule for `MetaExtractor` to skip private DICOM elements (odd
+ group number).'''
+ if elem.tag.group % 2 == 1:
+ return True
+ return False
+
+def ignore_pixel_data(elem):
+ return elem.tag == dicom.tag.Tag(0x7fe0, 0x10)
+
+def ignore_overlay_data(elem):
+ return elem.tag.group & 0xff00 == 0x6000 and elem.tag.elem == 0x3000
+
+def ignore_color_lut_data(elem):
+ return (elem.tag.group == 0x28 and
+ elem.tag.elem in (0x1201, 0x1202, 0x1203, 0x1221, 0x1222, 0x1223))
+
+default_ignore_rules = (ignore_private,
+ ignore_pixel_data,
+ ignore_overlay_data,
+ ignore_color_lut_data)
+'''The default tuple of ignore rules for `MetaExtractor`.'''
+
+
+Translator = namedtuple('Translator', ['name',
+ 'tag',
+ 'priv_creator',
+ 'trans_func']
+ )
+'''A namedtuple for storing the four elements of a translator: a name, the
+dicom.tag.Tag that can be translated, the private creator string (optional), and
+the function which takes the DICOM element and returns a dictionary.'''
+
+def simplify_csa_dict(csa_dict):
+ '''Simplify the result of nibabel.nicom.csareader.
+
+ Parameters
+ ----------
+ csa_dict : dict
+ The result from nibabel.nicom.csareader
+
+ Returns
+ -------
+ result : OrderedDict
+ Result where the keys come from the 'tags' sub dictionary of `csa_dict`.
+ The values come from the 'items' within that tags sub sub dictionary.
+ If items has only one element it will be unpacked from the list.
+ '''
+ if csa_dict is None:
+ return None
+
+ result = OrderedDict()
+ for tag in csa_dict['tags']:
+ items = csa_dict['tags'][tag]['items']
+ if len(items) == 0:
+ continue
+ elif len(items) == 1:
+ result[tag] = items[0]
+ else:
+ result[tag] = items
+ return result
+
+def csa_image_trans_func(elem):
+ '''Function for translating the CSA image sub header.'''
+ return simplify_csa_dict(csareader.read(elem.value))
+
+csa_image_trans = Translator('CsaImage',
+ dicom.tag.Tag(0x29, 0x1010),
+ 'SIEMENS CSA HEADER',
+ csa_image_trans_func)
+'''Translator for the CSA image sub header.'''
+
+class PhoenixParseError(Exception):
+ def __init__(self, line):
+ '''Exception indicating a error parsing a line from the Phoenix
+ Protocol.
+ '''
+ self.line = line
+
+ def __str__(self):
+ return 'Unable to parse phoenix protocol line: %s' % self.line
+
+def _parse_phoenix_line(line, str_delim='""'):
+ delim_len = len(str_delim)
+ #Handle most comments (not always when string literal involved)
+ comment_idx = line.find('#')
+ if comment_idx != -1:
+ #Check if the pound sign is in a string literal
+ if line[:comment_idx].count(str_delim) == 1:
+ if line[comment_idx:].find(str_delim) == -1:
+ raise PhoenixParseError(line)
+ else:
+ line = line[:comment_idx]
+
+ #Allow empty lines
+ if line.strip() == '':
+ return None
+
+ #Find the first equals sign and use that to split key/value
+ equals_idx = line.find('=')
+ if equals_idx == -1:
+ raise PhoenixParseError(line)
+ key = line[:equals_idx].strip()
+ val_str = line[equals_idx + 1:].strip()
+
+ #If there is a string literal, pull that out
+ if val_str.startswith(str_delim):
+ end_quote = val_str[delim_len:].find(str_delim) + delim_len
+ if end_quote == -1:
+ raise PhoenixParseError(line)
+ elif not end_quote == len(val_str) - delim_len:
+ #Make sure remainder is just comment
+ if not val_str[end_quote+delim_len:].strip().startswith('#'):
+ raise PhoenixParseError(line)
+
+ return (key, val_str[2:end_quote])
+
+ else: #Otherwise try to convert to an int or float
+ val = None
+ try:
+ val = int(val_str)
+ except ValueError:
+ pass
+ else:
+ return (key, val)
+
+ try:
+ val = int(val_str, 16)
+ except ValueError:
+ pass
+ else:
+ return (key, val)
+
+ try:
+ val = float(val_str)
+ except ValueError:
+ pass
+ else:
+ return (key, val)
+
+ raise PhoenixParseError(line)
+
+def parse_phoenix_prot(prot_key, prot_val):
+ '''Parse the MrPheonixProtocol string.
+
+ Parameters
+ ----------
+ prot_str : str
+ The 'MrPheonixProtocol' string from the CSA Series sub header.
+
+ Returns
+ -------
+ prot_dict : OrderedDict
+ Meta data pulled from the ASCCONV section.
+
+ Raises
+ ------
+ PhoenixParseError : A line of the ASCCONV section could not be parsed.
+ '''
+ if prot_key == 'MrPhoenixProtocol':
+ str_delim = '""'
+ elif prot_key == 'MrProtocol':
+ str_delim = '"'
+ else:
+ raise ValueError('Unknown protocol key: %s' % prot_key)
+ ascconv_start = prot_val.find('### ASCCONV BEGIN ')
+ ascconv_end = prot_val.find('### ASCCONV END ###')
+ ascconv = prot_val[ascconv_start:ascconv_end].split('\n')[1:-1]
+
+ result = OrderedDict()
+ for line in ascconv:
+ parse_result = _parse_phoenix_line(line, str_delim)
+ if parse_result:
+ result[parse_result[0]] = parse_result[1]
+
+ return result
+
+def csa_series_trans_func(elem):
+ '''Function for parsing the CSA series sub header.'''
+ csa_dict = simplify_csa_dict(csareader.read(elem.value))
+
+ #If there is a phoenix protocol, parse it and dump it into the csa_dict
+ phx_src = None
+ if 'MrPhoenixProtocol' in csa_dict:
+ phx_src = 'MrPhoenixProtocol'
+ elif 'MrProtocol' in csa_dict:
+ phx_src = 'MrProtocol'
+
+ if not phx_src is None:
+ phoenix_dict = parse_phoenix_prot(phx_src, csa_dict[phx_src])
+ del csa_dict[phx_src]
+ for key, val in phoenix_dict.items():
+ new_key = '%s.%s' % ('MrPhoenixProtocol', key)
+ csa_dict[new_key] = val
+
+ return csa_dict
+
+csa_series_trans = Translator('CsaSeries',
+ dicom.tag.Tag(0x29, 0x1020),
+ 'SIEMENS CSA HEADER',
+ csa_series_trans_func)
+'''Translator for parsing the CSA series sub header.'''
+
+
+default_translators = (csa_image_trans,
+ csa_series_trans,
+ )
+'''Default translators for MetaExtractor.'''
+
+def tag_to_str(tag):
+ '''Convert a DICOM tag to a string representation using the group and
+ element hex values seprated by an underscore.'''
+ return '%#X_%#X' % (tag.group, tag.elem)
+
+unpack_vr_map = {'SL' : 'i',
+ 'UL' : 'I',
+ 'FL' : 'f',
+ 'FD' : 'd',
+ 'SS' : 'h',
+ 'US' : 'H',
+ 'US or SS' : 'H',
+ }
+'''Dictionary mapping value representations to corresponding format strings for
+the struct.unpack function.'''
+
+def tm_to_seconds(time_str):
+ '''Convert a DICOM time value (value representation of 'TM') to the number
+ of seconds past midnight.
+
+ Parameters
+ ----------
+ time_str : str
+ The DICOM time value string
+
+ Returns
+ -------
+ A floating point representing the number of seconds past midnight
+ '''
+ #Allow ACR/NEMA style format by removing any colon chars
+ time_str = time_str.replace(':', '')
+
+ #Only the hours portion is required
+ result = int(time_str[:2]) * 3600
+
+ str_len = len(time_str)
+ if str_len > 2:
+ result += int(time_str[2:4]) * 60
+ if str_len > 4:
+ result += float(time_str[4:])
+
+ return float(result)
+
+def get_text(byte_str):
+ '''If the given byte string contains text data return it as unicode,
+ otherwise return None.
+
+ If the 'chardet' package is installed, this will be used to detect the
+ text encoding. Otherwise the input will only be decoded if it is ASCII.
+ '''
+ if have_chardet:
+ match = chardet.detect(byte_str)
+ if match['encoding'] is None:
+ return None
+ else:
+ return byte_str.decode(match['encoding'])
+ else:
+ if not is_ascii(byte_str):
+ return None
+ else:
+ return byte_str.decode('ascii')
+
+default_conversions = {'DS' : float,
+ 'IS' : int,
+ 'AT' : str,
+ 'OW' : get_text,
+ 'OB' : get_text,
+ 'OW or OB' : get_text,
+ 'OB or OW' : get_text,
+ 'UN' : get_text,
+ 'PN' : unicode_str,
+ 'UI' : unicode_str,
+ }
+
+class MetaExtractor(object):
+ '''Callable object for extracting meta data from a dicom dataset.
+ Initialize with a set of ignore rules, translators, and type
+ conversions.
+
+ Parameters
+ ----------
+ ignore_rules : sequence
+ A sequence of callables, each of which should take a DICOM element
+ and return True if it should be ignored. If None the module
+ default is used.
+
+ translators : sequence
+ A sequence of `Translator` objects each of which can convert a
+ DICOM element into a dictionary. Overrides any ignore rules. If
+ None the module default is used.
+
+ conversions : dict
+ Mapping of DICOM value representation (VR) strings to callables
+ that perform some conversion on the value
+
+ warn_on_trans_except : bool
+ Convert any exceptions from translators into warnings.
+ '''
+
+ def __init__(self, ignore_rules=None, translators=None, conversions=None,
+ warn_on_trans_except=True):
+ if ignore_rules is None:
+ self.ignore_rules = default_ignore_rules
+ else:
+ self.ignore_rules = ignore_rules
+ if translators is None:
+ self.translators = default_translators
+ else:
+ self.translators = translators
+ if conversions is None:
+ self.conversions = default_conversions
+ else:
+ self.conversions = conversions
+ self.warn_on_trans_except = warn_on_trans_except
+
+ def _get_elem_key(self, elem):
+ '''Get the key for any non-translated elements.'''
+ #Use standard DICOM keywords if possible
+ key = keyword_for_tag(elem.tag)
+
+ #For private tags we take elem.name and convert to camel case
+ if key == '':
+ key = elem.name
+ if key.startswith('[') and key.endswith(']'):
+ key = key[1:-1]
+ tokens = [token[0].upper() + token[1:]
+ for token in key.split()]
+ key = ''.join(tokens)
+
+ return key
+
+ def _get_elem_value(self, elem):
+ '''Get the value for any non-translated elements'''
+ #If the VR is implicit, we may need to unpack the values from a byte
+ #string. This may require us to make an assumption about whether the
+ #value is signed or not, but this is unavoidable.
+ if elem.VR in unpack_vr_map and isinstance(elem.value, str):
+ n_vals = len(elem.value)/struct.calcsize(unpack_vr_map[elem.VR])
+ if n_vals != elem.VM:
+ warnings.warn("The element's VM and the number of values do "
+ "not match.")
+ if n_vals == 1:
+ value = struct.unpack(unpack_vr_map[elem.VR], elem.value)[0]
+ else:
+ value = list(struct.unpack(unpack_vr_map[elem.VR]*n_vals,
+ elem.value)
+ )
+ else:
+ #Otherwise, just take a copy if the value is a list
+ n_vals = elem.VM
+ if n_vals > 1:
+ value = elem.value[:]
+ else:
+ value = elem.value
+
+ #Handle any conversions
+ if elem.VR in self.conversions:
+ if n_vals == 1:
+ value = self.conversions[elem.VR](value)
+ else:
+ value = [self.conversions[elem.VR](val) for val in value]
+
+ return value
+
+ def __call__(self, dcm):
+ '''Extract the meta data from a DICOM dataset.
+
+ Parameters
+ ----------
+ dcm : dicom.dataset.Dataset
+ The DICOM dataset to extract the meta data from.
+
+ Returns
+ -------
+ meta : dict
+ A dictionary of extracted meta data.
+
+ Notes
+ -----
+ Non-private tags use the DICOM keywords as keys. Translators have their
+ name, followed by a dot, prepended to the keys of any meta elements
+ they produce. Values are unchanged, except when the value
+ representation is 'DS' or 'IS' (decimal/integer strings) they are
+ converted to float and int types.
+ '''
+ standard_meta = []
+ trans_meta_dicts = OrderedDict()
+
+ #Make dict to track which tags map to which translators
+ trans_map = {}
+
+ # Convert text elements to unicode
+ dcm.decode()
+
+ for elem in dcm:
+ if isinstance(elem.value, str) and elem.value.strip() == '':
+ continue
+
+ #Get the name for non-translated elements
+ name = self._get_elem_key(elem)
+
+ #If it is a private creator element, setup any corresponding
+ #translators
+ if elem.name == "Private Creator":
+ for translator in self.translators:
+ if translator.priv_creator == elem.value:
+ new_elem = ((translator.tag.elem & 0xff) |
+ (elem.tag.elem * 16**2))
+ new_tag = dicom.tag.Tag(elem.tag.group, new_elem)
+ if new_tag in trans_map:
+ raise ValueError('More than one translator '
+ 'for tag: %s' % new_tag)
+ trans_map[new_tag] = translator
+
+ #If there is a translator for this element, use it
+ if elem.tag in trans_map:
+ try:
+ meta = trans_map[elem.tag].trans_func(elem)
+ except Exception as e:
+ if self.warn_on_trans_except:
+ warnings.warn("Exception from translator %s: %s" %
+ (trans_map[elem.tag].name,
+ repr(str(e))))
+ else:
+ raise
+ else:
+ if meta:
+ trans_meta_dicts[trans_map[elem.tag].name] = meta
+ #Otherwise see if we are supposed to ignore the element
+ elif any(rule(elem) for rule in self.ignore_rules):
+ continue
+ #Handle elements that are sequences with recursion
+ elif isinstance(elem.value, dicom.sequence.Sequence):
+ value = []
+ for val in elem.value:
+ value.append(self(val))
+ if all(x is None for x in value):
+ continue
+ standard_meta.append((name, value, elem.tag))
+ #Otherwise just make sure the value is unpacked
+ else:
+ value = self._get_elem_value(elem)
+ if value is None:
+ continue
+ standard_meta.append((name, value, elem.tag))
+
+ #Handle name collisions
+ name_counts = defaultdict(int)
+ for elem in standard_meta:
+ name_counts[elem[0]] += 1
+ result = OrderedDict()
+ for name, value, tag in standard_meta:
+ if name_counts[name] > 1:
+ name = name + '_' + tag_to_str(tag)
+ result[name] = value
+
+ #Inject translator results
+ for trans_name, meta in trans_meta_dicts.items():
+ for name, value in meta.items():
+ name = '%s.%s' % (trans_name, name)
+ result[name] = value
+
+ return result
+
+def minimal_extractor(dcm):
+ '''Meta data extractor that just extracts the minimal set of keys needed
+ by DicomStack objects.
+ '''
+ result = {}
+ for key in DicomStack.minimal_keys:
+ try:
+ result[key] = dcm.__getattr__(key)
+ except AttributeError:
+ pass
+ return result
+
+default_extractor = MetaExtractor()
+'''The default `MetaExtractor`.'''
diff --git a/fsl/data/dcmstack/info.py b/fsl/data/dcmstack/info.py
new file mode 100644
index 000000000..a39ab07fd
--- /dev/null
+++ b/fsl/data/dcmstack/info.py
@@ -0,0 +1,54 @@
+""" Information for setup.py that we may also want to access in dcmstack. Can
+not import dcmstack.
+"""
+import sys
+
+_version_major = 0
+_version_minor = 7
+_version_micro = 0
+_version_extra = 'dev'
+__version__ = "%s.%s.%s%s" % (_version_major,
+ _version_minor,
+ _version_micro,
+ _version_extra)
+
+CLASSIFIERS = ["Development Status :: 3 - Alpha",
+ "Environment :: Console",
+ "Intended Audience :: Science/Research",
+ "License :: OSI Approved :: MIT License",
+ "Operating System :: OS Independent",
+ "Programming Language :: Python",
+ "Topic :: Scientific/Engineering"]
+
+description = 'Stack DICOM images into volumes and convert to Nifti'
+
+# Hard dependencies
+install_requires = ['pydicom >= 0.9.7',
+ 'nibabel >= 2.0.0',
+ ]
+# Add version specific dependencies
+if sys.version_info < (2, 6):
+ raise Exception("must use python 2.6 or greater")
+elif sys.version_info < (2, 7):
+ install_requires.append('ordereddict')
+
+# Extra requirements for building documentation and testing
+extras_requires = {'doc': ["sphinx", "numpydoc"],
+ 'test': ["nose"],
+ }
+
+
+NAME = 'dcmstack'
+AUTHOR = "Brendan Moloney"
+AUTHOR_EMAIL = "moloney@ohsu.edu"
+MAINTAINER = "Brendan Moloney"
+MAINTAINER_EMAIL = "moloney@ohsu.edu"
+DESCRIPTION = description
+LICENSE = "MIT license"
+CLASSIFIERS = CLASSIFIERS
+PLATFORMS = "OS Independent"
+ISRELEASE = _version_extra == ''
+VERSION = __version__
+INSTALL_REQUIRES = install_requires
+EXTRAS_REQUIRES = extras_requires
+PROVIDES = ["dcmstack"]
\ No newline at end of file
diff --git a/fsl/data/dcmstack/nitool_cli.py b/fsl/data/dcmstack/nitool_cli.py
new file mode 100644
index 000000000..5e83da519
--- /dev/null
+++ b/fsl/data/dcmstack/nitool_cli.py
@@ -0,0 +1,251 @@
+"""
+Command line interface for nitool.
+
+@author: moloney
+"""
+from __future__ import print_function
+
+import os, sys, argparse
+import nibabel as nb
+from .dcmmeta import NiftiWrapper, DcmMetaExtension, MissingExtensionError
+
+prog_descrip = """Work with extended Nifti files created by dcmstack"""
+
+def main(argv=sys.argv):
+ #Setup the top level parser
+ arg_parser = argparse.ArgumentParser(description=prog_descrip)
+ sub_parsers = arg_parser.add_subparsers(title="Subcommands")
+
+ #Split command
+ split_help = ("Split src_nii file along a dimension. Defaults to the slice "
+ "dimension if 3D, otherwise the last dimension.")
+ split_parser = sub_parsers.add_parser('split', help=split_help)
+ split_parser.add_argument('src_nii', nargs=1)
+ split_parser.add_argument('-d', '--dimension', default=None, type=int,
+ help=("The dimension to split along. Must be in "
+ "the range [0, 5)"))
+ split_parser.add_argument('-o', '--output-format', default=None,
+ help=("Format string used to create the output "
+ "file names. Default is to prepend the index "
+ "number to the src_nii filename."))
+ split_parser.set_defaults(func=split)
+
+ #Merge Command
+ merge_help = ("Merge the provided Nifti files along a dimension. Defaults "
+ "to slice, then time, and then vector.")
+ merge_parser = sub_parsers.add_parser('merge', help=merge_help)
+ merge_parser.add_argument('output', nargs=1)
+ merge_parser.add_argument('src_niis', nargs='+')
+ merge_parser.add_argument('-d', '--dimension', default=None, type=int,
+ help=("The dimension to merge along. Must be "
+ "in the range [0, 5)"))
+ merge_parser.add_argument('-s', '--sort', default=None,
+ help=("Sort the source files using the provided "
+ "meta data key before merging"))
+ merge_parser.add_argument('-c', '--clear-slices', action='store_true',
+ help="Clear all per slice meta data")
+ merge_parser.set_defaults(func=merge)
+
+ #Dump Command
+ dump_help = "Dump the JSON meta data extension from the provided Nifti."
+ dump_parser = sub_parsers.add_parser('dump', help=dump_help)
+ dump_parser.add_argument('src_nii', nargs=1)
+ dump_parser.add_argument('dest_json', nargs='?',
+ type=argparse.FileType('w'),
+ default=sys.stdout)
+ dump_parser.add_argument('-m', '--make-empty', default=False,
+ action='store_true',
+ help="Make an empty extension if none exists")
+ dump_parser.add_argument('-r', '--remove', default=False,
+ action='store_true',
+ help="Remove the extension from the Nifti file")
+ dump_parser.set_defaults(func=dump)
+
+ #Embed Command
+ embed_help = "Embed a JSON extension into the Nifti file."
+ embed_parser = sub_parsers.add_parser('embed', help=embed_help)
+ embed_parser.add_argument('src_json', nargs='?', type=argparse.FileType('r'),
+ default=sys.stdin)
+ embed_parser.add_argument('dest_nii', nargs=1)
+ embed_parser.add_argument('-f', '--force-overwrite', action='store_true',
+ help="Overwrite any existing dcmmeta extension")
+ embed_parser.set_defaults(func=embed)
+
+ #Lookup command
+ lookup_help = "Lookup the value for the given meta data key."
+ lookup_parser = sub_parsers.add_parser('lookup', help=lookup_help)
+ lookup_parser.add_argument('key', nargs=1)
+ lookup_parser.add_argument('src_nii', nargs=1)
+ lookup_parser.add_argument('-i', '--index',
+ help=("Use the given voxel index. The index "
+ "must be provided as a comma separated list of "
+ "integers (one for each dimension)."))
+ lookup_parser.set_defaults(func=lookup)
+
+ #Inject command
+ inject_help = "Inject meta data into the JSON extension."
+ inject_parser = sub_parsers.add_parser('inject', help=inject_help)
+ inject_parser.add_argument('dest_nii', nargs=1)
+ inject_parser.add_argument('classification', nargs=2)
+ inject_parser.add_argument('key', nargs=1)
+ inject_parser.add_argument('values', nargs='+')
+ inject_parser.add_argument('-f', '--force-overwrite',
+ action='store_true',
+ help=("Overwrite any existing values "
+ "for the key"))
+ inject_parser.add_argument('-t', '--type', default=None,
+ help="Interpret the value as this type instead "
+ "of trying to determine the type automatically")
+ inject_parser.set_defaults(func=inject)
+
+ #Parse the arguments and call the appropriate function
+ args = arg_parser.parse_args(argv[1:])
+ return args.func(args)
+
+def split(args):
+ src_path = args.src_nii[0]
+ src_fn = os.path.basename(src_path)
+ src_dir = os.path.dirname(src_path)
+
+ src_nii = nb.load(src_path)
+ try:
+ src_wrp = NiftiWrapper(src_nii)
+ except MissingExtensionError:
+ print("No dcmmeta extension found, making empty one...")
+ src_wrp = NiftiWrapper(src_nii, make_empty=True)
+ for split_idx, split in enumerate(src_wrp.split(args.dimension)):
+ if args.output_format:
+ out_name = (args.output_format %
+ split.meta_ext.get_class_dict(('global', 'const'))
+ )
+ else:
+ out_name = os.path.join(src_dir, '%03d-%s' % (split_idx, src_fn))
+ nb.save(split, out_name)
+ return 0
+
+def make_key_func(meta_key, index=None):
+ def key_func(src_nii):
+ result = src_nii.get_meta(meta_key, index)
+ if result is None:
+ raise ValueError('Key not found: %s' ) % meta_key
+ return result
+
+ return key_func
+
+def merge(args):
+ src_wrps = []
+ for src_path in args.src_niis:
+ src_nii = nb.load(src_path)
+ try:
+ src_wrp = NiftiWrapper(src_nii)
+ except MissingExtensionError:
+ print("No dcmmeta extension found, making empty one...")
+ src_wrp = NiftiWrapper(src_nii, make_empty=True)
+ src_wrps.append(src_wrp)
+
+ if args.sort:
+ src_wrps.sort(key=make_key_func(args.sort))
+
+ result_wrp = NiftiWrapper.from_sequence(src_wrps, args.dimension)
+
+ if args.clear_slices:
+ result_wrp.meta_ext.clear_slice_meta()
+
+ out_name = (args.output[0] %
+ result_wrp.meta_ext.get_class_dict(('global', 'const')))
+ result_wrp.to_filename(out_name)
+ return 0
+
+def dump(args):
+ src_nii = nb.load(args.src_nii[0])
+ src_wrp = NiftiWrapper(src_nii, args.make_empty)
+ meta_str = src_wrp.meta_ext.to_json()
+ args.dest_json.write(meta_str)
+ args.dest_json.write('\n')
+
+ if args.remove:
+ src_wrp.remove_extension()
+ src_wrp.to_filename(args.src_nii[0])
+ return 0
+
+def check_overwrite():
+ usr_input = ''
+ while not usr_input in ('y', 'n'):
+ usr_input = input('Existing DcmMeta extension found, overwrite? '
+ '[y/n]').lower()
+ return usr_input == 'y'
+
+def embed(args):
+ dest_nii = nb.load(args.dest_nii[0])
+ hdr = dest_nii.get_header()
+ try:
+ src_wrp = NiftiWrapper(dest_nii, False)
+ except MissingExtensionError:
+ pass
+ else:
+ if not args.force_overwrite:
+ if not check_overwrite():
+ return
+ src_wrp.remove_extension()
+
+ hdr.extensions.append(DcmMetaExtension.from_json(args.src_json.read()))
+ nb.save(dest_nii, args.dest_nii[0])
+ return 0
+
+def lookup(args):
+ src_wrp = NiftiWrapper.from_filename(args.src_nii[0])
+ index = None
+ if args.index:
+ index = tuple(int(idx.strip()) for idx in args.index.split(','))
+ meta = src_wrp.get_meta(args.key[0], index)
+ if not meta is None:
+ print(meta)
+ return 0
+
+def convert_values(values, type_str=None):
+ if type_str is None:
+ for conv_type in (int, float):
+ try:
+ values = [conv_type(val) for val in values]
+ except ValueError:
+ pass
+ else:
+ break
+ else:
+ if type_str not in ('str', 'int', 'float'):
+ raise ValueError("Unrecognized type: %s" % type_str)
+ conv_type = eval(type_str)
+ values = [conv_type(val) for val in values]
+ if len(values) == 1:
+ return values[0]
+ return values
+
+def inject(args):
+ dest_nii = nb.load(args.dest_nii[0])
+ dest_wrp = NiftiWrapper(dest_nii, make_empty=True)
+ classification = tuple(args.classification)
+ if not classification in dest_wrp.meta_ext.get_valid_classes():
+ print("Invalid classification: %s" % (classification,))
+ return 1
+ n_vals = len(args.values)
+ mult = dest_wrp.meta_ext.get_multiplicity(classification)
+ if n_vals != mult:
+ print(("Invalid number of values for classification. Expected "
+ "%d but got %d") % (mult, n_vals))
+ return 1
+ key = args.key[0]
+ if key in dest_wrp.meta_ext.get_keys():
+ if not args.force_overwrite:
+ print("Key already exists, must pass --force-overwrite")
+ return 1
+ else:
+ curr_class = dest_wrp.meta_ext.get_classification(key)
+ curr_dict = dest_wrp.meta_ext.get_class_dict(curr_class)
+ del curr_dict[key]
+ class_dict = dest_wrp.meta_ext.get_class_dict(classification)
+ class_dict[key] = convert_values(args.values, args.type)
+ nb.save(dest_nii, args.dest_nii[0])
+ return 0
+
+if __name__ == '__main__':
+ sys.exit(main())
--
GitLab