diff --git a/.gitignore b/.gitignore index fcb7ca3..2089c33 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ __pycache__ .DS_Store *.png prototyping +.idea/ diff --git a/CHANGELOG.md b/CHANGELOG.md index 53d8458..35a2321 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,9 @@ This document contains the Spec2nii release history in reverse chronological order. +0.7.4 (WIP) +--------------------------------- +- Refinements and improvements to the GE SVS pipeline from Mark Mikkelsen + 0.7.3 (Tuesday 12th March 2024) ------------------------------- - Siemens .rda format now had corrected and validated orientations (tested on VE11 baseline). diff --git a/spec2nii/GE/ge_hdr_fields.py b/spec2nii/GE/ge_hdr_fields.py index f6ae330..ec24d08 100644 --- a/spec2nii/GE/ge_hdr_fields.py +++ b/spec2nii/GE/ge_hdr_fields.py @@ -1,4 +1,4 @@ -'''List of GE p-file offsets. +"""List of GE p-file offsets. This code is taken from the VESPA project https://scion.duhs.duke.edu/vespa/project. I therefore include their BSD statement here. @@ -60,7 +60,7 @@ ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' +""" import ctypes as ct import numpy as np @@ -754,10 +754,16 @@ def get_pfile_hdr_fields(version): plist.append(('rhr_rh_npasses', ct.c_short)) plist.append(('pad_xx', ct.c_char * 2)) plist.append(('rhr_rh_nslices', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 10)) + plist.append(('rhr_rh_nechoes', ct.c_short)) + plist.append(('rhr_rh_navs', ct.c_short)) + plist.append(('rhr_rh_nframes', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 4)) plist.append(('rhr_rh_frame_size', ct.c_ushort)) plist.append(('rhr_rh_point_size', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 32)) + plist.append(('pad_xx', ct.c_char * 18)) + plist.append(('rhr_rh_da_xres', ct.c_short)) + plist.append(('rhr_rh_da_yres', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 10)) plist.append(('rhr_rh_raw_pass_size', ct.c_uint)) plist.append(('pad_xx', ct.c_char * 80)) plist.append(('rhr_rh_dab[0]_start_rcv', ct.c_short)) @@ -978,10 +984,16 @@ def get_pfile_hdr_fields(version): plist.append(('rhr_rh_npasses', ct.c_short)) plist.append(('pad_xx', ct.c_char * 2)) plist.append(('rhr_rh_nslices', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 10)) + plist.append(('rhr_rh_nechoes', ct.c_short)) + plist.append(('rhr_rh_navs', ct.c_short)) + plist.append(('rhr_rh_nframes', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 4)) plist.append(('rhr_rh_frame_size', ct.c_ushort)) plist.append(('rhr_rh_point_size', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 32)) + plist.append(('pad_xx', ct.c_char * 18)) + plist.append(('rhr_rh_da_xres', ct.c_short)) + plist.append(('rhr_rh_da_yres', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 10)) plist.append(('rhr_rh_raw_pass_size', ct.c_uint)) plist.append(('pad_xx', ct.c_char * 80)) plist.append(('rhr_rh_dab[0]_start_rcv', ct.c_short)) @@ -1202,10 +1214,16 @@ def get_pfile_hdr_fields(version): plist.append(('rhr_rh_npasses', ct.c_short)) plist.append(('pad_xx', ct.c_char * 2)) plist.append(('rhr_rh_nslices', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 10)) + plist.append(('rhr_rh_nechoes', ct.c_short)) + plist.append(('rhr_rh_navs', ct.c_short)) + plist.append(('rhr_rh_nframes', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 4)) plist.append(('rhr_rh_frame_size', ct.c_ushort)) plist.append(('rhr_rh_point_size', ct.c_short)) - plist.append(('pad_xx', ct.c_char * 32)) + plist.append(('pad_xx', ct.c_char * 18)) + plist.append(('rhr_rh_da_xres', ct.c_short)) + plist.append(('rhr_rh_da_yres', ct.c_short)) + plist.append(('pad_xx', ct.c_char * 10)) plist.append(('rhr_rh_raw_pass_size', ct.c_uint)) plist.append(('pad_xx', ct.c_char * 80)) plist.append(('rhr_rh_dab[0]_start_rcv', ct.c_short)) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index d4472b1..acd42a8 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -1,4 +1,4 @@ -'''spec2nii module containing functions specific to interpreting the GE pfile format +"""spec2nii module containing functions specific to interpreting the GE pfile format Author: William Clarke Copyright (C) 2020 University of Oxford @@ -33,7 +33,7 @@ For portions of this code, copyright and license information differs from the above. In these cases, copyright and/or license information is inline. -''' +""" # Python modules from datetime import datetime @@ -83,25 +83,29 @@ def read_pfile(filename, fname_out): def _process_svs_pfile(pfile): - '''Handle SVS data + """Handle SVS data :param Pfile pfile: Pfile object :return: List of NIFTI MRS data objects :return: List of file name suffixes - ''' + """ psd = pfile.hdr.rhi_psdname.decode('utf-8').lower() - if psd in ('probe-p', 'probe-s'): + # MM: Some 'gaba' psd strings contain full path names, so truncate to the end of the path + if psd.endswith('gaba'): + psd = 'gaba' + + numecho = pfile.hdr.rhi_numecho + + if psd in ('probe-p', 'probe-s', 'probe-p_ach'): data, meta, dwelltime, fname_suffix = _process_probe_p(pfile) - elif psd in ('oslaser', 'slaser_cni'): + elif psd in ('oslaser', 'slaser_cni') and numecho == 1: # MM: If non-edited data, use _process_oslaser data, meta, dwelltime, fname_suffix = _process_oslaser(pfile) + elif psd == 'oslaser' and numecho > 1: # MM: If edited data, use _process_gaba + data, meta, dwelltime, fname_suffix = _process_gaba(pfile) elif psd == 'slaser': data, meta, dwelltime, fname_suffix = _process_slaser(pfile) - elif psd in ('gaba', 'hbcd'): - data, meta, dwelltime, fname_suffix = _process_gaba(pfile) - elif 'jpress_ac' in psd: # Bergen patch - data, meta, dwelltime, fname_suffix = _process_gaba(pfile) - elif psd == 'jpress': + elif psd in ('jpress', 'jpress_ac', 'gaba', 'hbcd', 'probe-p-mega_rml', 'repress7'): data, meta, dwelltime, fname_suffix = _process_gaba(pfile) else: raise UnsupportedPulseSequenceError(f'Unrecognised sequence {psd}.') @@ -116,12 +120,12 @@ def _process_svs_pfile(pfile): def _process_probe_p(pfile): - '''Extract metabolite and reference data from a prob_p format pfile + """Extract metabolite and reference data from a prob_p format pfile :param Pfile pfile: Pfile object :return: List numpy data arrays :return: List of file name suffixes - ''' + """ metab = pfile.map.raw_suppressed # typically (1,1,1,navg,ncoil,npts) metab = np.transpose(metab, [0, 1, 2, 5, 4, 3]) # swap to (1,1,1,npts,ncoil,navg) @@ -144,7 +148,7 @@ def _process_probe_p(pfile): def _process_oslaser(pfile): - '''Extract metabolite and reference data from a oslaser format pfile + """Extract metabolite and reference data from a oslaser format pfile I think this is like the CMRR sLASER sequence with 2 ecc and 2 water scaling scans at the start and end of each acquisition. @@ -152,7 +156,7 @@ def _process_oslaser(pfile): :param Pfile pfile: Pfile object :return: List numpy data arrays :return: List of file name suffixes - ''' + """ water = pfile.map.raw_unsuppressed # typically (1,1,1,navg,ncoil,npts) metab = pfile.map.raw_suppressed @@ -179,14 +183,14 @@ def _process_oslaser(pfile): def _process_slaser(pfile): - '''Extract metabolite and reference data from a slaser format pfile + """Extract metabolite and reference data from a slaser format pfile This seems to be like a standard probe-p. Maybe slaser is the canonical vendor implementation. :param Pfile pfile: Pfile object :return: List numpy data arrays :return: List of file name suffixes - ''' + """ metab = pfile.map.raw_suppressed # typically (1,1,1,navg,ncoil,npts) metab = np.transpose(metab, [0, 1, 2, 5, 4, 3]) # swap to (1,1,1,npts,ncoil,navg) @@ -202,13 +206,46 @@ def _process_slaser(pfile): return [metab, water], [meta, meta_ref], dwelltime, ['', '_ref'] +def _add_editing_info(pfile, meta, data): + """Add editing information to dimension tags and headers + + :param pfile: p-file object + :type pfile: Pfile + :param meta: Header extension object + :type meta: Hdr_Ext + :param data: Shaped complex data + :type data: np.ndarray + """ + edit_rf_waveform = pfile.hdr.rhi_user19 + # edit_rf_waveform == 19.0 is used by HERMES and HERCULES + if data.shape[-1] == 2 and not edit_rf_waveform == 19.0: + edit_rf_freq_off1 = pfile.hdr.rhi_user20 + edit_rf_freq_off2 = pfile.hdr.rhi_user21 + edit_rf_ppm_off1 = edit_rf_freq_off1 / float(pfile.hdr.rhr_rh_ps_mps_freq * 1E-7) + edit_rf_ppm_off2 = edit_rf_freq_off2 / float(pfile.hdr.rhr_rh_ps_mps_freq * 1E-7) + edit_rf_dur = pfile.hdr.rhi_user22 + # check for default value (-1) of pulse length + if edit_rf_dur <= 0: + edit_rf_dur = 16000 + dim_info = "MEGA-EDITED j-difference editing, two conditions" + dim_header = {"EditCondition": ["ON", "OFF"]} + edit_pulse_val = { + "ON": {"PulseOffset": edit_rf_ppm_off1, "PulseDuration": edit_rf_dur / 1E6}, + "OFF": {"PulseOffset": edit_rf_ppm_off2, "PulseDuration": edit_rf_dur / 1E6}} + + meta.set_dim_info(2, 'DIM_EDIT', hdr=dim_header, info=dim_info) + meta.set_standard_def("EditPulse", edit_pulse_val) + else: + meta.set_dim_info(2, 'DIM_EDIT') + + def _process_gaba(pfile): - '''Extract metabolite and reference data from a gaba (MPRESS) format pfile + """Extract metabolite and reference data from a gaba (MPRESS) format pfile :param Pfile pfile: Pfile object :return: List numpy data arrays :return: List of file name suffixes - ''' + """ # Note that custom mapper sorts dimensions already metab = pfile.map.raw_suppressed @@ -221,22 +258,26 @@ def _process_gaba(pfile): meta.set_dim_info(0, 'DIM_COIL') meta.set_dim_info(1, 'DIM_DYN') - meta.set_dim_info(2, 'DIM_EDIT') + # Only set an EDIT dim if there is an editing dimension + if metab.ndim == 7: + _add_editing_info(pfile, meta, metab) meta_ref.set_dim_info(0, 'DIM_COIL') meta_ref.set_dim_info(1, 'DIM_DYN') - meta_ref.set_dim_info(2, 'DIM_EDIT') + # Only set an EDIT dim if there is an editing dimension + if water.ndim == 7: + _add_editing_info(pfile, meta_ref, water) return [metab, water], [meta, meta_ref], dwelltime, ['', '_ref'] def _process_mrsi_pfile(pfile): - '''Handle MRSI data + """Handle MRSI data :param Pfile pfile: Pfile object :return: List of NIFTI MRS data objects :return: List of file name suffixes - ''' + """ psd = pfile.hdr.rhi_psdname.decode('utf-8').lower() known_formats = ('probe-p', 'probe-sl', 'slaser_cni', 'presscsi') @@ -270,7 +311,7 @@ def fft_and_shift(x, axis): def _calculate_affine_mrsi(pfile): - '''Calculate the 4x4 affine matrix for mrsi''' + """Calculate the 4x4 affine matrix for mrsi""" dcos = pfile.map.get_dcos.T dcos[dcos == 0.0] = 0.0 # remove -0.0 values @@ -292,7 +333,7 @@ def _calculate_affine_mrsi(pfile): def _calculate_affine(pfile): - '''Calculate the 4x4 affine matrix''' + """Calculate the 4x4 affine matrix""" dcos = pfile.map.get_dcos.T dcos[dcos == 0.0] = 0.0 # remove -0.0 values @@ -321,7 +362,7 @@ def _populate_metadata(pfile, water_suppressed=True, data_dimensions=None): :type pfile: pfile map object :param water_suppressed: Set water suppression header field, defaults to True :type water_suppressed: bool, optional - :param data_dimensions: If set to 5,6, or 7 will inlcude default dim tags for those diemnsions, defaults to None + :param data_dimensions: If set to 5,6, or 7 will include default dim tags for those dimensions, defaults to None :type data_dimensions: int, optional :return: Header extension object :rtype: nifti_mrs.hdr_ext diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index bc1db11..77054eb 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1,4 +1,4 @@ -'''Reader for GE p-files. +"""Reader for GE p-files. This code is taken from the VESPA project https://scion.duhs.duke.edu/vespa/project. I therefore include their BSD statement here. @@ -60,7 +60,7 @@ ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -''' +""" # Python modules import math import sys @@ -178,42 +178,42 @@ def get_mapper(self): psd = self.hdr.rhi_psdname.decode('utf-8').lower() - if psd in ('probe-p', 'probe-s'): + # MM: Some 'gaba' psd strings contain full path names, so truncate to the end of the path + if psd.endswith('gaba'): + psd = 'gaba' + + numecho = self.hdr.rhi_numecho + + if psd in \ + ( + 'probe-p', + 'probe-s', + 'probe-p_ach', # MM - added for Calgary PROBE-P sequence + 'presscsi', + 'fidcsi', # bjs - added for Pom's fidcsi 13C data + 'ia/stable/fidcsi', # bjs - added for Kearny's 13C data + 'presscsi_nfl', # bjs - added for Govind's SVS data off v25 + 'epsi_3d_24', # bjs - added for soher check of MIDAS Browndyke data + 'fidall' # WTC - added for JG's Hyperpolarised 13C data + ): mapper = PfileMapper - elif psd in ('oslaser', 'slaser_cni', 'slaser'): - mapper = PfileMapperSlaser - elif psd == 'presscsi': - mapper = PfileMapper - elif psd == 'fidcsi': - # bjs - added for Pom's fidcsi 13C data - mapper = PfileMapper - elif psd == 'ia/stable/fidcsi': - # bjs - added for Kearny's 13C data - mapper = PfileMapper - elif psd == 'presscsi_nfl': - # bjs - added for Govind's SVS data off v25 - mapper = PfileMapper - elif psd == 'epsi_3d_24': - # bjs - added for soher check of MIDAS Browndyke data - mapper = PfileMapper - elif psd == 'gaba': - # wtc - added for Nottingham MEGA-PRESS sequence. - mapper = PfileMapperGaba - elif psd == 'hbcd': - # ATG - added HBCD - reuse GABA mapper + elif psd in ('oslaser', 'slaser_cni', 'slaser') and numecho == 1: + mapper = PfileMapperSlaser # MM: If non-edited data, use PfileMapperSlaser + elif psd == 'oslaser' and numecho > 1: + mapper = PfileMapperGaba # MM: If edited data, use PfileMapperGaba + elif psd in \ + ( + 'jpress', # wtc - added for J-edited data. + 'jpress_ac', # ARC - added for Bergen jpress patch + 'gaba', # wtc - added for Nottingham MEGA-PRESS sequence + 'hbcd', # ATG - added HBCD - reuse GABA mapper + 'probe-p-mega_rml', # MM - added for Calgary MEGA-PRESS sequence + 'repress7' # MM - added for old CUBRIC data + ): mapper = PfileMapperGaba elif psd == 'probe-sl': - # wtc - added for CSI sequence from Manchester. + # wtc - added for CSI sequence from Manchester mapper = PfileMapperProbeSL - elif 'jpress_ac' in psd: - # ARC : Added for Bergen jpress patch - mapper = PfileMapperGaba - elif psd == 'jpress': - # wtc - Added for HURCULES data. - mapper = PfileMapperGaba - elif psd == 'fidall': - # WTC - added for JG's Hyperpolarised 13C data - mapper = PfileMapper else: raise UnknownPfile("No Pfile mapper for pulse sequence = %s" % psd) @@ -311,7 +311,7 @@ def _version(self, filelike): known_revisions = \ [ 7, 8, 9, 10, 11, - 14.0, 14.1, 14.2, + 14.0, 14.1, 14.2, 14.3, 15.000, 15.001, 16.000, 20.001, 20.002, 20.003, 20.004, 20.005, 20.006, 20.007, 24.000, @@ -711,7 +711,7 @@ def is_chop_on(self): @property def get_frequency_offset(self): - """ Returns the spectral frquency offset """ + """ Returns the spectral frequency offset """ if self.version > 9: return 0.0 else: @@ -825,8 +825,8 @@ def get_num_voxels_in_vol(self): def get_num_kspace_points(self): """ Determine the number of sampled k-space points in the data set. - This may differ from the number of voxels in the rectalinear grid, - for example if elliptical or another non rectangular acquisition + This may differ from the number of voxels in the rectilinear grid, + for example if elliptical or another non-rectangular acquisition sampling strategy was employed. GE product sequences pad the reduced k-space data with zeros so the number of k-space points is the same as the number of voxels, but that may not be true for @@ -1079,7 +1079,7 @@ def read_data(self): self.raw_data = data # Modify the data loading behavior. For single voxel multi-acq data - # this means return the averaged (suppresssed data, if applicable). + # this means return the averaged (suppressed data, if applicable). numUnsuppressed = self.get_number_unsuppressed_acquisitions numSuppressed = self.get_number_suppressed_acquisitions @@ -1186,10 +1186,10 @@ def __init__(self, file_name, hdr, version, endian): Without some more knowledgeable input I can't currently square the logic of BJS's mapper classes with the logic from Gannet/Osprey etc. - Therefore this is really a hack that uses the PfileMapper orientation logic + Therefore, this is really a hack that uses the PfileMapper orientation logic but otherwise uses the logic from Gannet/Osprey. - Thus most work is done in the overloaded read_data method and functions like + Thus, most work is done in the overloaded read_data method and functions like get_num_voxels_in_vol are currently meaningless. I would like to square both methods in the future. """ @@ -1199,8 +1199,8 @@ def read_data(self): """Function that contains all the data loading logic for the 'gaba' sequence. - This is currently a reimplementation of the Gannert/Osprey GELoad.m - function. Therefore the logic is unlike the other mappers. + This is currently a reimplementation of the Gannet/Osprey GELoad.m + function. Therefore, the logic is unlike the other mappers. The suppressed and unsuppressed data can be fetched from the raw_suppressed and raw_unsuppressed property. """ @@ -1211,16 +1211,19 @@ def read_data(self): npoints = self.hdr.rhr_rh_da_xres nrows = self.hdr.rhr_rh_da_yres - dataframes = self.hdr.rhi_user4 / nex - refframes = self.hdr.rhi_user19 + dataframes = self.hdr.rhr_rh_user4 / nex + refframes = int(self.hdr.rhr_rh_user19) nreceivers = self.get_num_coils - dataWordSize = self.hdr.rhr_rh_point_size + dataWordSize = self.hdr.rhr_rh_point_size + + # MM: CV24 is a control variable (CV) that defines several advanced options not stored in the other CVs + cv24 = self.hdr.rhi_user24 # Check if single voxel numVoxels = self.get_num_voxels self.is_svs = False - if (numVoxels[0] * numVoxels[1] * numVoxels[2] == 1): + if numVoxels[0] * numVoxels[1] * numVoxels[2] == 1: self.is_svs = True # Data loading @@ -1246,18 +1249,17 @@ def read_data(self): tempData = tempData.view(np.complex64) if nechoes == 1: - tempData = tempData.reshape((1, 1, 1, nreceivers, totalframes, npoints)) - self.raw_data = np.swapaxes(tempData, -1, -3) - - if (dataframes + refframes) != nframes: + if int(dataframes + refframes) != nframes: mult = 1 dataframes *= nex refframes = int(nframes - dataframes) else: - mult = 1.0 / nex + mult = 1 / nex - self.raw_unsuppressed = self.raw_data[:, :, :, :, 1:(refframes + 1), :] - self.raw_suppressed = self.raw_data[:, :, :, :, (refframes + 1):, :] + tempData = tempData.reshape((1, 1, 1, nreceivers, totalframes, npoints)) + self.raw_data = np.swapaxes(tempData, -1, -3) + self.raw_unsuppressed = self.raw_data[:, :, :, :, 1:(refframes + 1), :] * mult + self.raw_suppressed = self.raw_data[:, :, :, :, (refframes + 1):, :] * mult / 2 else: if int(dataframes + refframes) != nframes: @@ -1267,7 +1269,11 @@ def read_data(self): dataframes *= nex refframes = nframes - dataframes else: - mult = nex / 2 + # MM: Change mult to match latest code in Gannet (2020) + # Previous factors: + # mult: 1 (RTN 2017), nex / 2 (RTN 2016) + # multw: 1 / nex (RTN 2017), 1 (RTN 2016) + mult = 1 / nex multw = 1 noadd = 0 @@ -1281,7 +1287,11 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(refframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - Y1 = (-1)**(noadd * X1) + # Do not apply any phase cycling correction if the receiver phase toggle in sLASER was set + if cv24 >= 16384: + Y1 = np.ones_like(X1) + else: + Y1 = (-1) ** (noadd * X1) Y1 = np.moveaxis(np.broadcast_to(Y1, (1, 1, 1, npoints, nreceivers, Y1.size)), -1, -2) Y2 = (1 + (totalframes / nechoes) * X2 + X1).astype(int) @@ -1290,17 +1300,35 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(dataframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - Y1 = (-1)**(noadd * X1) + # Do not apply any phase cycling correction if the receiver phase toggle in sLASER was set + if cv24 >= 16384: + Y1 = np.ones_like(X1) + else: + Y1 = (-1) ** (noadd * X1) Y1 = np.moveaxis(np.broadcast_to(Y1, (1, 1, 1, npoints, nreceivers, Y1.size)), -1, -2) Y2 = (1 + refframes + (totalframes / nechoes) * X2 + X1).astype(int) self.raw_suppressed = self.raw_data[:, :, :, :, Y2, :] * Y1 * mult + # Test if this is the case of the jpress like sequence being used for a non-editing condition + edit_waveform = self.hdr.rhi_user19 + if nechoes == 1 and edit_waveform == 0: + # Editing conditions = 1 and there is no editing waveform + self.raw_suppressed[:, :, :, :, 0::2, :] *= np.exp(1j * np.pi) + self.raw_unsuppressed[:, :, :, :, 0::2, :] *= np.exp(1j * np.pi) + + # Rearrange axes to (x, y, z, t, coils, dynamics) + self.raw_suppressed = np.moveaxis(self.raw_suppressed, (4, 5), (5, 4)) + self.raw_unsuppressed = np.moveaxis(self.raw_unsuppressed, (4, 5), (5, 4)) + else: + # Editing condition, update nechoes if needed for special cases like probe-p-mega_rml + if nechoes == 1: + nechoes = 2 + # Up to this point we have simply replicated the logic of the GELoad function. - # Now reorganise dimensions to give a editing dimension. + # Now reorganise dimensions to give an editing dimension. # This means that this is done in a not particularly clear order, but it enables testing against # the matlab code. - reorg_suppressed = [] reorg_unsuppressed = [] for ne in range(nechoes): @@ -1308,9 +1336,9 @@ def read_data(self): reorg_unsuppressed.append(self.raw_unsuppressed[:, :, :, :, ne::nechoes, :]) reorg_suppressed = np.stack(reorg_suppressed, axis=-1) - # Rearange axes to (x, y, z, t, coils, dynamics, edit) + # Rearrange axes to (x, y, z, t, coils, dynamics, edit) self.raw_suppressed = np.moveaxis(reorg_suppressed, (4, 5), (5, 4)) reorg_unsuppressed = np.stack(reorg_unsuppressed, axis=-1) - # Rearange axes to (x, y, z, t, coils, dynamics, edit) + # Rearrange axes to (x, y, z, t, coils, dynamics, edit) self.raw_unsuppressed = np.moveaxis(reorg_unsuppressed, (4, 5), (5, 4)) diff --git a/tests/io_for_tests.py b/tests/io_for_tests.py index dcdd842..f546e2f 100644 --- a/tests/io_for_tests.py +++ b/tests/io_for_tests.py @@ -1,4 +1,5 @@ import nibabel as nib +import json def read_nifti_mrs(file_path): @@ -7,3 +8,14 @@ def read_nifti_mrs(file_path): img = nib.load(str(file_path)) return img + + +def read_nifti_mrs_with_hdr(file_path): + '''Read NIFTI MRS files using nibabel. + Return img and header extracted.''' + + img = nib.load(str(file_path)) + hdr_ext_codes = img.header.extensions.get_codes() + hdr_ext = json.loads(img.header.extensions[hdr_ext_codes.index(44)].get_content()) + + return img, hdr_ext diff --git a/tests/spec2nii_test_data b/tests/spec2nii_test_data index 74415ef..2181603 160000 --- a/tests/spec2nii_test_data +++ b/tests/spec2nii_test_data @@ -1 +1 @@ -Subproject commit 74415ef6bb9556384cd685e766c0e4b47998902b +Subproject commit 2181603035ef4647a338ad4950a21aaaa7c9929c diff --git a/tests/test_ge_pfile.py b/tests/test_ge_pfile.py index 793631c..fe37e66 100644 --- a/tests/test_ge_pfile.py +++ b/tests/test_ge_pfile.py @@ -11,6 +11,8 @@ import numpy as np from .io_for_tests import read_nifti_mrs +from .io_for_tests import read_nifti_mrs_with_hdr + # Data paths ge_path = Path(__file__).parent / 'spec2nii_test_data' / 'ge' @@ -18,6 +20,18 @@ mrsi_path = ge_path / 'pFiles' / 'MRSI' / 'P18432.7' mpress_path = ge_path / 'pFiles' / 'big_gaba' / 'S01_GABA_68.7' +# Test set from Mark Mikkelsen +mm_herc = [ + ge_path / 'pFiles' / 'HERCULES' / 'sub-01_ses-01_acq-herculespress_svs_noID.7', + ge_path / 'pFiles' / 'HERCULES' / 'sub-01_ses-01_acq-herculesslaser_svs_noID.7'] +mm_hermes = [ + ge_path / 'pFiles' / 'HERMES' / 'sub-01_ses-01_acq-hermespress_svs_noID.7', + ge_path / 'pFiles' / 'HERMES' / 'sub-01_ses-01_acq-hermesslaser_svs_noID.7'] +mm_mega = ge_path / 'pFiles' / 'MEGA' / 'Big_GABA' +mm_press = ge_path / 'pFiles' / 'PRESS' / 'Big_GABA' +mm_press_noid = ge_path / 'pFiles' / 'PRESS' / 'sub-01_ses-01_acq-press_svs_noID.7' +mm_slaser = ge_path / 'pFiles' / 'sLASER' / 'sub-01_ses-01_acq-slaser_svs_noID.7' + def test_svs(tmp_path): @@ -104,3 +118,207 @@ def test_mpress(tmp_path): assert np.isclose(hdr_ext['EchoTime'], 0.068) assert np.isclose(hdr_ext['RepetitionTime'], 2.0) + + +def test_mm_hercules(tmp_path): + for path in mm_herc: + subprocess.run([ + 'spec2nii', 'ge', + '-f', path.stem, + '-o', tmp_path, + '-j', + path]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / path.with_suffix('.nii.gz').name) + + assert img.shape == (1, 1, 1, 4096, 32, 56, 4) + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] == 5000.0 + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + assert hdr_ext['dim_7'] == 'DIM_EDIT' + assert np.isclose(127.771, hdr_ext['SpectrometerFrequency'][0], atol=1E-3) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == path.name + assert hdr_ext['WaterSuppressed'] + + # Test that the MEGA editing header info is not added + assert 'dim_7_header' not in hdr_ext + assert 'EditPulse' not in hdr_ext + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{path.stem}_ref.nii.gz') + assert img_ref.shape[:5] == (1, 1, 1, 4096, 32) + assert img_ref.shape[6] == 4 + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL' + assert hdr_ext_ref['dim_6'] == 'DIM_DYN' + assert hdr_ext_ref['dim_7'] == 'DIM_EDIT' + + +def test_mm_hermes(tmp_path): + for path in mm_hermes: + subprocess.run([ + 'spec2nii', 'ge', + '-f', path.stem, + '-o', tmp_path, + '-j', + path]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / path.with_suffix('.nii.gz').name) + + assert img.shape == (1, 1, 1, 4096, 32, 56, 4) + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] == 5000.0 + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + assert hdr_ext['dim_7'] == 'DIM_EDIT' + assert np.isclose(127.771, hdr_ext['SpectrometerFrequency'][0], atol=1E-3) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == path.name + assert hdr_ext['WaterSuppressed'] + + assert 'dim_7_header' not in hdr_ext + assert 'EditPulse' not in hdr_ext + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{path.stem}_ref.nii.gz') + assert img_ref.shape[:5] == (1, 1, 1, 4096, 32) + assert img_ref.shape[6] == 4 + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL' + assert hdr_ext_ref['dim_6'] == 'DIM_DYN' + assert hdr_ext_ref['dim_7'] == 'DIM_EDIT' + + +def test_mm_mega(tmp_path): + for path in mm_mega.rglob('*.7'): + subprocess.run([ + 'spec2nii', 'ge', + '-f', path.stem, + '-o', tmp_path, + '-j', + path]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / path.with_suffix('.nii.gz').name) + + assert img.shape[:3] == (1, 1, 1) + assert img.shape[3] in (2048, 4096) + assert img.shape[6] == 2 + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] in (2000.0, 5000.0) + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + assert hdr_ext['dim_7'] == 'DIM_EDIT' + assert np.isclose(127.7, hdr_ext['SpectrometerFrequency'][0], atol=1E-1) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == path.name + assert hdr_ext['WaterSuppressed'] + + assert 'dim_7_header' in hdr_ext + assert 'EditPulse' in hdr_ext + assert 'dim_7_info' in hdr_ext + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{path.stem}_ref.nii.gz') + assert img_ref.shape[:3] == (1, 1, 1) + assert img_ref.shape[6] == 2 + # Check reference has same number of points and coils as main without checking number of coils + assert img_ref.shape[3] == img.shape[3] + assert img_ref.shape[4] == img.shape[4] + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL' + assert hdr_ext_ref['dim_6'] == 'DIM_DYN' + assert hdr_ext_ref['dim_7'] == 'DIM_EDIT' + + +def test_mm_press(tmp_path): + for path in mm_press.rglob('*.7'): + subprocess.run([ + 'spec2nii', 'ge', + '-f', path.stem, + '-o', tmp_path, + '-j', + path]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / path.with_suffix('.nii.gz').name) + + assert img.shape[:3] == (1, 1, 1) + assert img.shape[3] in (2048, 4096) + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] in (2000.0, 5000.0) + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + + assert np.isclose(127.7, hdr_ext['SpectrometerFrequency'][0], atol=1E-1) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == path.name + assert hdr_ext['WaterSuppressed'] + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{path.stem}_ref.nii.gz') + assert img_ref.shape[:3] == (1, 1, 1) + # Check reference has same number of points and coils as main without checking actual number + assert img_ref.shape[3] == img.shape[3] + assert img_ref.shape[4] == img.shape[4] + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL' + assert hdr_ext_ref['dim_6'] == 'DIM_DYN' + + +def test_mm_press_noid(tmp_path): + subprocess.run([ + 'spec2nii', 'ge', + '-f', mm_press_noid.stem, + '-o', tmp_path, + '-j', + mm_press_noid]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / mm_press_noid.with_suffix('.nii.gz').name) + + assert img.shape == (1, 1, 1, 4096, 32, 64) + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] == 5000.0 + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + assert np.isclose(127.771, hdr_ext['SpectrometerFrequency'][0], atol=1E-3) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == mm_press_noid.name + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{mm_press_noid.stem}_ref.nii.gz') + assert img_ref.shape[:5] == (1, 1, 1, 4096, 32) + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL' + + +def test_mm_slaser(tmp_path): + subprocess.run([ + 'spec2nii', 'ge', + '-f', mm_slaser.stem, + '-o', tmp_path, + '-j', + mm_slaser]) + + img, hdr_ext = read_nifti_mrs_with_hdr( + tmp_path / mm_slaser.with_suffix('.nii.gz').name) + + assert img.shape == (1, 1, 1, 4096, 32, 64) + assert np.iscomplexobj(img.dataobj) + assert 1 / img.header['pixdim'][4] == 5000.0 + assert hdr_ext['dim_5'] == 'DIM_COIL' + assert hdr_ext['dim_6'] == 'DIM_DYN' + assert np.isclose(127.771, hdr_ext['SpectrometerFrequency'][0], atol=1E-3) + assert hdr_ext['ResonantNucleus'][0] == '1H' + assert hdr_ext['OriginalFile'][0] == mm_slaser.name + + img_ref, hdr_ext_ref = read_nifti_mrs_with_hdr( + tmp_path / f'{mm_slaser.stem}_ref.nii.gz') + assert img_ref.shape[:5] == (1, 1, 1, 4096, 32) + assert not hdr_ext_ref['WaterSuppressed'] + assert hdr_ext_ref['dim_5'] == 'DIM_COIL'