From 64ed1882457e6867305a4ba010c4b53c1f705f32 Mon Sep 17 00:00:00 2001 From: Gregory Lee Date: Thu, 18 Oct 2018 11:41:53 -0400 Subject: [PATCH 01/16] add Philips XML/REC reader. Philips now encourages the XML-style headers over the older PAR format. This initial implementation reads the XML file, converting each field to its PAR file equivalent and then inherits from the PARRECHeader and PARRECImage to avoid excessive code duplication. --- nibabel/xmlrec.py | 584 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 584 insertions(+) create mode 100755 nibabel/xmlrec.py diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py new file mode 100755 index 0000000000..9ed9aedf67 --- /dev/null +++ b/nibabel/xmlrec.py @@ -0,0 +1,584 @@ +from copy import copy, deepcopy +import warnings +import xml.etree.ElementTree as ET + +import numpy as np +from nibabel.parrec import one_line +from nibabel.parrec import PARRECImage, PARRECHeader + + +class XMLRECError(Exception): + """Exception for XML/REC format related problems. + + To be raised whenever XML/REC is not happy, or we are not happy with + XML/REC. + """ + +# Dictionary of conversions from enumerated types to integer values. +# The following named Enums are used in XML files. +# I do not know thefull definitions for these. +# This was partially adapted from: +# https://github.com/xiangruili/dicm2nii/blob/6855498a702b06bbc7d1de716f4e9bdfe50c74fc/xml2par.m#L356 +enums_dict = { + 'Label Type': {'-': 1}, + 'Type': {'M': 0, 'R': 1, 'I': 2, 'P': 3}, # TODO: 'B0', etc... + 'Sequence': {'SE': 1, 'FFE': 2}, + 'Image Type Ed Es': {'U': 2}, # TODO: others? + 'Display Orientation': {'-': 0, 'NONE': 0}, # TODO: correct value for None + 'Slice Orientation': {'Transversal': 1, 'Sagittal': 2, 'Coronal': 3}, + 'Contrast Type': {'DIFFUSION': 0, 'T2': 4, 'TAGGING': 6, 'T1': 7}, + 'Diffusion Anisotropy Type': {'-': 0}} + +# Dict for converting XML/REC types to appropriate python types +# could convert the Enumeration strings to int, but enums_dict is incomplete +# and the strings are more easily interpretable +xml_type_dict = {'Float': np.float32, + 'Double': np.float64, + 'Int16': np.int16, + 'Int32': np.int32, + 'UInt16': np.uint16, + 'Enumeration': str, + 'Boolean': bool, + 'String': str, + 'Date': str, + 'Time': str} + +# for use with parrec.py, the Enumeration strings must be converted to ints +par_type_dict = copy(xml_type_dict) +par_type_dict['Enumeration'] = int + +supported_xml_versions = ['PRIDE_V5'] + +# convert XML names to expected property names as used in nibabel.parrec +# Note: additional special cases are handled in general_info_xml_to_par() +general_info_XML_to_nibabel = { + 'Patient Name': 'patient_name', + 'Examination Name': 'exam_name', + 'Protocol Name': 'protocol_name', + 'Aquisition Number': 'acq_nr', + 'Reconstruction Number': 'recon_nr', + 'Scan Duration': 'scan_duration', + 'Max No Phases': 'max_cardiac_phases', + 'Max No Echoes': 'max_echoes', + 'Max No Slices': 'max_slices', + 'Max No Dynamics': 'max_dynamics', + 'Max No Mixes': 'max_mixes', + 'Patient Position': 'patient_position', + 'Preparation Direction': 'prep_direction', + 'Technique': 'tech', + 'Scan Mode': 'scan_mode', + 'Repetition Times': 'repetition_time', + 'Water Fat Shift': 'water_fat_shift', + 'Flow Compensation': 'flow_compensation', + 'Presaturation': 'presaturation', + 'Phase Encoding Velocity': 'phase_enc_velocity', + 'MTC': 'mtc', + 'SPIR': 'spir', + 'EPI factor': 'epi_factor', + 'Dynamic Scan': 'dyn_scan', + 'Diffusion': 'diffusion', + 'Diffusion Echo Time': 'diffusion_echo_time', + 'Max No B Values': 'max_diffusion_values', + 'Max No Gradient Orients': 'max_gradient_orient', + 'No Label Types': 'nr_label_types', + 'Series Data Type': 'series_type'} + + +def general_info_xml_to_par(xml_info): + """Convert general_info from XML-style names to PAR-style names.""" + xml_info_init = xml_info + xml_info = deepcopy(xml_info) + general_info = {} + for k in xml_info_init.keys(): + # convert all keys with a simple 1-1 name conversion + if k in general_info_XML_to_nibabel: + general_info[general_info_XML_to_nibabel[k]] = xml_info.pop(k) + try: + general_info['exam_date'] = '{} / {}'.format( + xml_info.pop('Examination Date'), + xml_info.pop('Examination Time')) + except KeyError: + pass + + try: + general_info['angulation'] = np.asarray( + [xml_info.pop('Angulation AP'), + xml_info.pop('Angulation FH'), + xml_info.pop('Angulation RL')]) + except KeyError: + pass + + try: + general_info['off_center'] = np.asarray( + [xml_info.pop('Off Center AP'), + xml_info.pop('Off Center FH'), + xml_info.pop('Off Center RL')]) + except KeyError: + pass + + try: + general_info['fov'] = np.asarray( + [xml_info.pop('FOV AP'), + xml_info.pop('FOV FH'), + xml_info.pop('FOV RL')]) + except KeyError: + pass + + try: + general_info['scan_resolution'] = np.asarray( + [xml_info.pop('Scan Resolution X'), + xml_info.pop('Scan Resolution Y')], + dtype=int) + except KeyError: + pass + + # copy any excess keys not normally in the .PARREC general info + # These will not be used by the PARREC code, but are kept for completeness + general_info.update(xml_info) + + return general_info + + +# TODO: remove this function? It is currently unused, but may be useful for +# testing roundtrip convesion. +def general_info_par_to_xml(par_info): + """Convert general_info from PAR-style names to XML-style names.""" + general_info_nibabel_to_XML = { + v: k for k, v in general_info_XML_to_nibabel.items()} + par_info_init = par_info + par_info = deepcopy(par_info) + general_info = {} + for k in par_info_init.keys(): + # convert all keys with a simple 1-1 name conversion + if k in general_info_nibabel_to_XML: + general_info[general_info_nibabel_to_XML[k]] = par_info.pop(k) + try: + tmp = par_info['exam_date'].split('/') + general_info['Examination Date'] = tmp[0].strip() + general_info['Examination Time'] = tmp[1].strip() + except KeyError: + pass + + try: + general_info['Angulation AP'] = par_info['angulation'][0] + general_info['Angulation FH'] = par_info['angulation'][1] + general_info['Angulation RL'] = par_info['angulation'][2] + par_info.pop('angulation') + except KeyError: + pass + + try: + general_info['Off Center AP'] = par_info['off_center'][0] + general_info['Off Center FH'] = par_info['off_center'][1] + general_info['Off Center RL'] = par_info['off_center'][2] + par_info.pop('off_center') + except KeyError: + pass + + try: + general_info['FOV AP'] = par_info['fov'][0] + general_info['FOV FH'] = par_info['fov'][1] + general_info['FOV RL'] = par_info['fov'][2] + par_info.pop('fov') + except KeyError: + pass + + try: + general_info['Scan Resolution X'] = par_info['scan_resolution'][0] + general_info['Scan Resolution Y'] = par_info['scan_resolution'][1] + par_info.pop('scan_resolution') + except KeyError: + pass + + # copy any unrecognized keys as is. + # known keys found in XML PRIDE_V5, but not in PAR v4.2 are: + # 'Samples Per Pixel' and 'Image Planar Configuration' + general_info.update(par_info) + + return general_info + +# dictionary mapping fieldnames in the XML header to the corresponding names +# in a PAR V4.2 file +image_def_XML_to_PAR = { + 'Image Type Ed Es': 'image_type_ed_es', + 'No Averages': 'number of averages', + 'Max RR Interval': 'maximum RR-interval', + 'Type': 'image_type_mr', + 'Display Orientation': 'image_display_orientation', + 'Image Flip Angle': 'image_flip_angle', + 'Rescale Slope': 'rescale slope', + 'Label Type': 'label type', + 'fMRI Status Indication': 'fmri_status_indication', + 'TURBO Factor': 'TURBO factor', + 'Min RR Interval': 'minimum RR-interval', + 'Scale Slope': 'scale slope', + 'Inversion Delay': 'Inversion delay', + 'Window Width': 'window width', + 'Sequence': 'scanning sequence', + 'Diffusion Anisotropy Type': 'diffusion anisotropy type', + 'Index': 'index in REC file', + 'Rescale Intercept': 'rescale intercept', + 'Diffusion B Factor': 'diffusion_b_factor', + 'Trigger Time': 'trigger_time', + 'Echo': 'echo number', + 'Echo Time': 'echo_time', + 'Pixel Spacing': 'pixel spacing', + 'Slice Gap': 'slice gap', + 'Dyn Scan Begin Time': 'dyn_scan_begin_time', + 'Window Center': 'window center', + 'Contrast Type': 'contrast type', + 'Slice': 'slice number', + 'BValue': 'diffusion b value number', + 'Scan Percentage': 'scan percentage', + 'Phase': 'cardiac phase number', + 'Slice Thickness': 'slice thickness', + 'Slice Orientation': 'slice orientation', + 'Dynamic': 'dynamic scan number', + 'Pixel Size': 'image pixel size', + 'Grad Orient': 'gradient orientation number', + 'Cardiac Frequency': 'cardiac frequency'} + +# copy of enums_dict but with key names converted to their PAR equivalents +enums_dict_PAR = {image_def_XML_to_PAR[k]: v for k, v in enums_dict.items()} + + +# TODO?: The following values have different names in the XML vs. PAR header +rename_XML_to_PAR = { + 'HFS': 'Head First Supine', + 'LR': 'Left-Right', + 'RL': 'Right-Left', + 'AP': 'Anterior-Posterior', + 'PA': 'Posterior-Anterior', + 'N': 0, + 'Y': 1, +} + + +def _process_gen_dict_XML(xml_root): + """Read the general_info from an XML file. + + This is the equivalent of _process_gen_dict() for .PAR files + """ + info = xml_root.find('Series_Info') + if info is None: + raise RuntimeError("No 'Series_Info' found in the XML file") + general_info = {} + for e in info: + a = e.attrib + if 'Name' in a: + entry_type = xml_type_dict[a['Type']] + if 'ArraySize' in a: + val = [entry_type(i) for i in e.text.strip().split()] + else: + val = entry_type(e.text) + general_info[a['Name']] = val + return general_info + + +def _get_image_def_attributes(xml_root, dtype_format='xml'): + """Get names and dtypes for all attributes defined for each image. + + called by _process_image_lines_xml + + Paramters + --------- + xml_root : + TODO + dtype_format : {'xml', 'par'} + If 'par'', XML paramter names and dtypes are converted to their PAR + equivalents. + + Returns + ------- + key_attributes : list of tuples + The attributes that are used when sorting data from the REC file. + other_attributes : list of tuples + Additional attributes that are not considered when sorting. + """ + image_defs_array = xml_root.find('Image_Array') + if image_defs_array is None: + raise XMLRECError("No 'Image_Array' found in the XML file") + + # can get all of the fields from the first entry + first_def = image_defs_array[0] + + # Key element contains attributes corresponding to image keys + # e.g. (Slice, Echo, Dynamic, ...) that can be used for sorting the images. + img_keys = first_def.find('Key') + if img_keys is None: + raise XMLRECError( + "Expected to find a Key attribute for each element in Image_Array") + if not np.all(['Name' in k.attrib for k in img_keys]): + raise XMLRECError("Expected each Key attribute to have a Name") + + dtype_format = dtype_format.lower() + if dtype_format == 'xml': + # xml_type_dict keeps enums as their string representation + key_attributes = [ + (k.get('Name'), xml_type_dict[k.get('Type')]) for k in img_keys] + elif dtype_format == 'par': + # par_type_dict converts enums to int as used in PAR/REC + key_attributes = [ + (image_def_XML_to_PAR[k.get('Name')], par_type_dict[k.get('Type')]) + for k in img_keys] + else: + raise XMLRECError("dtype_format must be 'par' or 'xml'") + + # Process other attributes that are not considered image keys + other_attributes = [] + for element in first_def: + a = element.attrib + if 'Name' in a: + # if a['Type'] == 'Enumeration': + # enum_type = a['EnumType'] + # print("enum_type = {}".format(enum_type)) + name = a['Name'] + if dtype_format == 'xml': + entry_type = xml_type_dict[a['Type']] + else: + entry_type = par_type_dict[a['Type']] + if name in image_def_XML_to_PAR: + name = image_def_XML_to_PAR[name] + if 'ArraySize' in a: + # handle vector entries (e.g. 'Pixel Size' is length 2) + entry = (name, entry_type, int(a['ArraySize'])) + else: + entry = (name, entry_type) + other_attributes.append(entry) + + return key_attributes, other_attributes + + +# these keys are elements of a multi-valued key in the PAR/REC image_defs +composite_PAR_keys = ['Diffusion AP', 'Diffusion FH', 'Diffusion RL', + 'Angulation AP', 'Angulation FH', 'Angulation RL', + 'Offcenter AP', 'Offcenter FH', 'Offcenter RL', + 'Resolution X', 'Resolution Y'] + + +def _composite_attributes_xml_to_par(other_attributes): + """utility used in conversion from XML-style to PAR-style image_defs. + + called by _process_image_lines_xml + """ + if ('Diffusion AP', np.float32) in other_attributes: + other_attributes.remove(('Diffusion AP', np.float32)) + other_attributes.remove(('Diffusion FH', np.float32)) + other_attributes.remove(('Diffusion RL', np.float32)) + other_attributes.append(('diffusion', float, (3, ))) + + if ('Angulation AP', np.float64) in other_attributes: + other_attributes.remove(('Angulation AP', np.float64)) + other_attributes.remove(('Angulation FH', np.float64)) + other_attributes.remove(('Angulation RL', np.float64)) + other_attributes.append(('image angulation', float, (3, ))) + + if ('Offcenter AP', np.float64) in other_attributes: + other_attributes.remove(('Offcenter AP', np.float64)) + other_attributes.remove(('Offcenter FH', np.float64)) + other_attributes.remove(('Offcenter RL', np.float64)) + other_attributes.append(('image offcentre', float, (3, ))) + + if ('Resolution X', np.uint16) in other_attributes: + other_attributes.remove(('Resolution X', np.uint16)) + other_attributes.remove(('Resolution Y', np.uint16)) + other_attributes.append(('recon resolution', int, (2, ))) + + return other_attributes + + +# TODO: remove? this one is currently unused, but perhaps useful for testing +def _composite_attributes_par_to_xml(other_attributes): + if ('diffusion', float, (3, )) in other_attributes: + other_attributes.append(('Diffusion AP', np.float32)) + other_attributes.append(('Diffusion FH', np.float32)) + other_attributes.append(('Diffusion RL', np.float32)) + other_attributes.remove(('diffusion', float, (3, ))) + + if ('image angulation', float, (3, )) in other_attributes: + other_attributes.append(('Angulation AP', np.float64)) + other_attributes.append(('Angulation FH', np.float64)) + other_attributes.append(('Angulation RL', np.float64)) + other_attributes.remove(('image angulation', float, (3, ))) + + if ('image offcentre', float, (3, )) in other_attributes: + other_attributes.append(('Offcenter AP', np.float64)) + other_attributes.append(('Offcenter FH', np.float64)) + other_attributes.append(('Offcenter RL', np.float64)) + other_attributes.remove(('image offcentre', float, (3, ))) + + if ('recon resolution', int, (2, )) in other_attributes: + other_attributes.append(('Resolution X', np.uint16)) + other_attributes.append(('Resolution Y', np.uint16)) + other_attributes.remove(('recon resolution', int, (2, ))) + return other_attributes + + +def _get_composite_key_index(key): + """utility used in conversion from XML-style to PAR-style image_defs. + + called by _process_image_lines_xml + """ + if 'Diffusion' in key: + name = 'diffusion' + elif 'Angulation' in key: + name = 'image angulation' + elif 'Offcenter' in key: + name = 'image offcentre' + elif 'Resolution' in key: + name = 'recon resolution' + else: + raise ValueError("unrecognized composite element name: {}".format(key)) + + if key in ['Diffusion AP', 'Angulation AP', 'Offcenter AP', + 'Resolution X']: + index = 0 + elif key in ['Diffusion FH', 'Angulation FH', 'Offcenter FH', + 'Resolution Y']: + index = 1 + elif key in ['Diffusion RL', 'Angulation RL', 'Offcenter RL']: + index = 2 + else: + raise ValueError("unrecognized composite element name: {}".format(key)) + return (name, index) + + +def _process_image_lines_xml(xml_root, dtype_format='xml'): + """Build image_defs by parsing the XML file. + + Parameters + ---------- + xml_root : + TODO + dtype_format : {'xml', 'par'} + If 'par'', XML paramter names and dtypes are converted to their PAR file + equivalents. + """ + image_defs_array = xml_root.find('Image_Array') + if image_defs_array is None: + raise RuntimeError("No 'Image_Array' found in the XML file") + + key_attributes, other_attributes = _get_image_def_attributes(xml_root, dtype_format=dtype_format) + + image_def_dtd = key_attributes + other_attributes + # dtype dict based on the XML attribute names + dtype_dict = {a[0]: a[1] for a in image_def_dtd} + + if dtype_format == 'par': + # image_defs based on conversion of types in composite_PAR_keys + image_def_dtd = _composite_attributes_xml_to_par(image_def_dtd) + + image_defs = np.zeros(len(image_defs_array), dtype=image_def_dtd) + for i, image_def in enumerate(image_defs_array): + + if image_def.find('Key') != image_def[0]: + raise RuntimeError("Expected first element of image_def to be Key") + + key_def = image_def[0] + for key in key_def: + name = key.get('Name') + if name in image_def_XML_to_PAR: + name = image_def_XML_to_PAR[name] + if dtype_format == 'par' and name in enums_dict_PAR: + # string -> int + val = enums_dict_PAR[name][key.text] + else: + val = key.text + image_defs[name][i] = dtype_dict[name](val) + + # for all image properties we know about + for element in image_def[1:]: + a = element.attrib + text = element.text + if 'Name' in a: + name = a['Name'] + if dtype_format == 'par' and name in image_def_XML_to_PAR: + name = image_def_XML_to_PAR[name] + # composite_PAR_keys + entry_dtype = dtype_dict[name] + if 'ArraySize' in a: + val = [entry_dtype(i) for i in text.strip().split()] + else: + if dtype_format == 'par' and name in enums_dict_PAR: + # string -> int + val = enums_dict_PAR[name][text] + else: + val = entry_dtype(text) + if dtype_format == 'par' and name in composite_PAR_keys: + # conversion of types in composite_PAR_keys + name, vec_index = _get_composite_key_index(name) + image_defs[name][i][vec_index] = val + else: + image_defs[name][i] = val + return image_defs + + +def parse_XML_header(fobj, dtype_format='xml'): + """Parse a XML header and aggregate all information into useful containers. + + Parameters + ---------- + fobj : file-object + The PAR header file object. + dtype_format : {'xml', 'par'} + If 'par' the image_defs will be converted to a format matching that + found in PARRECHeader. + + Returns + ------- + general_info : dict + Contains all "General Information" from the header file + image_info : ndarray + Structured array with fields giving all "Image information" in the + header + """ + # single pass through the header + tree = ET.parse(fobj) + root = tree.getroot() + + # _split_header() equivalent + + version = root.tag # e.g. PRIDE_V5 + if version not in supported_xml_versions: + warnings.warn(one_line( + """XML version '{0}' is currently not supported. Only PRIDE_V5 XML + files have been tested. --making an attempt to read nevertheless. + Please email the NiBabel mailing list, if you are interested in + adding support for this version. + """.format(version))) + general_info = _process_gen_dict_XML(root) + image_defs = _process_image_lines_xml( + root, dtype_format=dtype_format) + + return general_info, image_defs + + +class XMLRECHeader(PARRECHeader): + @classmethod + def from_fileobj(klass, fileobj, permit_truncated=False, + strict_sort=False): + info, image_defs = parse_XML_header(fileobj, dtype_format='par') + # convert to PAR/REC format general_info + info = general_info_xml_to_par(info) + return klass(info, image_defs, permit_truncated, strict_sort) + @classmethod + def from_header(klass, header=None): + if header is None: + raise XMLRECError('Cannot create XMLRECHeader from air.') + if type(header) == klass: + return header.copy() + raise XMLRECError('Cannot create XMLREC header from ' + 'non-XMLREC header.') + def copy(self): + return XMLRECHeader(deepcopy(self.general_info), + self.image_defs.copy(), + self.permit_truncated, + self.strict_sort) + +class XMLRECImage(PARRECImage): + """XML/REC image""" + header_class = XMLRECHeader + valid_exts = ('.REC', '.xml') + files_types = (('image', '.REC'), ('header', '.xml')) + +load = XMLRECImage.load From 1ff1826ae45254b5b867675fa2bbe8e78e5a0312 Mon Sep 17 00:00:00 2001 From: Gregory Lee Date: Thu, 18 Oct 2018 11:42:28 -0400 Subject: [PATCH 02/16] add XMLRECImage support to imageclasses.py --- nibabel/imageclasses.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/nibabel/imageclasses.py b/nibabel/imageclasses.py index c1a0b7133a..71313adea0 100644 --- a/nibabel/imageclasses.py +++ b/nibabel/imageclasses.py @@ -18,6 +18,7 @@ from .nifti1 import Nifti1Pair, Nifti1Image from .nifti2 import Nifti2Pair, Nifti2Image from .parrec import PARRECImage +from .xmlrec import XMLRECImage from .spm99analyze import Spm99AnalyzeImage from .spm2analyze import Spm2AnalyzeImage from .volumeutils import Recoder @@ -32,7 +33,7 @@ Cifti2Image, Nifti2Image, # Cifti2 before Nifti2 Spm2AnalyzeImage, Spm99AnalyzeImage, AnalyzeImage, Minc1Image, Minc2Image, MGHImage, - PARRECImage, GiftiImage, AFNIImage] + PARRECImage, XMLRECImage, GiftiImage, AFNIImage] # DEPRECATED: mapping of names to classes and class functionality @@ -90,6 +91,11 @@ def __getitem__(self, *args, **kwargs): 'has_affine': True, 'makeable': False, 'rw': False}, + xml={'class': XMLRECImage, + 'ext': '.xml', + 'has_affine': True, + 'makeable': False, + 'rw': False}, afni={'class': AFNIImage, 'ext': '.brik', 'has_affine': True, @@ -113,6 +119,7 @@ def __getitem__(self, *args, **kwargs): ('mgh', '.mgh'), ('mgz', '.mgz'), ('par', '.par'), + ('xml', '.xml'), ('brik', '.brik') )) @@ -121,7 +128,7 @@ def __getitem__(self, *args, **kwargs): # here. KNOWN_SPATIAL_FIRST = (Nifti1Pair, Nifti1Image, Nifti2Pair, Nifti2Image, Spm2AnalyzeImage, Spm99AnalyzeImage, AnalyzeImage, - MGHImage, PARRECImage, AFNIImage) + MGHImage, PARRECImage, XMLRECImage, AFNIImage) def spatial_axes_first(img): From df2751ef6f9d80bcb20dfd2fabadfdb5a47f927e Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 18 Oct 2018 22:59:31 -0400 Subject: [PATCH 03/16] Assign an enumerated value of -1 in the PAR image_defs header for unrecognized enum strings. add an error message with a hint regarding truncated XML files --- nibabel/xmlrec.py | 47 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 39 insertions(+), 8 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 9ed9aedf67..6dd4b24ec4 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -469,6 +469,7 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): image_def_dtd = _composite_attributes_xml_to_par(image_def_dtd) image_defs = np.zeros(len(image_defs_array), dtype=image_def_dtd) + already_warned = [] for i, image_def in enumerate(image_defs_array): if image_def.find('Key') != image_def[0]: @@ -477,11 +478,23 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): key_def = image_def[0] for key in key_def: name = key.get('Name') - if name in image_def_XML_to_PAR: + if dtype_format == 'par' and name in image_def_XML_to_PAR: name = image_def_XML_to_PAR[name] if dtype_format == 'par' and name in enums_dict_PAR: # string -> int - val = enums_dict_PAR[name][key.text] + if key.text in enums_dict_PAR[name]: + val = enums_dict_PAR[name][key.text] + else: + if (name, key.text) not in already_warned: + warnings.warn( + ("Unknown enumerated value for {} with name {}. " + "Setting the value to -1. Please contact the " + "nibabel developers about adding support for " + "this to the XML/REC reader.").format( + name, key.text)) + val = -1 + # avoid repeated warnings for this enum + already_warned.append((name, key.text)) else: val = key.text image_defs[name][i] = dtype_dict[name](val) @@ -501,7 +514,19 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): else: if dtype_format == 'par' and name in enums_dict_PAR: # string -> int - val = enums_dict_PAR[name][text] + if text in enums_dict_PAR[name]: + val = enums_dict_PAR[name][text] + else: + if (name, text) not in already_warned: + warnings.warn( + ("Unknown enumerated value for {} with " + "name {}. Setting the value to -1. " + "Please contact the nibabel developers " + "about adding support for this to the " + "XML/REC reader.").format(name, text)) + val = -1 + # avoid repeated warnings for this enum + already_warned.append((name, text)) else: val = entry_dtype(text) if dtype_format == 'par' and name in composite_PAR_keys: @@ -518,8 +543,8 @@ def parse_XML_header(fobj, dtype_format='xml'): Parameters ---------- - fobj : file-object - The PAR header file object. + fobj : file-object or str + The XML header file object or file name. dtype_format : {'xml', 'par'} If 'par' the image_defs will be converted to a format matching that found in PARRECHeader. @@ -546,9 +571,15 @@ def parse_XML_header(fobj, dtype_format='xml'): Please email the NiBabel mailing list, if you are interested in adding support for this version. """.format(version))) - general_info = _process_gen_dict_XML(root) - image_defs = _process_image_lines_xml( - root, dtype_format=dtype_format) + try: + general_info = _process_gen_dict_XML(root) + image_defs = _process_image_lines_xml( + root, dtype_format=dtype_format) + except ET.ParseError: + raise XMLRECError( + "A ParseError occured in the ElementTree library while " + "reading the XML file. This may be due to a truncated XML " + "file.") return general_info, image_defs From 8c405c118b872ace941760e734c734e87815ae41 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Thu, 18 Oct 2018 23:14:44 -0400 Subject: [PATCH 04/16] expand enums_dict using empirically determined values --- nibabel/xmlrec.py | 31 ++++++++++++++++++++----------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 6dd4b24ec4..9042cff534 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -14,19 +14,28 @@ class XMLRECError(Exception): XML/REC. """ -# Dictionary of conversions from enumerated types to integer values. -# The following named Enums are used in XML files. -# I do not know thefull definitions for these. -# This was partially adapted from: -# https://github.com/xiangruili/dicm2nii/blob/6855498a702b06bbc7d1de716f4e9bdfe50c74fc/xml2par.m#L356 +# Dictionary of conversions from enumerated types to integer value for use in +# converting from XML enum names to PAR-style integers. +# The keys in enums_dict are the names of the enum keys used in XML files. +# The present conversions for strings to enumerated values was determined +# empirically via comparison of simultaneously exported .PAR and .xml files +# from a range of different scan types. Any enum labels that are not recognized +# will result in a warning encouraging the user to report the unkown case to +# the nibabel developers. enums_dict = { - 'Label Type': {'-': 1}, - 'Type': {'M': 0, 'R': 1, 'I': 2, 'P': 3}, # TODO: 'B0', etc... - 'Sequence': {'SE': 1, 'FFE': 2}, - 'Image Type Ed Es': {'U': 2}, # TODO: others? - 'Display Orientation': {'-': 0, 'NONE': 0}, # TODO: correct value for None + 'Label Type': {'CONTROL': 1, 'LABEL': 2, '-': 1}, + 'Type': {'M': 0, 'R': 1, 'I': 2, 'P': 3, 'T1': 6, 'T2': 7, 'ADC': 11, + 'EADC': 17, 'B0': 18, 'PERFUSION': 30, 'F': 31, 'IP': 32, + 'FF': 34, 'R2': -1, 'R2_STAR': -1, 'T2_STAR': -1, 'W': -1, + 'STIFF': -1, 'WAVE': -1, 'SW_M': -1, 'SW_P': -1}, + 'Sequence': {'IR': 0, 'SE': 1, 'FFE': 2, 'PCA': 4, 'UNSPECIFIED': 5, + 'DERIVED': 7, 'B1': 9, 'MRE': 10}, + 'Image Type Ed Es': {'U': 2}, + 'Display Orientation': {'-': 0, 'NONE': 0}, 'Slice Orientation': {'Transversal': 1, 'Sagittal': 2, 'Coronal': 3}, - 'Contrast Type': {'DIFFUSION': 0, 'T2': 4, 'TAGGING': 6, 'T1': 7}, + 'Contrast Type': {'DIFFUSION': 0, 'FLOW_ENCODED': 1, 'PERFUSION': 3, + 'PROTON_DENSITY': 4, 'TAGGING': 6, 'T1': 7, 'T2': 8, + 'UNKNOWN': 11}, 'Diffusion Anisotropy Type': {'-': 0}} # Dict for converting XML/REC types to appropriate python types From a4636b272c508e7243b4b30e9dcc3a1ed8acd285 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 00:03:10 -0400 Subject: [PATCH 05/16] Fix reading strings into image_defs when dtype_format='xml' --- nibabel/xmlrec.py | 39 +++++++++++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 6 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 9042cff534..ff674f0ac1 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -52,6 +52,12 @@ class XMLRECError(Exception): 'Date': str, 'Time': str} +# choose appropriate string length for use within the structured array dtype +str_length_dict = {'Enumeration': '|S32', + 'String': '|S128', + 'Date': '|S16', + 'Time': '|S16'} + # for use with parrec.py, the Enumeration strings must be converted to ints par_type_dict = copy(xml_type_dict) par_type_dict['Enumeration'] = int @@ -276,6 +282,8 @@ def _process_gen_dict_XML(xml_root): a = e.attrib if 'Name' in a: entry_type = xml_type_dict[a['Type']] + if entry_type in ['S16', 'S32', 'S128']: + entry_type = str if 'ArraySize' in a: val = [entry_type(i) for i in e.text.strip().split()] else: @@ -320,15 +328,23 @@ def _get_image_def_attributes(xml_root, dtype_format='xml'): if not np.all(['Name' in k.attrib for k in img_keys]): raise XMLRECError("Expected each Key attribute to have a Name") + def _get_type(name, type_dict, str_length_dict=str_length_dict): + t = type_dict[name] + if t is str: # convert from str to a specific length such as |S32 + t = str_length_dict[name] + return t + dtype_format = dtype_format.lower() if dtype_format == 'xml': # xml_type_dict keeps enums as their string representation key_attributes = [ - (k.get('Name'), xml_type_dict[k.get('Type')]) for k in img_keys] + (k.get('Name'), _get_type(k.get('Type'), xml_type_dict)) + for k in img_keys] elif dtype_format == 'par': # par_type_dict converts enums to int as used in PAR/REC key_attributes = [ - (image_def_XML_to_PAR[k.get('Name')], par_type_dict[k.get('Type')]) + (image_def_XML_to_PAR[k.get('Name')], + _get_type(k.get('Type'), par_type_dict)) for k in img_keys] else: raise XMLRECError("dtype_format must be 'par' or 'xml'") @@ -343,9 +359,9 @@ def _get_image_def_attributes(xml_root, dtype_format='xml'): # print("enum_type = {}".format(enum_type)) name = a['Name'] if dtype_format == 'xml': - entry_type = xml_type_dict[a['Type']] + entry_type = _get_type(a['Type'], xml_type_dict) else: - entry_type = par_type_dict[a['Type']] + entry_type = _get_type(a['Type'], par_type_dict) if name in image_def_XML_to_PAR: name = image_def_XML_to_PAR[name] if 'ArraySize' in a: @@ -477,6 +493,17 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): # image_defs based on conversion of types in composite_PAR_keys image_def_dtd = _composite_attributes_xml_to_par(image_def_dtd) + def _get_val(entry_dtype, text): + if entry_dtype == '|S16': + val = text[:16] + elif entry_dtype == '|S32': + val = text[:32] + elif entry_dtype == '|S128': + val = text[:128] + else: + val = entry_dtype(text) + return val + image_defs = np.zeros(len(image_defs_array), dtype=image_def_dtd) already_warned = [] for i, image_def in enumerate(image_defs_array): @@ -506,7 +533,7 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): already_warned.append((name, key.text)) else: val = key.text - image_defs[name][i] = dtype_dict[name](val) + image_defs[name][i] = _get_val(dtype_dict[name], val) # for all image properties we know about for element in image_def[1:]: @@ -537,7 +564,7 @@ def _process_image_lines_xml(xml_root, dtype_format='xml'): # avoid repeated warnings for this enum already_warned.append((name, text)) else: - val = entry_dtype(text) + val = _get_val(entry_dtype, text) if dtype_format == 'par' and name in composite_PAR_keys: # conversion of types in composite_PAR_keys name, vec_index = _get_composite_key_index(name) From 570a7edafaac47ba819f8d94d1e358de1a8af515 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 08:57:33 -0400 Subject: [PATCH 06/16] XMLRECImage: file extensions should be specified in lower case --- nibabel/xmlrec.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index ff674f0ac1..bb3a3705b2 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -621,6 +621,7 @@ def parse_XML_header(fobj, dtype_format='xml'): class XMLRECHeader(PARRECHeader): + @classmethod def from_fileobj(klass, fileobj, permit_truncated=False, strict_sort=False): @@ -628,6 +629,7 @@ def from_fileobj(klass, fileobj, permit_truncated=False, # convert to PAR/REC format general_info info = general_info_xml_to_par(info) return klass(info, image_defs, permit_truncated, strict_sort) + @classmethod def from_header(klass, header=None): if header is None: @@ -636,6 +638,7 @@ def from_header(klass, header=None): return header.copy() raise XMLRECError('Cannot create XMLREC header from ' 'non-XMLREC header.') + def copy(self): return XMLRECHeader(deepcopy(self.general_info), self.image_defs.copy(), @@ -645,7 +648,7 @@ def copy(self): class XMLRECImage(PARRECImage): """XML/REC image""" header_class = XMLRECHeader - valid_exts = ('.REC', '.xml') - files_types = (('image', '.REC'), ('header', '.xml')) + valid_exts = ('.rec', '.xml') + files_types = (('image', '.rec'), ('header', '.xml')) load = XMLRECImage.load From bbeb710cf11f8606b440f9c10bba47d735888d9d Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 08:58:31 -0400 Subject: [PATCH 07/16] loading from .REC: special case to find .xml when .PAR not present --- nibabel/loadsave.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/nibabel/loadsave.py b/nibabel/loadsave.py index 8c3041e73c..036433c635 100644 --- a/nibabel/loadsave.py +++ b/nibabel/loadsave.py @@ -19,6 +19,7 @@ from .arrayproxy import is_proxy from .py3k import FileNotFoundError from .deprecated import deprecate_with_version +from .parrec import PARRECImage def load(filename, **kwargs): @@ -46,6 +47,15 @@ def load(filename, **kwargs): for image_klass in all_image_classes: is_valid, sniff = image_klass.path_maybe_image(filename, sniff) if is_valid: + if image_klass is PARRECImage and '.REC' in filename: + # a .REC file can have either a .PAR of .xml header. + # This skip case assumes PARRECImage is beforeXMLRECImage in + # all_image_classes. + par_exists = os.path.exists(filename.replace('.REC', '.PAR')) + xml_exists = os.path.exists(filename.replace('.REC', '.xml')) + if not par_exists and xml_exists: + continue # skip trying .PAR and proceed to .xml + print(image_klass) img = image_klass.from_filename(filename, **kwargs) return img From a7a9edce7e27b08e942ddb07ee653ee0be082fd2 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 09:00:43 -0400 Subject: [PATCH 08/16] remove explicit mention of PAR from warning/errors that could also occur for XML/REC --- nibabel/parrec.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabel/parrec.py b/nibabel/parrec.py index 87e1ac81e6..72047bda5a 100644 --- a/nibabel/parrec.py +++ b/nibabel/parrec.py @@ -865,7 +865,7 @@ def get_data_offset(self): def set_data_offset(self, offset): """ PAR header always has 0 data offset (into REC file) """ if offset != 0: - raise PARRECError("PAR header assumes offset 0") + raise PARRECError("header assumes offset 0") def _calc_zooms(self): """Compute image zooms from header data. @@ -894,7 +894,7 @@ def _calc_zooms(self): # If 4D dynamic scan, convert time from milliseconds to seconds if len(zooms) > 3 and self.general_info['dyn_scan']: if len(self.general_info['repetition_time']) > 1: - warnings.warn("multiple TRs found in .PAR file") + warnings.warn("multiple TRs found in header file") zooms[3] = self.general_info['repetition_time'][0] / 1000. return zooms From ed43e97467efe65ddda994a836e1fda46d4e7812 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 10:14:30 -0400 Subject: [PATCH 09/16] Add XML/REC support to parrec2nii.py Add an xmlrec2nii command line script. This currently is identical to parrec2nii because the underlying function now supports both formats. --- bin/xmlrec2nii | 10 ++++++++++ nibabel/cmdline/parrec2nii.py | 22 +++++++++++++++++----- 2 files changed, 27 insertions(+), 5 deletions(-) create mode 100755 bin/xmlrec2nii diff --git a/bin/xmlrec2nii b/bin/xmlrec2nii new file mode 100755 index 0000000000..ae04adf024 --- /dev/null +++ b/bin/xmlrec2nii @@ -0,0 +1,10 @@ +#!python +"""XML/REC to NIfTI converter +""" + +# parrec2nii reads XML/REC as well +from nibabel.cmdline.parrec2nii import main + + +if __name__ == '__main__': + main() diff --git a/nibabel/cmdline/parrec2nii.py b/nibabel/cmdline/parrec2nii.py index dc26870ff6..7ce947cc99 100644 --- a/nibabel/cmdline/parrec2nii.py +++ b/nibabel/cmdline/parrec2nii.py @@ -10,6 +10,7 @@ import csv import nibabel import nibabel.parrec as pr +import nibabel.xmlrec as xr from nibabel.parrec import one_line from nibabel.mriutils import calculate_dwell_time, MRIError import nibabel.nifti1 as nifti1 @@ -147,7 +148,18 @@ def error(msg, exit_code): def proc_file(infile, opts): # figure out the output filename, and see if it exists - basefilename = splitext_addext(os.path.basename(infile))[0] + basefilename, ext = splitext_addext(os.path.basename(infile))[:2] + + ext = ext.lower() + if (ext == '.xml' or + (ext == '.rec' and os.path.exists(basefilename + '.xml'))): + pr_module = xr + if opts.permit_truncated: + raise ValueError("The permit_truncated option is not currently " + "supported for .xml/.REC data") + else: + pr_module = pr + if opts.outdir is not None: # set output path basefilename = os.path.join(opts.outdir, basefilename) @@ -165,10 +177,10 @@ def proc_file(infile, opts): # load the PAR header and data scaling = 'dv' if opts.scaling == 'off' else opts.scaling infile = fname_ext_ul_case(infile) - pr_img = pr.load(infile, - permit_truncated=opts.permit_truncated, - scaling=scaling, - strict_sort=opts.strict_sort) + pr_img = pr_module.load(infile, + permit_truncated=opts.permit_truncated, + scaling=scaling, + strict_sort=opts.strict_sort) pr_hdr = pr_img.header affine = pr_hdr.get_affine(origin=opts.origin) slope, intercept = pr_hdr.get_data_scaling(scaling) From 5fa6810529004dd89ad6c032fbd7fcd080466c91 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 13:10:04 -0400 Subject: [PATCH 10/16] refactor/simplify. Do not convert to PAR/REC style names --- nibabel/xmlrec.py | 947 +++++++++++++++++++++++++--------------------- 1 file changed, 526 insertions(+), 421 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index bb3a3705b2..b3e649ba50 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -1,10 +1,16 @@ from copy import copy, deepcopy import warnings import xml.etree.ElementTree as ET +from collections import OrderedDict import numpy as np -from nibabel.parrec import one_line -from nibabel.parrec import PARRECImage, PARRECHeader + +from .affines import from_matvec, dot_reduce, apply_affine +from .deprecated import deprecate_with_version +from .nifti1 import unit_codes +from .eulerangles import euler2mat +from .parrec import (one_line, PARRECImage, PARRECHeader, ACQ_TO_PSL, + PSL_TO_RAS, vol_is_full, vol_numbers, DEG2RAD) class XMLRECError(Exception): @@ -14,29 +20,10 @@ class XMLRECError(Exception): XML/REC. """ -# Dictionary of conversions from enumerated types to integer value for use in -# converting from XML enum names to PAR-style integers. -# The keys in enums_dict are the names of the enum keys used in XML files. -# The present conversions for strings to enumerated values was determined -# empirically via comparison of simultaneously exported .PAR and .xml files -# from a range of different scan types. Any enum labels that are not recognized -# will result in a warning encouraging the user to report the unkown case to -# the nibabel developers. -enums_dict = { - 'Label Type': {'CONTROL': 1, 'LABEL': 2, '-': 1}, - 'Type': {'M': 0, 'R': 1, 'I': 2, 'P': 3, 'T1': 6, 'T2': 7, 'ADC': 11, - 'EADC': 17, 'B0': 18, 'PERFUSION': 30, 'F': 31, 'IP': 32, - 'FF': 34, 'R2': -1, 'R2_STAR': -1, 'T2_STAR': -1, 'W': -1, - 'STIFF': -1, 'WAVE': -1, 'SW_M': -1, 'SW_P': -1}, - 'Sequence': {'IR': 0, 'SE': 1, 'FFE': 2, 'PCA': 4, 'UNSPECIFIED': 5, - 'DERIVED': 7, 'B1': 9, 'MRE': 10}, - 'Image Type Ed Es': {'U': 2}, - 'Display Orientation': {'-': 0, 'NONE': 0}, - 'Slice Orientation': {'Transversal': 1, 'Sagittal': 2, 'Coronal': 3}, - 'Contrast Type': {'DIFFUSION': 0, 'FLOW_ENCODED': 1, 'PERFUSION': 3, - 'PROTON_DENSITY': 4, 'TAGGING': 6, 'T1': 7, 'T2': 8, - 'UNKNOWN': 11}, - 'Diffusion Anisotropy Type': {'-': 0}} +slice_orientation_rename = dict(Transversal='transverse', + Sagittal='sagittal', + Coronal='coronal') + # Dict for converting XML/REC types to appropriate python types # could convert the Enumeration strings to int, but enums_dict is incomplete @@ -58,216 +45,8 @@ class XMLRECError(Exception): 'Date': '|S16', 'Time': '|S16'} -# for use with parrec.py, the Enumeration strings must be converted to ints -par_type_dict = copy(xml_type_dict) -par_type_dict['Enumeration'] = int - supported_xml_versions = ['PRIDE_V5'] -# convert XML names to expected property names as used in nibabel.parrec -# Note: additional special cases are handled in general_info_xml_to_par() -general_info_XML_to_nibabel = { - 'Patient Name': 'patient_name', - 'Examination Name': 'exam_name', - 'Protocol Name': 'protocol_name', - 'Aquisition Number': 'acq_nr', - 'Reconstruction Number': 'recon_nr', - 'Scan Duration': 'scan_duration', - 'Max No Phases': 'max_cardiac_phases', - 'Max No Echoes': 'max_echoes', - 'Max No Slices': 'max_slices', - 'Max No Dynamics': 'max_dynamics', - 'Max No Mixes': 'max_mixes', - 'Patient Position': 'patient_position', - 'Preparation Direction': 'prep_direction', - 'Technique': 'tech', - 'Scan Mode': 'scan_mode', - 'Repetition Times': 'repetition_time', - 'Water Fat Shift': 'water_fat_shift', - 'Flow Compensation': 'flow_compensation', - 'Presaturation': 'presaturation', - 'Phase Encoding Velocity': 'phase_enc_velocity', - 'MTC': 'mtc', - 'SPIR': 'spir', - 'EPI factor': 'epi_factor', - 'Dynamic Scan': 'dyn_scan', - 'Diffusion': 'diffusion', - 'Diffusion Echo Time': 'diffusion_echo_time', - 'Max No B Values': 'max_diffusion_values', - 'Max No Gradient Orients': 'max_gradient_orient', - 'No Label Types': 'nr_label_types', - 'Series Data Type': 'series_type'} - - -def general_info_xml_to_par(xml_info): - """Convert general_info from XML-style names to PAR-style names.""" - xml_info_init = xml_info - xml_info = deepcopy(xml_info) - general_info = {} - for k in xml_info_init.keys(): - # convert all keys with a simple 1-1 name conversion - if k in general_info_XML_to_nibabel: - general_info[general_info_XML_to_nibabel[k]] = xml_info.pop(k) - try: - general_info['exam_date'] = '{} / {}'.format( - xml_info.pop('Examination Date'), - xml_info.pop('Examination Time')) - except KeyError: - pass - - try: - general_info['angulation'] = np.asarray( - [xml_info.pop('Angulation AP'), - xml_info.pop('Angulation FH'), - xml_info.pop('Angulation RL')]) - except KeyError: - pass - - try: - general_info['off_center'] = np.asarray( - [xml_info.pop('Off Center AP'), - xml_info.pop('Off Center FH'), - xml_info.pop('Off Center RL')]) - except KeyError: - pass - - try: - general_info['fov'] = np.asarray( - [xml_info.pop('FOV AP'), - xml_info.pop('FOV FH'), - xml_info.pop('FOV RL')]) - except KeyError: - pass - - try: - general_info['scan_resolution'] = np.asarray( - [xml_info.pop('Scan Resolution X'), - xml_info.pop('Scan Resolution Y')], - dtype=int) - except KeyError: - pass - - # copy any excess keys not normally in the .PARREC general info - # These will not be used by the PARREC code, but are kept for completeness - general_info.update(xml_info) - - return general_info - - -# TODO: remove this function? It is currently unused, but may be useful for -# testing roundtrip convesion. -def general_info_par_to_xml(par_info): - """Convert general_info from PAR-style names to XML-style names.""" - general_info_nibabel_to_XML = { - v: k for k, v in general_info_XML_to_nibabel.items()} - par_info_init = par_info - par_info = deepcopy(par_info) - general_info = {} - for k in par_info_init.keys(): - # convert all keys with a simple 1-1 name conversion - if k in general_info_nibabel_to_XML: - general_info[general_info_nibabel_to_XML[k]] = par_info.pop(k) - try: - tmp = par_info['exam_date'].split('/') - general_info['Examination Date'] = tmp[0].strip() - general_info['Examination Time'] = tmp[1].strip() - except KeyError: - pass - - try: - general_info['Angulation AP'] = par_info['angulation'][0] - general_info['Angulation FH'] = par_info['angulation'][1] - general_info['Angulation RL'] = par_info['angulation'][2] - par_info.pop('angulation') - except KeyError: - pass - - try: - general_info['Off Center AP'] = par_info['off_center'][0] - general_info['Off Center FH'] = par_info['off_center'][1] - general_info['Off Center RL'] = par_info['off_center'][2] - par_info.pop('off_center') - except KeyError: - pass - - try: - general_info['FOV AP'] = par_info['fov'][0] - general_info['FOV FH'] = par_info['fov'][1] - general_info['FOV RL'] = par_info['fov'][2] - par_info.pop('fov') - except KeyError: - pass - - try: - general_info['Scan Resolution X'] = par_info['scan_resolution'][0] - general_info['Scan Resolution Y'] = par_info['scan_resolution'][1] - par_info.pop('scan_resolution') - except KeyError: - pass - - # copy any unrecognized keys as is. - # known keys found in XML PRIDE_V5, but not in PAR v4.2 are: - # 'Samples Per Pixel' and 'Image Planar Configuration' - general_info.update(par_info) - - return general_info - -# dictionary mapping fieldnames in the XML header to the corresponding names -# in a PAR V4.2 file -image_def_XML_to_PAR = { - 'Image Type Ed Es': 'image_type_ed_es', - 'No Averages': 'number of averages', - 'Max RR Interval': 'maximum RR-interval', - 'Type': 'image_type_mr', - 'Display Orientation': 'image_display_orientation', - 'Image Flip Angle': 'image_flip_angle', - 'Rescale Slope': 'rescale slope', - 'Label Type': 'label type', - 'fMRI Status Indication': 'fmri_status_indication', - 'TURBO Factor': 'TURBO factor', - 'Min RR Interval': 'minimum RR-interval', - 'Scale Slope': 'scale slope', - 'Inversion Delay': 'Inversion delay', - 'Window Width': 'window width', - 'Sequence': 'scanning sequence', - 'Diffusion Anisotropy Type': 'diffusion anisotropy type', - 'Index': 'index in REC file', - 'Rescale Intercept': 'rescale intercept', - 'Diffusion B Factor': 'diffusion_b_factor', - 'Trigger Time': 'trigger_time', - 'Echo': 'echo number', - 'Echo Time': 'echo_time', - 'Pixel Spacing': 'pixel spacing', - 'Slice Gap': 'slice gap', - 'Dyn Scan Begin Time': 'dyn_scan_begin_time', - 'Window Center': 'window center', - 'Contrast Type': 'contrast type', - 'Slice': 'slice number', - 'BValue': 'diffusion b value number', - 'Scan Percentage': 'scan percentage', - 'Phase': 'cardiac phase number', - 'Slice Thickness': 'slice thickness', - 'Slice Orientation': 'slice orientation', - 'Dynamic': 'dynamic scan number', - 'Pixel Size': 'image pixel size', - 'Grad Orient': 'gradient orientation number', - 'Cardiac Frequency': 'cardiac frequency'} - -# copy of enums_dict but with key names converted to their PAR equivalents -enums_dict_PAR = {image_def_XML_to_PAR[k]: v for k, v in enums_dict.items()} - - -# TODO?: The following values have different names in the XML vs. PAR header -rename_XML_to_PAR = { - 'HFS': 'Head First Supine', - 'LR': 'Left-Right', - 'RL': 'Right-Left', - 'AP': 'Anterior-Posterior', - 'PA': 'Posterior-Anterior', - 'N': 0, - 'Y': 1, -} - def _process_gen_dict_XML(xml_root): """Read the general_info from an XML file. @@ -292,7 +71,7 @@ def _process_gen_dict_XML(xml_root): return general_info -def _get_image_def_attributes(xml_root, dtype_format='xml'): +def _get_image_def_attributes(xml_root): """Get names and dtypes for all attributes defined for each image. called by _process_image_lines_xml @@ -300,10 +79,7 @@ def _get_image_def_attributes(xml_root, dtype_format='xml'): Paramters --------- xml_root : - TODO - dtype_format : {'xml', 'par'} - If 'par'', XML paramter names and dtypes are converted to their PAR - equivalents. + The root of the XML tree Returns ------- @@ -334,20 +110,10 @@ def _get_type(name, type_dict, str_length_dict=str_length_dict): t = str_length_dict[name] return t - dtype_format = dtype_format.lower() - if dtype_format == 'xml': - # xml_type_dict keeps enums as their string representation - key_attributes = [ - (k.get('Name'), _get_type(k.get('Type'), xml_type_dict)) - for k in img_keys] - elif dtype_format == 'par': - # par_type_dict converts enums to int as used in PAR/REC - key_attributes = [ - (image_def_XML_to_PAR[k.get('Name')], - _get_type(k.get('Type'), par_type_dict)) - for k in img_keys] - else: - raise XMLRECError("dtype_format must be 'par' or 'xml'") + # xml_type_dict keeps enums as their string representation + key_attributes = [ + (k.get('Name'), _get_type(k.get('Type'), xml_type_dict)) + for k in img_keys] # Process other attributes that are not considered image keys other_attributes = [] @@ -358,12 +124,7 @@ def _get_type(name, type_dict, str_length_dict=str_length_dict): # enum_type = a['EnumType'] # print("enum_type = {}".format(enum_type)) name = a['Name'] - if dtype_format == 'xml': - entry_type = _get_type(a['Type'], xml_type_dict) - else: - entry_type = _get_type(a['Type'], par_type_dict) - if name in image_def_XML_to_PAR: - name = image_def_XML_to_PAR[name] + entry_type = _get_type(a['Type'], xml_type_dict) if 'ArraySize' in a: # handle vector entries (e.g. 'Pixel Size' is length 2) entry = (name, entry_type, int(a['ArraySize'])) @@ -374,125 +135,24 @@ def _get_type(name, type_dict, str_length_dict=str_length_dict): return key_attributes, other_attributes -# these keys are elements of a multi-valued key in the PAR/REC image_defs -composite_PAR_keys = ['Diffusion AP', 'Diffusion FH', 'Diffusion RL', - 'Angulation AP', 'Angulation FH', 'Angulation RL', - 'Offcenter AP', 'Offcenter FH', 'Offcenter RL', - 'Resolution X', 'Resolution Y'] - - -def _composite_attributes_xml_to_par(other_attributes): - """utility used in conversion from XML-style to PAR-style image_defs. - - called by _process_image_lines_xml - """ - if ('Diffusion AP', np.float32) in other_attributes: - other_attributes.remove(('Diffusion AP', np.float32)) - other_attributes.remove(('Diffusion FH', np.float32)) - other_attributes.remove(('Diffusion RL', np.float32)) - other_attributes.append(('diffusion', float, (3, ))) - - if ('Angulation AP', np.float64) in other_attributes: - other_attributes.remove(('Angulation AP', np.float64)) - other_attributes.remove(('Angulation FH', np.float64)) - other_attributes.remove(('Angulation RL', np.float64)) - other_attributes.append(('image angulation', float, (3, ))) - - if ('Offcenter AP', np.float64) in other_attributes: - other_attributes.remove(('Offcenter AP', np.float64)) - other_attributes.remove(('Offcenter FH', np.float64)) - other_attributes.remove(('Offcenter RL', np.float64)) - other_attributes.append(('image offcentre', float, (3, ))) - - if ('Resolution X', np.uint16) in other_attributes: - other_attributes.remove(('Resolution X', np.uint16)) - other_attributes.remove(('Resolution Y', np.uint16)) - other_attributes.append(('recon resolution', int, (2, ))) - - return other_attributes - - -# TODO: remove? this one is currently unused, but perhaps useful for testing -def _composite_attributes_par_to_xml(other_attributes): - if ('diffusion', float, (3, )) in other_attributes: - other_attributes.append(('Diffusion AP', np.float32)) - other_attributes.append(('Diffusion FH', np.float32)) - other_attributes.append(('Diffusion RL', np.float32)) - other_attributes.remove(('diffusion', float, (3, ))) - - if ('image angulation', float, (3, )) in other_attributes: - other_attributes.append(('Angulation AP', np.float64)) - other_attributes.append(('Angulation FH', np.float64)) - other_attributes.append(('Angulation RL', np.float64)) - other_attributes.remove(('image angulation', float, (3, ))) - - if ('image offcentre', float, (3, )) in other_attributes: - other_attributes.append(('Offcenter AP', np.float64)) - other_attributes.append(('Offcenter FH', np.float64)) - other_attributes.append(('Offcenter RL', np.float64)) - other_attributes.remove(('image offcentre', float, (3, ))) - - if ('recon resolution', int, (2, )) in other_attributes: - other_attributes.append(('Resolution X', np.uint16)) - other_attributes.append(('Resolution Y', np.uint16)) - other_attributes.remove(('recon resolution', int, (2, ))) - return other_attributes - - -def _get_composite_key_index(key): - """utility used in conversion from XML-style to PAR-style image_defs. - - called by _process_image_lines_xml - """ - if 'Diffusion' in key: - name = 'diffusion' - elif 'Angulation' in key: - name = 'image angulation' - elif 'Offcenter' in key: - name = 'image offcentre' - elif 'Resolution' in key: - name = 'recon resolution' - else: - raise ValueError("unrecognized composite element name: {}".format(key)) - - if key in ['Diffusion AP', 'Angulation AP', 'Offcenter AP', - 'Resolution X']: - index = 0 - elif key in ['Diffusion FH', 'Angulation FH', 'Offcenter FH', - 'Resolution Y']: - index = 1 - elif key in ['Diffusion RL', 'Angulation RL', 'Offcenter RL']: - index = 2 - else: - raise ValueError("unrecognized composite element name: {}".format(key)) - return (name, index) - - -def _process_image_lines_xml(xml_root, dtype_format='xml'): +def _process_image_lines_xml(xml_root): """Build image_defs by parsing the XML file. Parameters ---------- xml_root : - TODO - dtype_format : {'xml', 'par'} - If 'par'', XML paramter names and dtypes are converted to their PAR file - equivalents. + The root of the XML tree """ image_defs_array = xml_root.find('Image_Array') if image_defs_array is None: raise RuntimeError("No 'Image_Array' found in the XML file") - key_attributes, other_attributes = _get_image_def_attributes(xml_root, dtype_format=dtype_format) + key_attributes, other_attributes = _get_image_def_attributes(xml_root) image_def_dtd = key_attributes + other_attributes # dtype dict based on the XML attribute names dtype_dict = {a[0]: a[1] for a in image_def_dtd} - if dtype_format == 'par': - # image_defs based on conversion of types in composite_PAR_keys - image_def_dtd = _composite_attributes_xml_to_par(image_def_dtd) - def _get_val(entry_dtype, text): if entry_dtype == '|S16': val = text[:16] @@ -514,25 +174,7 @@ def _get_val(entry_dtype, text): key_def = image_def[0] for key in key_def: name = key.get('Name') - if dtype_format == 'par' and name in image_def_XML_to_PAR: - name = image_def_XML_to_PAR[name] - if dtype_format == 'par' and name in enums_dict_PAR: - # string -> int - if key.text in enums_dict_PAR[name]: - val = enums_dict_PAR[name][key.text] - else: - if (name, key.text) not in already_warned: - warnings.warn( - ("Unknown enumerated value for {} with name {}. " - "Setting the value to -1. Please contact the " - "nibabel developers about adding support for " - "this to the XML/REC reader.").format( - name, key.text)) - val = -1 - # avoid repeated warnings for this enum - already_warned.append((name, key.text)) - else: - val = key.text + val = key.text image_defs[name][i] = _get_val(dtype_dict[name], val) # for all image properties we know about @@ -541,49 +183,22 @@ def _get_val(entry_dtype, text): text = element.text if 'Name' in a: name = a['Name'] - if dtype_format == 'par' and name in image_def_XML_to_PAR: - name = image_def_XML_to_PAR[name] - # composite_PAR_keys entry_dtype = dtype_dict[name] if 'ArraySize' in a: val = [entry_dtype(i) for i in text.strip().split()] else: - if dtype_format == 'par' and name in enums_dict_PAR: - # string -> int - if text in enums_dict_PAR[name]: - val = enums_dict_PAR[name][text] - else: - if (name, text) not in already_warned: - warnings.warn( - ("Unknown enumerated value for {} with " - "name {}. Setting the value to -1. " - "Please contact the nibabel developers " - "about adding support for this to the " - "XML/REC reader.").format(name, text)) - val = -1 - # avoid repeated warnings for this enum - already_warned.append((name, text)) - else: - val = _get_val(entry_dtype, text) - if dtype_format == 'par' and name in composite_PAR_keys: - # conversion of types in composite_PAR_keys - name, vec_index = _get_composite_key_index(name) - image_defs[name][i][vec_index] = val - else: - image_defs[name][i] = val + val = _get_val(entry_dtype, text) + image_defs[name][i] = val return image_defs -def parse_XML_header(fobj, dtype_format='xml'): +def parse_XML_header(fobj): """Parse a XML header and aggregate all information into useful containers. Parameters ---------- fobj : file-object or str The XML header file object or file name. - dtype_format : {'xml', 'par'} - If 'par' the image_defs will be converted to a format matching that - found in PARRECHeader. Returns ------- @@ -609,8 +224,7 @@ def parse_XML_header(fobj, dtype_format='xml'): """.format(version))) try: general_info = _process_gen_dict_XML(root) - image_defs = _process_image_lines_xml( - root, dtype_format=dtype_format) + image_defs = _process_image_lines_xml(root) except ET.ParseError: raise XMLRECError( "A ParseError occured in the ElementTree library while " @@ -620,15 +234,80 @@ def parse_XML_header(fobj, dtype_format='xml'): return general_info, image_defs -class XMLRECHeader(PARRECHeader): +def _truncation_checks(general_info, image_defs, permit_truncated): + """ Check for presence of truncation in PAR file parameters + + Raise error if truncation present and `permit_truncated` is False. + """ + def _err_or_warn(msg): + if not permit_truncated: + raise XMLRECError(msg) + warnings.warn(msg) + + def _chk_trunc(idef_name, gdef_max_name): + if gdef_max_name not in general_info: + return + id_values = image_defs[idef_name] + n_have = len(set(id_values)) + n_expected = general_info[gdef_max_name] + if n_have != n_expected: + _err_or_warn( + "Header inconsistency: Found {0} {1} values, " + "but expected {2}".format(n_have, idef_name, n_expected)) + + _chk_trunc('Slice', 'Max No Slices') + _chk_trunc('Echo', 'Max No Echoes') + _chk_trunc('Dynamic', 'Max No Dynamics') + _chk_trunc('BValue', 'Max No B Values') + _chk_trunc('Grad Orient', 'Max No Gradient Orients') + + # Final check for partial volumes + if not np.all(vol_is_full(image_defs['Slice'], + general_info['Max No Slices'])): + _err_or_warn("Found one or more partial volume(s)") - @classmethod - def from_fileobj(klass, fileobj, permit_truncated=False, - strict_sort=False): - info, image_defs = parse_XML_header(fileobj, dtype_format='par') - # convert to PAR/REC format general_info - info = general_info_xml_to_par(info) - return klass(info, image_defs, permit_truncated, strict_sort) + +class XMLRECHeader(PARRECHeader): + """PAR/REC header""" + + def __init__(self, info, image_defs, permit_truncated=False, + strict_sort=False): + """ + Parameters + ---------- + info : dict + "General information" from the PAR file (as returned by + `parse_PAR_header()`). + image_defs : array + Structured array with image definitions from the PAR file (as + returned by `parse_PAR_header()`). + permit_truncated : bool, optional + If True, a warning is emitted instead of an error when a truncated + recording is detected. + strict_sort : bool, optional, keyword-only + If True, a larger number of header fields are used while sorting + the REC data array. This may produce a different sort order than + `strict_sort=False`, where volumes are sorted by the order in which + the slices appear in the .PAR file. + """ + self.general_info = info.copy() + self.image_defs = image_defs.copy() + self.permit_truncated = permit_truncated + self.strict_sort = strict_sort + _truncation_checks(info, image_defs, permit_truncated) + # charge with basic properties to be able to use base class + # functionality + # dtype + bitpix = self._get_unique_image_prop('Pixel Size') + if bitpix not in (8, 16): + raise XMLRECError('Only 8- and 16-bit data supported (not %s)' + 'please report this to the nibabel developers' + % bitpix) + # REC data always little endian + dt = np.dtype('uint' + str(bitpix)).newbyteorder('<') + super(PARRECHeader, self).__init__(data_dtype=dt, + shape=self._calc_data_shape(), + zooms=self._calc_zooms()) @classmethod def from_header(klass, header=None): @@ -639,12 +318,438 @@ def from_header(klass, header=None): raise XMLRECError('Cannot create XMLREC header from ' 'non-XMLREC header.') + @classmethod + def from_fileobj(klass, fileobj, permit_truncated=False, + strict_sort=False): + info, image_defs = parse_XML_header(fileobj) + return klass(info, image_defs, permit_truncated, strict_sort) + def copy(self): return XMLRECHeader(deepcopy(self.general_info), self.image_defs.copy(), self.permit_truncated, self.strict_sort) + def as_analyze_map(self): + """Convert PAR parameters to NIFTI1 format""" + # Entries in the dict correspond to the parameters found in + # the NIfTI1 header, specifically in nifti1.py `header_dtd` defs. + # Here we set the parameters we can to simplify PAR/REC + # to NIfTI conversion. + exam_date = '{}/{}'.format( + self.general_info['Examination Date'], + self.general_info['Examination Time']) + descr = ("%s;%s;%s;%s" + % (self.general_info['Examination Name'], + self.general_info['Patient Name'], + exam_date, + self.general_info['Protocol Name']))[:80] # max len + is_fmri = (self.general_info['Max No Dynamics'] > 1) + t = 'msec' if is_fmri else 'unknown' + xyzt_units = unit_codes['mm'] + unit_codes[t] + return dict(descr=descr, xyzt_units=xyzt_units) # , pixdim=pixdim) + + def get_water_fat_shift(self): + """Water fat shift, in pixels""" + return self.general_info['Water Fat Shift'] + + def get_echo_train_length(self): + """Echo train length of the recording""" + return self.general_info['EPI factor'] + + def get_bvals_bvecs(self): + """Get bvals and bvecs from data + + Returns + ------- + b_vals : None or array + Array of b values, shape (n_directions,), or None if not a + diffusion acquisition. + b_vectors : None or array + Array of b vectors, shape (n_directions, 3), or None if not a + diffusion acquisition. + """ + if self.general_info['Diffusion'] == 0: + return None, None + reorder = self.get_sorted_slice_indices() + n_slices, n_vols = self.get_data_shape()[-2:] + bvals = self.image_defs['Diffusion B Factor'][reorder].reshape( + (n_slices, n_vols), order='F') + # All bvals within volume should be the same + assert not np.any(np.diff(bvals, axis=0)) + bvals = bvals[0] + if 'Diffusion' not in self.image_defs.dtype.names: + return bvals, None + bvecs = self.image_defs['Diffusion'][reorder].reshape( + (n_slices, n_vols, 3), order='F') + # All 3 values of bvecs should be same within volume + assert not np.any(np.diff(bvecs, axis=0)) + bvecs = bvecs[0] + # rotate bvecs to match stored image orientation + permute_to_psl = ACQ_TO_PSL[self.get_slice_orientation()] + bvecs = apply_affine(np.linalg.inv(permute_to_psl), bvecs) + return bvals, bvecs + + @deprecate_with_version('get_voxel_size deprecated. ' + 'Please use "get_zooms" instead.', + '2.0', '4.0') + + def _get_unique_resolution(self): + """Return the 2D image plane shape. + + An error is raised if the shape is not unique. + """ + resx = self.image_defs['Resolution X'] + resy = self.image_defs['Resolution Y'] + if len(set(zip(resx, resy))) > 1: + raise XMLRECError('Varying resolution in image sequence. This is ' + 'not suppported.') + return (resx[0], resy[0]) + + def get_voxel_size(self): + """Returns the spatial extent of a voxel. + + Does not include the slice gap in the slice extent. + + If you need the slice thickness not including the slice gap, use + ``self.image_defs['slice thickness']``. + + Returns + ------- + vox_size: shape (3,) ndarray + """ + # slice orientation for the whole image series + slice_thickness = self._get_unique_image_prop('Slice Thickness') + voxsize_inplane = self._get_unique_image_prop('Pixel Spacing') + voxsize = np.array((voxsize_inplane[0], + voxsize_inplane[1], + slice_thickness)) + return voxsize + + + def _calc_zooms(self): + """Compute image zooms from header data. + + Spatial axis are first three. + + Returns + ------- + zooms : array + Length 3 array for 3D image, length 4 array for 4D image. + + Notes + ----- + This routine gets called in ``__init__``, so may not be able to use + some attributes available in the fully initialized object. + """ + # slice orientation for the whole image series + slice_gap = self._get_unique_image_prop('Slice Gap') + # scaling per image axis + n_dim = 4 if self._get_n_vols() > 1 else 3 + zooms = np.ones(n_dim) + # spatial sizes are inplane X mm, inplane Y mm + inter slice gap + zooms[:2] = self._get_unique_image_prop('Pixel Spacing') + slice_thickness = self._get_unique_image_prop('Slice Thickness') + zooms[2] = slice_thickness + slice_gap + # If 4D dynamic scan, convert time from milliseconds to seconds + if len(zooms) > 3 and self.general_info['Dynamic Scan']: + if len(self.general_info['Repetition Times']) > 1: + warnings.warn("multiple TRs found in header file") + zooms[3] = self.general_info['Repetition Times'][0] / 1000. + return zooms + + def get_affine(self, origin='scanner'): + """Compute affine transformation into scanner space. + + The method only considers global rotation and offset settings in the + header and ignores potentially deviating information in the image + definitions. + + Parameters + ---------- + origin : {'scanner', 'fov'} + Transformation origin. By default the transformation is computed + relative to the scanner's iso center. If 'fov' is requested the + transformation origin will be the center of the field of view + instead. + + Returns + ------- + aff : (4, 4) array + 4x4 array, with output axis order corresponding to RAS or (x,y,z) + or (lr, pa, fh). + + Notes + ----- + Transformations appear to be specified in (ap, fh, rl) axes. The + orientation of data is recorded in the "slice orientation" field of the + PAR header "General Information". + + We need to: + + * translate to coordinates in terms of the center of the FOV + * apply voxel size scaling + * reorder / flip the data to Philips' PSL axes + * apply the rotations + * apply any isocenter scaling offset if `origin` == "scanner" + * reorder and flip to RAS axes + """ + # shape, zooms in original data ordering (ijk ordering) + ijk_shape = np.array(self.get_data_shape()[:3]) + to_center = from_matvec(np.eye(3), -(ijk_shape - 1) / 2.) + zoomer = np.diag(list(self.get_zooms()[:3]) + [1]) + slice_orientation = self.get_slice_orientation() + permute_to_psl = ACQ_TO_PSL.get(slice_orientation) + if permute_to_psl is None: + raise XMLRECError( + "Unknown slice orientation ({0}).".format(slice_orientation)) + # hdr has deg, we need radians + # Order is [ap, fh, rl] + ap_rot = self.general_info['Angulation AP'] * DEG2RAD + fh_rot = self.general_info['Angulation FH'] * DEG2RAD + rl_rot = self.general_info['Angulation RL'] * DEG2RAD + Mx = euler2mat(x=ap_rot) + My = euler2mat(y=fh_rot) + Mz = euler2mat(z=rl_rot) + # By trial and error, this unexpected order of rotations seem to give + # the closest to the observed (converted NIfTI) affine. + rot = from_matvec(dot_reduce(Mz, Mx, My)) + # compose the PSL affine + psl_aff = dot_reduce(rot, permute_to_psl, zoomer, to_center) + if origin == 'scanner': + # offset to scanner's isocenter (in ap, fh, rl) + iso_offset = np.asarray([self.general_info['Off Center AP'], + self.general_info['Off Center FH'], + self.general_info['Off Center RL']]) + psl_aff[:3, 3] += iso_offset + # Currently in PSL; apply PSL -> RAS + return np.dot(PSL_TO_RAS, psl_aff) + + def _get_n_slices(self): + """ Get number of slices for output data """ + return len(set(self.image_defs['Slice'])) + + def _get_n_vols(self): + """ Get number of volumes for output data """ + slice_nos = self.image_defs['Slice'] + vol_nos = vol_numbers(slice_nos) + is_full = vol_is_full(slice_nos, self.general_info['Max No Slices']) + return len(set(np.array(vol_nos)[is_full])) + + def _calc_data_shape(self): + """ Calculate the output shape of the image data + + Returns length 3 tuple for 3D image, length 4 tuple for 4D. + + Returns + ------- + n_inplaneX : int + number of voxels in X direction. + n_inplaneY : int + number of voxels in Y direction. + n_slices : int + number of slices. + n_vols : int + number of volumes or absent for 3D image. + + Notes + ----- + This routine gets called in ``__init__``, so may not be able to use + some attributes available in the fully initialized object. + """ + inplane_shape = self._get_unique_resolution() + shape = inplane_shape + (self._get_n_slices(),) + n_vols = self._get_n_vols() + return shape + (n_vols,) if n_vols > 1 else shape + + def get_data_scaling(self, method="dv"): + """Returns scaling slope and intercept. + + Parameters + ---------- + method : {'fp', 'dv'} + Scaling settings to be reported -- see notes below. + + Returns + ------- + slope : array + scaling slope + intercept : array + scaling intercept + + Notes + ----- + The PAR header contains two different scaling settings: 'dv' (value on + console) and 'fp' (floating point value). Here is how they are defined: + + DV = PV * RS + RI + FP = DV / (RS * SS) + + where: + + PV: value in REC + RS: rescale slope + RI: rescale intercept + SS: scale slope + """ + # These will be 3D or 4D + scale_slope = self.image_defs['Scale Slope'] + rescale_slope = self.image_defs['Rescale Slope'] + rescale_intercept = self.image_defs['Rescale Intercept'] + if method == 'dv': + slope, intercept = rescale_slope, rescale_intercept + elif method == 'fp': + slope = 1.0 / scale_slope + intercept = rescale_intercept / (rescale_slope * scale_slope) + else: + raise ValueError("Unknown scaling method '%s'." % method) + reorder = self.get_sorted_slice_indices() + slope = slope[reorder] + intercept = intercept[reorder] + shape = (1, 1) + self.get_data_shape()[2:] + slope = slope.reshape(shape, order='F') + intercept = intercept.reshape(shape, order='F') + return slope, intercept + + def get_slice_orientation(self): + """Returns the slice orientation label. + + Returns + ------- + orientation : {'transverse', 'sagittal', 'coronal'} + """ + orientations = list(self.image_defs['Slice Orientation']) + if len(set(orientations)) > 1: + raise XMLRECError( + 'Varying slice orientation found in the image sequence. This ' + ' is not suppported.') + return slice_orientation_rename[orientations[0].decode()] + + def get_rec_shape(self): + inplane_shape = self._get_unique_resolution() + return inplane_shape + (len(self.image_defs),) + + def _strict_sort_order(self): + """ Determine the sort order based on several image definition fields. + + The fields taken into consideration, if present, are (in order from + slowest to fastest variation after sorting): + + - image_defs['image_type_mr'] # Re, Im, Mag, Phase + - image_defs['dynamic scan number'] # repetition + - image_defs['label type'] # ASL tag/control + - image_defs['diffusion b value number'] # diffusion b value + - image_defs['gradient orientation number'] # diffusion directoin + - image_defs['cardiac phase number'] # cardiac phase + - image_defs['echo number'] # echo + - image_defs['slice number'] # slice + + Data sorting is done in two stages: + + 1. an initial sort using the keys described above + 2. a resort after generating two additional sort keys: + + * a key to assign unique volume numbers to any volumes that + didn't have a unique sort based on the keys above + (see :func:`vol_numbers`). + * a sort key based on `vol_is_full` to identify truncated + volumes + + A case where the initial sort may not create a unique label for each + volume is diffusion scans acquired in the older V4 .PAR format, where + diffusion direction info is not available. + """ + # sort keys present in all supported .PAR versions + idefs = self.image_defs + index_nos = idefs['Index'] + slice_nos = idefs['Slice'] + dynamics = idefs['Dynamic'] + phases = idefs['Phase'] + echos = idefs['Echo'] + image_type = idefs['Type'] + + # sort keys only present in a subset of .PAR files + asl_keys = (idefs['Label Type'], ) + diffusion_keys = (idefs['BValue'], idefs['Grad Orientation']) + + # initial sort (last key is highest precedence) + keys = (index_nos, slice_nos, echos, phases) + \ + diffusion_keys + asl_keys + (dynamics, image_type) + initial_sort_order = np.lexsort(keys) + + # sequentially number the volumes based on the initial sort + vol_nos = vol_numbers(slice_nos[initial_sort_order]) + # identify truncated volumes + is_full = vol_is_full(slice_nos[initial_sort_order], + self.general_info['Max No Slices']) + + # second stage of sorting + return initial_sort_order[np.lexsort((vol_nos, is_full))] + + def _lax_sort_order(self): + """ + Sorts by (fast to slow): slice number, volume number. + + We calculate volume number by looking for repeating slice numbers (see + :func:`vol_numbers`). + """ + slice_nos = self.image_defs['Slice'] + is_full = vol_is_full(slice_nos, self.general_info['Max No Slices']) + keys = (slice_nos, vol_numbers(slice_nos), np.logical_not(is_full)) + return np.lexsort(keys) + + def get_volume_labels(self): + """ Dynamic labels corresponding to the final data dimension(s). + + This is useful for custom data sorting. A subset of the info in + ``self.image_defs`` is returned in an order that matches the final + data dimension(s). Only labels that have more than one unique value + across the dataset will be returned. + + Returns + ------- + sort_info : dict + Each key corresponds to volume labels for a dynamically varying + sequence dimension. The ordering of the labels matches the volume + ordering determined via ``self.get_sorted_slice_indices``. + """ + sorted_indices = self.get_sorted_slice_indices() + image_defs = self.image_defs + + # define which keys which might vary across image volumes + dynamic_keys = ['Phase', + 'Echo', + 'Label Type', + 'Type', + 'Dynamic', + 'Sequence', + 'Grad Orient', + 'BValue'] + + # remove dynamic keys that may not be present in older .PAR versions + dynamic_keys = [d for d in dynamic_keys if d in + image_defs.dtype.fields] + + non_unique_keys = [] + for key in dynamic_keys: + ndim = image_defs[key].ndim + if ndim == 1: + num_unique = len(np.unique(image_defs[key])) + else: + raise ValueError("unexpected image_defs shape > 1D") + if num_unique > 1: + non_unique_keys.append(key) + + # each key in dynamic keys will be identical across slices, so use + # the value at slice 1. + sl1_indices = image_defs['Slice'][sorted_indices] == 1 + + sort_info = OrderedDict() + for key in non_unique_keys: + sort_info[key] = image_defs[key][sorted_indices][sl1_indices] + return sort_info + + class XMLRECImage(PARRECImage): """XML/REC image""" header_class = XMLRECHeader From 10bf7cca36690aea5c668e282e95ab710be2fd62 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Fri, 19 Oct 2018 13:15:13 -0400 Subject: [PATCH 11/16] remove permit_truncated restriction --- nibabel/cmdline/parrec2nii.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/nibabel/cmdline/parrec2nii.py b/nibabel/cmdline/parrec2nii.py index 7ce947cc99..7bc2b77826 100644 --- a/nibabel/cmdline/parrec2nii.py +++ b/nibabel/cmdline/parrec2nii.py @@ -154,9 +154,6 @@ def proc_file(infile, opts): if (ext == '.xml' or (ext == '.rec' and os.path.exists(basefilename + '.xml'))): pr_module = xr - if opts.permit_truncated: - raise ValueError("The permit_truncated option is not currently " - "supported for .xml/.REC data") else: pr_module = pr From 8d621199a43a6282784adb95df6f119ce2835091 Mon Sep 17 00:00:00 2001 From: Gregory Lee Date: Fri, 19 Oct 2018 14:12:17 -0400 Subject: [PATCH 12/16] minor cleanup --- nibabel/xmlrec.py | 77 +++++++++++++++++++++-------------------------- 1 file changed, 34 insertions(+), 43 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index b3e649ba50..6924096dc8 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -48,7 +48,7 @@ class XMLRECError(Exception): supported_xml_versions = ['PRIDE_V5'] -def _process_gen_dict_XML(xml_root): +def _process_gen_dict_xml(xml_root): """Read the general_info from an XML file. This is the equivalent of _process_gen_dict() for .PAR files @@ -79,7 +79,7 @@ def _get_image_def_attributes(xml_root): Paramters --------- xml_root : - The root of the XML tree + The root of the XML tree. Returns ------- @@ -141,7 +141,12 @@ def _process_image_lines_xml(xml_root): Parameters ---------- xml_root : - The root of the XML tree + The root of the XML tree. + + Returns + ------- + image_defs : np.ndarray + A structured array with the labels for each 2D image frame. """ image_defs_array = xml_root.find('Image_Array') if image_defs_array is None: @@ -165,7 +170,6 @@ def _get_val(entry_dtype, text): return val image_defs = np.zeros(len(image_defs_array), dtype=image_def_dtd) - already_warned = [] for i, image_def in enumerate(image_defs_array): if image_def.find('Key') != image_def[0]: @@ -208,12 +212,9 @@ def parse_XML_header(fobj): Structured array with fields giving all "Image information" in the header """ - # single pass through the header tree = ET.parse(fobj) root = tree.getroot() - # _split_header() equivalent - version = root.tag # e.g. PRIDE_V5 if version not in supported_xml_versions: warnings.warn(one_line( @@ -223,7 +224,7 @@ def parse_XML_header(fobj): adding support for this version. """.format(version))) try: - general_info = _process_gen_dict_XML(root) + general_info = _process_gen_dict_xml(root) image_defs = _process_image_lines_xml(root) except ET.ParseError: raise XMLRECError( @@ -235,7 +236,7 @@ def parse_XML_header(fobj): def _truncation_checks(general_info, image_defs, permit_truncated): - """ Check for presence of truncation in PAR file parameters + """ Check for presence of truncation in XML file parameters Raise error if truncation present and `permit_truncated` is False. """ @@ -268,7 +269,7 @@ def _chk_trunc(idef_name, gdef_max_name): class XMLRECHeader(PARRECHeader): - """PAR/REC header""" + """XML/REC header""" def __init__(self, info, image_defs, permit_truncated=False, strict_sort=False): @@ -276,11 +277,11 @@ def __init__(self, info, image_defs, permit_truncated=False, Parameters ---------- info : dict - "General information" from the PAR file (as returned by - `parse_PAR_header()`). + "General information" from the XML file (as returned by + `parse_XML_header()`). image_defs : array - Structured array with image definitions from the PAR file (as - returned by `parse_PAR_header()`). + Structured array with image definitions from the XML file (as + returned by `parse_XML_header()`). permit_truncated : bool, optional If True, a warning is emitted instead of an error when a truncated recording is detected. @@ -288,7 +289,7 @@ def __init__(self, info, image_defs, permit_truncated=False, If True, a larger number of header fields are used while sorting the REC data array. This may produce a different sort order than `strict_sort=False`, where volumes are sorted by the order in which - the slices appear in the .PAR file. + the slices appear in the .xml file. """ self.general_info = info.copy() self.image_defs = image_defs.copy() @@ -331,10 +332,10 @@ def copy(self): self.strict_sort) def as_analyze_map(self): - """Convert PAR parameters to NIFTI1 format""" + """Convert XML parameters to NIFTI1 format""" # Entries in the dict correspond to the parameters found in # the NIfTI1 header, specifically in nifti1.py `header_dtd` defs. - # Here we set the parameters we can to simplify PAR/REC + # Here we set the parameters we can to simplify XML/REC # to NIfTI conversion. exam_date = '{}/{}'.format( self.general_info['Examination Date'], @@ -390,10 +391,6 @@ def get_bvals_bvecs(self): bvecs = apply_affine(np.linalg.inv(permute_to_psl), bvecs) return bvals, bvecs - @deprecate_with_version('get_voxel_size deprecated. ' - 'Please use "get_zooms" instead.', - '2.0', '4.0') - def _get_unique_resolution(self): """Return the 2D image plane shape. @@ -406,6 +403,9 @@ def _get_unique_resolution(self): 'not suppported.') return (resx[0], resy[0]) + @deprecate_with_version('get_voxel_size deprecated. ' + 'Please use "get_zooms" instead.', + '2.0', '4.0') def get_voxel_size(self): """Returns the spatial extent of a voxel. @@ -426,7 +426,6 @@ def get_voxel_size(self): slice_thickness)) return voxsize - def _calc_zooms(self): """Compute image zooms from header data. @@ -482,8 +481,8 @@ def get_affine(self, origin='scanner'): Notes ----- Transformations appear to be specified in (ap, fh, rl) axes. The - orientation of data is recorded in the "slice orientation" field of the - PAR header "General Information". + orientation of data is recorded in the "Slice Orientation" field of the + XML header "General Information". We need to: @@ -579,7 +578,7 @@ def get_data_scaling(self, method="dv"): Notes ----- - The PAR header contains two different scaling settings: 'dv' (value on + The XML header contains two different scaling settings: 'dv' (value on console) and 'fp' (floating point value). Here is how they are defined: DV = PV * RS + RI @@ -635,14 +634,15 @@ def _strict_sort_order(self): The fields taken into consideration, if present, are (in order from slowest to fastest variation after sorting): - - image_defs['image_type_mr'] # Re, Im, Mag, Phase - - image_defs['dynamic scan number'] # repetition - - image_defs['label type'] # ASL tag/control - - image_defs['diffusion b value number'] # diffusion b value - - image_defs['gradient orientation number'] # diffusion directoin - - image_defs['cardiac phase number'] # cardiac phase - - image_defs['echo number'] # echo - - image_defs['slice number'] # slice + - image_defs['Type'] # Re, Im, Mag, Phase + - image_defs['Dynamicr'] # repetition + - image_defs['Label Type'] # ASL tag/control + - image_defs['BValue'] # diffusion b value + - image_defs['Grad Orientation'] # diffusion directoin + - image_defs['Phase'] # cardiac phase + - image_defs['Echo'] # echo + - image_defs['Slice'] # slice + - image_defs['Index'] # index in the REC file Data sorting is done in two stages: @@ -655,11 +655,8 @@ def _strict_sort_order(self): * a sort key based on `vol_is_full` to identify truncated volumes - A case where the initial sort may not create a unique label for each - volume is diffusion scans acquired in the older V4 .PAR format, where - diffusion direction info is not available. """ - # sort keys present in all supported .PAR versions + # sort keys present in all supported .xml versions idefs = self.image_defs index_nos = idefs['Index'] slice_nos = idefs['Slice'] @@ -667,8 +664,6 @@ def _strict_sort_order(self): phases = idefs['Phase'] echos = idefs['Echo'] image_type = idefs['Type'] - - # sort keys only present in a subset of .PAR files asl_keys = (idefs['Label Type'], ) diffusion_keys = (idefs['BValue'], idefs['Grad Orientation']) @@ -726,10 +721,6 @@ def get_volume_labels(self): 'Grad Orient', 'BValue'] - # remove dynamic keys that may not be present in older .PAR versions - dynamic_keys = [d for d in dynamic_keys if d in - image_defs.dtype.fields] - non_unique_keys = [] for key in dynamic_keys: ndim = image_defs[key].ndim From d7bf69325a927ef7bfcf37d8253e80824dcc7d0a Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 22 Oct 2018 14:57:22 -0400 Subject: [PATCH 13/16] XMLRECImage: fix for .REC expected case mismatch --- nibabel/xmlrec.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 6924096dc8..0c023df4f7 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -1,4 +1,5 @@ from copy import copy, deepcopy +import os import warnings import xml.etree.ElementTree as ET from collections import OrderedDict @@ -747,4 +748,15 @@ class XMLRECImage(PARRECImage): valid_exts = ('.rec', '.xml') files_types = (('image', '.rec'), ('header', '.xml')) + @classmethod + def filespec_to_file_map(klass, filespec): + file_map = super(PARRECImage, klass).filespec_to_file_map(filespec) + # fix case of .REC (.xml/.REC tends to have mixed case file extensions) + fname = file_map['image'].filename + fname_upper = fname.replace('.rec', '.REC') + if (not os.path.exists(fname) and + os.path.exists(fname_upper)): + file_map['image'].filename = fname_upper + return file_map + load = XMLRECImage.load From fd88ccd65530f389a7b6daa96b3147502bc4b0e4 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 22 Oct 2018 15:17:38 -0400 Subject: [PATCH 14/16] fix misnamed key within strict_sort --- nibabel/xmlrec.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 0c023df4f7..c9f73d8c99 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -639,7 +639,7 @@ def _strict_sort_order(self): - image_defs['Dynamicr'] # repetition - image_defs['Label Type'] # ASL tag/control - image_defs['BValue'] # diffusion b value - - image_defs['Grad Orientation'] # diffusion directoin + - image_defs['Grad Orient'] # diffusion directoin - image_defs['Phase'] # cardiac phase - image_defs['Echo'] # echo - image_defs['Slice'] # slice @@ -666,7 +666,7 @@ def _strict_sort_order(self): echos = idefs['Echo'] image_type = idefs['Type'] asl_keys = (idefs['Label Type'], ) - diffusion_keys = (idefs['BValue'], idefs['Grad Orientation']) + diffusion_keys = (idefs['BValue'], idefs['Grad Orient']) # initial sort (last key is highest precedence) keys = (index_nos, slice_nos, echos, phases) + \ From 7048640f2cd120584ba08359c46887c164973aca Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 22 Oct 2018 15:18:03 -0400 Subject: [PATCH 15/16] fix for conversion of DWI scans where only the trace was retained --- nibabel/xmlrec.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index c9f73d8c99..77eef950ec 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -374,7 +374,10 @@ def get_bvals_bvecs(self): if self.general_info['Diffusion'] == 0: return None, None reorder = self.get_sorted_slice_indices() - n_slices, n_vols = self.get_data_shape()[-2:] + if len(self.get_data_shape()) == 3: + n_slices, n_vols = self.get_data_shape()[-1], 1 + else: + n_slices, n_vols = self.get_data_shape()[-2:] bvals = self.image_defs['Diffusion B Factor'][reorder].reshape( (n_slices, n_vols), order='F') # All bvals within volume should be the same From 4068548b786c88528033e658edf63db6ba1d9fa2 Mon Sep 17 00:00:00 2001 From: "Gregory R. Lee" Date: Mon, 22 Oct 2018 15:53:43 -0400 Subject: [PATCH 16/16] fix bug in conversion of Boolean entries from string such as 'N' and 'Y' prior to this commit, both would get stored as true because bool('N') == True --- nibabel/xmlrec.py | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/nibabel/xmlrec.py b/nibabel/xmlrec.py index 77eef950ec..fb8655a6e7 100755 --- a/nibabel/xmlrec.py +++ b/nibabel/xmlrec.py @@ -28,6 +28,15 @@ class XMLRECError(Exception): # Dict for converting XML/REC types to appropriate python types # could convert the Enumeration strings to int, but enums_dict is incomplete + +def _xml_str_to_bool(s): + if isinstance(s, bytes): + s = s.decode() + if s == 'Y': + return True + else: + return False + # and the strings are more easily interpretable xml_type_dict = {'Float': np.float32, 'Double': np.float64, @@ -62,12 +71,15 @@ def _process_gen_dict_xml(xml_root): a = e.attrib if 'Name' in a: entry_type = xml_type_dict[a['Type']] - if entry_type in ['S16', 'S32', 'S128']: - entry_type = str + if entry_type == bool: + # convert string such as 'N' or 'Y' to boolean + cast = _xml_str_to_bool + else: + cast = entry_type if 'ArraySize' in a: - val = [entry_type(i) for i in e.text.strip().split()] + val = [cast(i) for i in e.text.strip().split()] else: - val = entry_type(e.text) + val = cast(e.text) general_info[a['Name']] = val return general_info @@ -166,6 +178,9 @@ def _get_val(entry_dtype, text): val = text[:32] elif entry_dtype == '|S128': val = text[:128] + elif entry_dtype == bool: + # convert string 'Y' or 'N' to boolean + val = _xml_str_to_bool(text) else: val = entry_dtype(text) return val @@ -371,7 +386,7 @@ def get_bvals_bvecs(self): Array of b vectors, shape (n_directions, 3), or None if not a diffusion acquisition. """ - if self.general_info['Diffusion'] == 0: + if not self.general_info['Diffusion']: return None, None reorder = self.get_sorted_slice_indices() if len(self.get_data_shape()) == 3: