From 01ea003020bff61cb5ed6598f37ef021d27e9197 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Sat, 2 Mar 2024 18:33:23 -0500 Subject: [PATCH 01/18] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index fcb7ca3..2089c33 100644 --- a/.gitignore +++ b/.gitignore @@ -8,3 +8,4 @@ __pycache__ .DS_Store *.png prototyping +.idea/ From 33efc63b9d8140fba085f0b3ae9382a2e300c04f Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 4 Mar 2024 20:41:57 -0500 Subject: [PATCH 02/18] Some cosmetic edits based on PEP 8 and grammar --- spec2nii/GE/ge_hdr_fields.py | 4 ++-- spec2nii/GE/ge_pfile.py | 32 ++++++++++++++++---------------- spec2nii/GE/ge_read_pfile.py | 8 ++++---- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/spec2nii/GE/ge_hdr_fields.py b/spec2nii/GE/ge_hdr_fields.py index f6ae330..c11497d 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 diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index d4472b1..9d80f4d 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,12 +83,12 @@ 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'): @@ -116,12 +116,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 +144,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 +152,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 +179,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) @@ -203,12 +203,12 @@ def _process_slaser(pfile): 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 @@ -231,12 +231,12 @@ def _process_gaba(pfile): 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 +270,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 +292,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 diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index bc1db11..69b3a6c 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 @@ -1200,7 +1200,7 @@ def read_data(self): sequence. This is currently a reimplementation of the Gannert/Osprey GELoad.m - function. Therefore the logic is unlike the other mappers. + 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. """ @@ -1297,7 +1297,7 @@ def read_data(self): self.raw_suppressed = self.raw_data[:, :, :, :, Y2, :] * Y1 * mult # 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. From fdb8ddf48f427270d0f581bfc129eb10397f484f Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 4 Mar 2024 20:50:45 -0500 Subject: [PATCH 03/18] Add cv24 variable that is important for advanced options in GE data --- spec2nii/GE/ge_read_pfile.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index 69b3a6c..f1b1c8c 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1215,7 +1215,9 @@ def read_data(self): refframes = self.hdr.rhi_user19 nreceivers = self.get_num_coils - dataWordSize = self.hdr.rhr_rh_point_size + dataWordSize = self.hdr.rhr_rh_point_size + + cv24 = self.hdr.rhi_user24 # MM: CV24 is a control variable (CV) that defines several advanced options not stored in the other CVs # Check if single voxel numVoxels = self.get_num_voxels @@ -1281,7 +1283,10 @@ 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) + if cv24 >= 16384: # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + 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,7 +1295,10 @@ 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) + if cv24 >= 16384: # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + 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) From fb49a56c7f75bc5a10b58031992487530363b61e Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Tue, 5 Mar 2024 16:09:23 -0500 Subject: [PATCH 04/18] Spelling corrections --- spec2nii/GE/ge_pfile.py | 2 +- spec2nii/GE/ge_read_pfile.py | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index 9d80f4d..a88ac26 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -321,7 +321,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 f1b1c8c..f88a38d 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -209,7 +209,7 @@ def get_mapper(self): # ARC : Added for Bergen jpress patch mapper = PfileMapperGaba elif psd == 'jpress': - # wtc - Added for HURCULES data. + # wtc - Added for J-edited data. mapper = PfileMapperGaba elif psd == 'fidall': # WTC - added for JG's Hyperpolarised 13C data @@ -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,7 +1199,7 @@ 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 + 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. @@ -1316,9 +1316,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)) From 2557ec589cce842c6601a77cb5998a6369e7fc8b Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Tue, 5 Mar 2024 21:37:12 -0500 Subject: [PATCH 05/18] Update ge_read_pfile.py Add an if-then statement for sLASER datasets where the receiver phase was set --- spec2nii/GE/ge_read_pfile.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index f88a38d..f7d7089 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1283,7 +1283,8 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(refframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - if cv24 >= 16384: # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + if cv24 >= 16384: Y1 = np.ones_like(X1) else: Y1 = (-1) ** (noadd * X1) @@ -1295,7 +1296,8 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(dataframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - if cv24 >= 16384: # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + if cv24 >= 16384: Y1 = np.ones_like(X1) else: Y1 = (-1) ** (noadd * X1) From 7ce2e2d287231a057d3b9e571d831b9f28b37aa4 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Tue, 5 Mar 2024 21:53:29 -0500 Subject: [PATCH 06/18] Refined handling of edited GE data --- spec2nii/GE/ge_pfile.py | 13 ++++++------- spec2nii/GE/ge_read_pfile.py | 9 ++++++--- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index a88ac26..9ff7333 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -90,18 +90,17 @@ def _process_svs_pfile(pfile): :return: List of file name suffixes """ psd = pfile.hdr.rhi_psdname.decode('utf-8').lower() + numecho = pfile.hdr.rhi_numecho if psd in ('probe-p', 'probe-s'): 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 == '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 + elif psd in 'oslaser' and numecho > 1: # MM: If edited data, use _process_gaba data, meta, dwelltime, fname_suffix = _process_gaba(pfile) - elif psd == 'jpress': + elif psd in 'slaser': + data, meta, dwelltime, fname_suffix = _process_slaser(pfile) + elif psd in ('jpress', 'jpress_ac', 'gaba', 'hbcd'): data, meta, dwelltime, fname_suffix = _process_gaba(pfile) else: raise UnsupportedPulseSequenceError(f'Unrecognised sequence {psd}.') diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index f7d7089..ae4d947 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -177,12 +177,15 @@ def get_mapper(self): return None psd = self.hdr.rhi_psdname.decode('utf-8').lower() + numecho = self.hdr.rhi_numecho if psd in ('probe-p', 'probe-s'): mapper = PfileMapper - elif psd in ('oslaser', 'slaser_cni', 'slaser'): - mapper = PfileMapperSlaser - elif psd == 'presscsi': + elif psd in ('oslaser', 'slaser_cni', 'slaser') and numecho == 1: + mapper = PfileMapperSlaser # MM: If non-edited data, use PfileMapperSlaser + elif psd in 'oslaser' and numecho > 1: + mapper = PfileMapperGaba # MM: If edited data, use PfileMapperGaba + elif psd in 'presscsi': mapper = PfileMapper elif psd == 'fidcsi': # bjs - added for Pom's fidcsi 13C data From 2bcb814e01dfb060a92d7f27ab5d1ed55a9e8d8e Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:29:59 -0400 Subject: [PATCH 07/18] Update `ge_hdr_fields.py` Added missing header fields for revisions 14, 15, and 16 --- spec2nii/GE/ge_hdr_fields.py | 30 ++++++++++++++++++++++++------ 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/spec2nii/GE/ge_hdr_fields.py b/spec2nii/GE/ge_hdr_fields.py index c11497d..ec24d08 100644 --- a/spec2nii/GE/ge_hdr_fields.py +++ b/spec2nii/GE/ge_hdr_fields.py @@ -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)) From 9ac5bf9986a4057f871cdd06fb2d5cb3940c0071 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:32:39 -0400 Subject: [PATCH 08/18] Update ge_read_pfile.py Cosmetic edits to code --- spec2nii/GE/ge_read_pfile.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index ae4d947..a928255 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -206,7 +206,7 @@ def get_mapper(self): # ATG - added HBCD - reuse GABA mapper 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 @@ -1220,12 +1220,13 @@ def read_data(self): nreceivers = self.get_num_coils dataWordSize = self.hdr.rhr_rh_point_size - cv24 = self.hdr.rhi_user24 # MM: CV24 is a control variable (CV) that defines several advanced options not stored in the other CVs + # 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 @@ -1259,7 +1260,7 @@ def read_data(self): 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):, :] @@ -1286,7 +1287,7 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(refframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + # 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: @@ -1299,7 +1300,7 @@ def read_data(self): X1, X2 = np.meshgrid(np.arange(dataframes), np.arange(nechoes)) X1 = X1.T.ravel() X2 = X2.T.ravel() - # Do not apply any phase cycling correction when the receiver phase toggle in sLASER has been set + # 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: From 7dd986c665fd24fef488413bb8e359eb3f1a7665 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:40:04 -0400 Subject: [PATCH 09/18] Better handling of GE psd strings - Cleaned up code a bit in `ge_read_pfile.py` - Added support for a few other psd's --- spec2nii/GE/ge_pfile.py | 11 +++++--- spec2nii/GE/ge_read_pfile.py | 55 ++++++++++++++++-------------------- 2 files changed, 32 insertions(+), 34 deletions(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index 9ff7333..578f390 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -90,17 +90,20 @@ def _process_svs_pfile(pfile): :return: List of file name suffixes """ psd = pfile.hdr.rhi_psdname.decode('utf-8').lower() + if psd.endswith('gaba'): # MM: some 'gaba' psd strings contain full path names, so truncate to the end of the path + psd = 'gaba' + numecho = pfile.hdr.rhi_numecho - if psd in ('probe-p', 'probe-s'): + 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') and numecho == 1: # MM: If non-edited data, use _process_oslaser data, meta, dwelltime, fname_suffix = _process_oslaser(pfile) - elif psd in 'oslaser' and numecho > 1: # MM: If edited data, use _process_gaba + elif psd == 'oslaser' and numecho > 1: # MM: If edited data, use _process_gaba data, meta, dwelltime, fname_suffix = _process_gaba(pfile) - elif psd in 'slaser': + elif psd == 'slaser': data, meta, dwelltime, fname_suffix = _process_slaser(pfile) - elif psd in ('jpress', 'jpress_ac', 'gaba', 'hbcd'): + 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}.') diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index a928255..d32506c 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -177,46 +177,41 @@ def get_mapper(self): return None psd = self.hdr.rhi_psdname.decode('utf-8').lower() + if psd.endswith('gaba'): + psd = 'gaba' + numecho = self.hdr.rhi_numecho - if psd in ('probe-p', 'probe-s'): + 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') and numecho == 1: mapper = PfileMapperSlaser # MM: If non-edited data, use PfileMapperSlaser - elif psd in 'oslaser' and numecho > 1: + elif psd == 'oslaser' and numecho > 1: mapper = PfileMapperGaba # MM: If edited data, use PfileMapperGaba - elif psd in '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 \ + ( + '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 mapper = PfileMapperProbeSL - elif 'jpress_ac' in psd: - # ARC : Added for Bergen jpress patch - mapper = PfileMapperGaba - elif psd == 'jpress': - # wtc - Added for J-edited 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) From f9dc7748dad3d7b72d73bd0ced8e637fb17de9d0 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:40:44 -0400 Subject: [PATCH 10/18] Add support of GE header revision 14.3 --- spec2nii/GE/ge_read_pfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index d32506c..c5054dc 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -309,7 +309,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, From 62b1158ad7d0b522507caca2d806baf604d15b59 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:43:50 -0400 Subject: [PATCH 11/18] Bug fix in `ge_read_pfile.py` The last block of lines was previously within the `else` statement just above, which prevented those lines from running when `nechoes == 1` --- spec2nii/GE/ge_read_pfile.py | 38 +++++++++++++++++++----------------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index c5054dc..6d9a942 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1305,21 +1305,23 @@ def read_data(self): self.raw_suppressed = self.raw_data[:, :, :, :, Y2, :] * Y1 * mult - # Up to this point we have simply replicated the logic of the GELoad function. - # 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): - reorg_suppressed.append(self.raw_suppressed[:, :, :, :, ne::nechoes, :]) - reorg_unsuppressed.append(self.raw_unsuppressed[:, :, :, :, ne::nechoes, :]) - - reorg_suppressed = np.stack(reorg_suppressed, axis=-1) - # 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) - # Rearrange axes to (x, y, z, t, coils, dynamics, edit) - self.raw_unsuppressed = np.moveaxis(reorg_unsuppressed, (4, 5), (5, 4)) + # Up to this point we have simply replicated the logic of the GELoad function. + # 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 = [] + if nechoes == 1: + nechoes = 2 + + for ne in range(nechoes): + reorg_suppressed.append(self.raw_suppressed[:, :, :, :, ne::nechoes, :]) + reorg_unsuppressed.append(self.raw_unsuppressed[:, :, :, :, ne::nechoes, :]) + + reorg_suppressed = np.stack(reorg_suppressed, axis=-1) + # 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) + # Rearrange axes to (x, y, z, t, coils, dynamics, edit) + self.raw_unsuppressed = np.moveaxis(reorg_unsuppressed, (4, 5), (5, 4)) From 284fe9a242556dec8f3acdf9bbe23e412e707a47 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 14:47:48 -0400 Subject: [PATCH 12/18] Bug fix in `ge_read_pfile.py` Incorrect header fields were being used to determine `dataframes` and `refframes` --- spec2nii/GE/ge_read_pfile.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index 6d9a942..0798d7b 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1209,8 +1209,8 @@ 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 From 88df41188d7daeeff742f64fc6ac9793d7b47f26 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Mon, 11 Mar 2024 15:05:47 -0400 Subject: [PATCH 13/18] Update ge_read_pfile.py - Updated the `mult` and `multw` factors to match the current ones used in Gannet [NOTE: These have still not been completely finalized and may need to be adjusted in the future] - Some cosmetic edits --- spec2nii/GE/ge_read_pfile.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index 0798d7b..fd8b89b 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1247,18 +1247,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 / 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: @@ -1268,7 +1267,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 From cffc9162a6dd24096ed5189065fae8ef93352a41 Mon Sep 17 00:00:00 2001 From: Mark Mikkelsen Date: Tue, 12 Mar 2024 11:36:58 -0400 Subject: [PATCH 14/18] Cosmetic edits --- spec2nii/GE/ge_pfile.py | 4 +++- spec2nii/GE/ge_read_pfile.py | 2 ++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index 578f390..121ecfc 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -90,7 +90,9 @@ def _process_svs_pfile(pfile): :return: List of file name suffixes """ psd = pfile.hdr.rhi_psdname.decode('utf-8').lower() - if psd.endswith('gaba'): # MM: some 'gaba' psd strings contain full path names, so truncate to the end of the path + + # 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 diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index fd8b89b..e119872 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -177,6 +177,8 @@ def get_mapper(self): return None psd = self.hdr.rhi_psdname.decode('utf-8').lower() + + # MM: Some 'gaba' psd strings contain full path names, so truncate to the end of the path if psd.endswith('gaba'): psd = 'gaba' From 82f535f98a55bb4315846b933e9f5e9899cf79df Mon Sep 17 00:00:00 2001 From: wtclarke Date: Wed, 20 Mar 2024 14:11:07 +0000 Subject: [PATCH 15/18] Add tests to MMs GE edits --- CHANGELOG.md | 4 + tests/io_for_tests.py | 12 +++ tests/spec2nii_test_data | 2 +- tests/test_ge_pfile.py | 214 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 231 insertions(+), 1 deletion(-) 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/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..91d2b5d 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,203 @@ 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'] + + 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'] + + 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] == 2000.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'] + + 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' + if hdr_ext['SequenceName'] == 'jpress': + assert img.shape[6] == 2 + 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'] + + 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 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' + if hdr_ext['SequenceName'] == 'jpress': + assert img.shape[6] == 2 + assert hdr_ext_ref['dim_7'] == 'DIM_EDIT' + + +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, 32, 2) + 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] == 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' From 7fb6d9257e53b41b529547750b70476f18f35570 Mon Sep 17 00:00:00 2001 From: wtclarke Date: Wed, 20 Mar 2024 14:28:03 +0000 Subject: [PATCH 16/18] Tweak test --- tests/test_ge_pfile.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_ge_pfile.py b/tests/test_ge_pfile.py index 91d2b5d..65d47f7 100644 --- a/tests/test_ge_pfile.py +++ b/tests/test_ge_pfile.py @@ -202,7 +202,7 @@ def test_mm_mega(tmp_path): assert img.shape[3] in (2048, 4096) assert img.shape[6] == 2 assert np.iscomplexobj(img.dataobj) - assert 1 / img.header['pixdim'][4] == 2000.0 + 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' From 07c85a4c16ab8c6c19c3084a4c1b016c90002c1b Mon Sep 17 00:00:00 2001 From: wtclarke Date: Wed, 20 Mar 2024 21:40:33 +0000 Subject: [PATCH 17/18] Fix for the editing disabled jpress case. --- spec2nii/GE/ge_pfile.py | 8 ++++-- spec2nii/GE/ge_read_pfile.py | 52 ++++++++++++++++++++++-------------- tests/test_ge_pfile.py | 11 ++------ 3 files changed, 40 insertions(+), 31 deletions(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index 121ecfc..7523ef1 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -225,11 +225,15 @@ 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: + meta.set_dim_info(2, 'DIM_EDIT') 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: + meta_ref.set_dim_info(2, 'DIM_EDIT') return [metab, water], [meta, meta_ref], dwelltime, ['', '_ref'] diff --git a/spec2nii/GE/ge_read_pfile.py b/spec2nii/GE/ge_read_pfile.py index e119872..77054eb 100644 --- a/spec2nii/GE/ge_read_pfile.py +++ b/spec2nii/GE/ge_read_pfile.py @@ -1310,23 +1310,35 @@ def read_data(self): self.raw_suppressed = self.raw_data[:, :, :, :, Y2, :] * Y1 * mult - # Up to this point we have simply replicated the logic of the GELoad function. - # 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 = [] - if nechoes == 1: - nechoes = 2 - - for ne in range(nechoes): - reorg_suppressed.append(self.raw_suppressed[:, :, :, :, ne::nechoes, :]) - reorg_unsuppressed.append(self.raw_unsuppressed[:, :, :, :, ne::nechoes, :]) - - reorg_suppressed = np.stack(reorg_suppressed, axis=-1) - # 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) - # Rearrange axes to (x, y, z, t, coils, dynamics, edit) - self.raw_unsuppressed = np.moveaxis(reorg_unsuppressed, (4, 5), (5, 4)) + # 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 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): + reorg_suppressed.append(self.raw_suppressed[:, :, :, :, ne::nechoes, :]) + reorg_unsuppressed.append(self.raw_unsuppressed[:, :, :, :, ne::nechoes, :]) + + reorg_suppressed = np.stack(reorg_suppressed, axis=-1) + # 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) + # 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/test_ge_pfile.py b/tests/test_ge_pfile.py index 65d47f7..56544cf 100644 --- a/tests/test_ge_pfile.py +++ b/tests/test_ge_pfile.py @@ -242,9 +242,6 @@ def test_mm_press(tmp_path): 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' - if hdr_ext['SequenceName'] == 'jpress': - assert img.shape[6] == 2 - 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' @@ -254,15 +251,12 @@ def test_mm_press(tmp_path): 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 number of coils + # 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' - if hdr_ext['SequenceName'] == 'jpress': - assert img.shape[6] == 2 - assert hdr_ext_ref['dim_7'] == 'DIM_EDIT' def test_mm_press_noid(tmp_path): @@ -276,12 +270,11 @@ def test_mm_press_noid(tmp_path): 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, 32, 2) + 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 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] == mm_press_noid.name From 072d1d19ebf177f6fd2dfeee45c21b824f01ed37 Mon Sep 17 00:00:00 2001 From: wtclarke Date: Thu, 21 Mar 2024 15:32:11 +0000 Subject: [PATCH 18/18] Add MEGA editing pulse information to the header. --- spec2nii/GE/ge_pfile.py | 37 +++++++++++++++++++++++++++++++++++-- tests/test_ge_pfile.py | 11 +++++++++++ 2 files changed, 46 insertions(+), 2 deletions(-) diff --git a/spec2nii/GE/ge_pfile.py b/spec2nii/GE/ge_pfile.py index 7523ef1..acd42a8 100644 --- a/spec2nii/GE/ge_pfile.py +++ b/spec2nii/GE/ge_pfile.py @@ -206,6 +206,39 @@ 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 @@ -227,13 +260,13 @@ def _process_gaba(pfile): meta.set_dim_info(1, 'DIM_DYN') # Only set an EDIT dim if there is an editing dimension if metab.ndim == 7: - meta.set_dim_info(2, 'DIM_EDIT') + _add_editing_info(pfile, meta, metab) meta_ref.set_dim_info(0, 'DIM_COIL') meta_ref.set_dim_info(1, 'DIM_DYN') # Only set an EDIT dim if there is an editing dimension if water.ndim == 7: - meta_ref.set_dim_info(2, 'DIM_EDIT') + _add_editing_info(pfile, meta_ref, water) return [metab, water], [meta, meta_ref], dwelltime, ['', '_ref'] diff --git a/tests/test_ge_pfile.py b/tests/test_ge_pfile.py index 56544cf..fe37e66 100644 --- a/tests/test_ge_pfile.py +++ b/tests/test_ge_pfile.py @@ -143,6 +143,10 @@ def test_mm_hercules(tmp_path): 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) @@ -176,6 +180,9 @@ def test_mm_hermes(tmp_path): 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) @@ -211,6 +218,10 @@ def test_mm_mega(tmp_path): 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)