Twix MRSI pathway - Generate empty NIfTI-MRS (#69)

* Enable empty file generation for twix MRSI.

* Lint

* Update fsleyes orientation render commands.

CHANGELOG.md

@@ -1,5 +1,9 @@
 This document contains the Spec2nii release history in reverse chronological order.
 
+0.6.6 (Thursday 9th March 2023)
+-------------------------------
+- Added in ability to generate empty NIfTI-MRS for twix-pathway MRSI scans.
+
 0.6.5 (Wednesday 8th February 2023)
 -----------------------------------
 - Fixed bug in philips and rda metadata translation functions.

README.md

@@ -39,9 +39,9 @@ Miniconda can be installed by following the instructions on the [Conda website](
 ## Currently supported formats
 This table lists the currently supported formats. I have very limited experience with Philips and GE formats. Please get in touch if you are willing to help add to this list and/or supply validation data.
 
-| Format        | File extension | SVS | CSI | Automatic orientation |
+| Format        | File extension | SVS | MRSI| Automatic orientation |
 |---------------|----------------|-----|-----|-----------------------|
-| Siemens Twix  | .dat           | Yes | No  | Yes                   |
+| Siemens Twix  | .dat           | Yes | †   | Yes                   |
 | Siemens DICOM | .ima / .dcm    | Yes | Yes | Yes                   |
 | Siemens RDA   | .rda           | Yes | No  | Yes (WIP)             |
 | Philips       | .SPAR/.SDAT    | Yes | No  | Yes                   |
@@ -57,6 +57,8 @@ This table lists the currently supported formats. I have very limited experience
 | jMRUI         | .mrui          | Yes | No  | No                    |
 | ASCII         | .txt           | Yes | No  | No                    |
 
+† Partial handling - see section on Twix pathway for MRSI handling.
+
 ## Instructions
 spec2nii is called on the command line, and the conversion file type is specified with a subcommand.
 
@@ -79,6 +81,8 @@ Call with the -e flag to specify which MDH flag to convert. e.g.
 
 Twix format loop variables (e.g. `Ave` or `ida`) can be assigned to specific NIfTI dimensions using the `-d{5,6,7}` command line options. NIfTI MRS dimension tags (e.g. `DIM_COIL`) can be specified using the `-t{5,6,7}` command line options.
 
+As `spec2nii` is __not__ a reconstruction program, it cannot convert MRSI data. Far too little information is held in the twix headers to reconstruct arbitrary k,t-space data. However, if passed a file containing MRSI data `spec2nii` will attempt to create an empty NIfTI-MRS file with the correct orientation information, data shape, and header information. This empty file can then have data inserted from an offline reconstruction routine.
+
 ### Siemens DICOM
 `spec2nii dicom DCM_FILE_or_DIR`
 

build_process.md

@@ -1,4 +1,9 @@
-# Manual build process
+# Build Process
+Package build and upload to Pypi is handled by `.github/workflows/publish.yml`. This is triggered by publication of a new tagged release on the [Github releases page](https://github.com/wtclarke/spec2nii/releases). Upload access to Pypi is via stored secret.
+
+[Conda-forge](https://github.com/conda-forge/spec2nii-feedstock) then picks up the new Pypi package automatically.
+
+## Old manual build process
 1. Commit, push and run CI tests
 2. git tag -m VX.X.X X.X.X
 3. git push github master --tags

tests/test_ge_pfile_orientation.py

@@ -133,7 +133,8 @@ def test_svs_orientation(tmp_path):
                                '-xc', '0', '0', '-yc', '0', '0', '-zc', '0', '0',
                                '-hc', ge_path / 'from_dicom' / 'T1.nii.gz',
                                '-dr', '-211', '7400',
-                               tmp_path / 'svs.nii.gz', '-a', '50', '-cm', 'blue'])
+                               tmp_path / 'svs.nii.gz', '-ot', 'complex',
+                               '-a', '50', '-cm', 'blue'])
 
         fsl_ss = Image.open(tmp_path / f'svs_{idx}.png')
         fsl_ss_cropped = crop_and_flip_first_third(fsl_ss)
@@ -176,7 +177,8 @@ def test_mrsi_orientation(tmp_path):
                                '-hc', ge_path / 'from_dicom' / 'T1.nii.gz',
                                '-dr', '-211', '7400',
                                str(dcm), '-a', '50', '-cm', 'red',
-                               tmp_path / 'mrsi.nii.gz', '-a', '50', '-cm', 'blue'])
+                               tmp_path / 'mrsi.nii.gz', '-ot', 'complex',
+                               '-a', '50', '-cm', 'blue'])
 
         fsl_ss = Image.open(tmp_path / f'mrsi_{idx}.png')
         fsl_ss_cropped = crop_and_flip_first_third(fsl_ss)

tests/test_philips_dicom_orientation.py

@@ -83,7 +83,8 @@ def test_svs_orientation(tmp_path):
                                '-vl', str(pos[0]), str(pos[1]), str(pos[2]),
                                '-xc', '0', '0', '-yc', '0', '0', '-zc', '0', '0',
                                '-hc', structural_data,
-                               tmp_path / 'svs.nii.gz', '-a', '50', '-cm', 'blue'])
+                               tmp_path / 'svs.nii.gz', '-ot', 'complex',
+                               '-a', '50', '-cm', 'blue'])
 
         fsl_ss = Image.open(tmp_path / f'svs_{idx}.png')
         width, height = fsl_ss.size

tests/test_philips_sdat_spar_orientation.py

@@ -81,7 +81,8 @@ def test_svs_orientation(tmp_path):
                                '-vl', str(pos[0]), str(pos[1]), str(pos[2]),
                                '-xc', '0', '0', '-yc', '0', '0', '-zc', '0', '0',
                                '-hc', p_path / 'T1.nii.gz', '-dr', '-8200', '204600',
-                               tmp_path / 'svs.nii.gz', '-a', '50', '-cm', 'blue'])
+                               tmp_path / 'svs.nii.gz', '-ot', 'complex',
+                               '-a', '50', '-cm', 'blue'])
 
         fsl_ss = Image.open(tmp_path / f'svs_{idx}.png')
         width, height = fsl_ss.size

tests/test_siemens_mrsi_orientation.py

@@ -98,11 +98,12 @@ def test_VB(tmp_path):
         else:
             x, y, z = '95', '88', '95'
 
-        subprocess.check_call(['pythonw', '/Users/wclarke/opt/miniconda3/envs/fsl_mrs/bin/fsleyes',
+        subprocess.check_call(['fsleyes',
                                'render', '-of', op.join(tmp_path, f'csi_{idx}.png'),
                                '-vl', x, y, z,
                                '-hc', op.join(vb_path, 'T1.nii.gz'),
-                               op.join(tmp_path, name + '_d.nii.gz'), '-a', '60', '-cm', 'red'])
+                               op.join(tmp_path, name + '_d.nii.gz'), '-ot', 'complex',
+                               '-a', '60', '-cm', 'red'])
 
         fsl_ss = Image.open(op.join(tmp_path, f'csi_{idx}.png'))
         width, height = fsl_ss.size
@@ -142,11 +143,12 @@ def test_VE(tmp_path):
         # Make fsleyes rendering
         x, y, z = '57', '63', '37'
 
-        subprocess.check_call(['pythonw', '/Users/wclarke/opt/miniconda3/envs/fsl_mrs/bin/fsleyes',
+        subprocess.check_call(['fsleyes',
                                'render', '-of', op.join(tmp_path, f'csi_{idx}.png'),
                                '-vl', x, y, z,
                                '-hc', op.join(ve_path, 'T1.nii.gz'),
-                               op.join(tmp_path, name + '_d.nii.gz'), '-a', '60', '-cm', 'red'])
+                               op.join(tmp_path, name + '_d.nii.gz'), '-ot', 'complex',
+                               '-a', '60', '-cm', 'red'])
 
         fsl_ss = Image.open(op.join(tmp_path, f'csi_{idx}.png'))
         width, height = fsl_ss.size

tests/test_siemens_svs_orientation.py

@@ -93,12 +93,14 @@ def test_VB(tmp_path):
                                '-j', op.join(vb_path, f_d)])
 
         # Make fsleyes rendering
-        subprocess.check_call(['pythonw', '/Users/wclarke/opt/miniconda3/envs/fsl_mrs/bin/fsleyes',
+        subprocess.check_call(['fsleyes',
                                'render', '-of', op.join(tmp_path, f'svs_{idx}.png'),
                                '-vl', '95', '89', '96',
                                '-hc', op.join(vb_path, 'T1.nii.gz'),
-                               op.join(tmp_path, name + '_t.nii.gz'), '-a', '50', '-cm', 'red',
-                               op.join(tmp_path, name + '_d.nii.gz'), '-a', '50', '-cm', 'blue'])
+                               op.join(tmp_path, name + '_t.nii.gz'),
+                               '-ot', 'complex', '-a', '50', '-cm', 'red',
+                               op.join(tmp_path, name + '_d.nii.gz'),
+                               '-ot', 'complex', '-a', '50', '-cm', 'blue'])
 
         img_t = read_nifti_mrs(op.join(tmp_path, name + '_t.nii.gz'))
         img_d = read_nifti_mrs(op.join(tmp_path, name + '_d.nii.gz'))
@@ -148,12 +150,14 @@ def test_VE(tmp_path):
             x, y, z = '57', '48', '37'
         else:
             x, y, z = '53', '67', '56'
-        subprocess.check_call(['pythonw', '/Users/wclarke/opt/miniconda3/envs/fsl_mrs/bin/fsleyes',
+        subprocess.check_call(['fsleyes',
                                'render', '-of', op.join(tmp_path, f'svs_{idx}.png'),
                                '-vl', x, y, z,
                                '-hc', op.join(ve_path, 'T1.nii.gz'),
-                               op.join(tmp_path, name + '_t.nii.gz'), '-a', '50', '-cm', 'red',
-                               op.join(tmp_path, name + '_d.nii.gz'), '-a', '50', '-cm', 'blue'])
+                               op.join(tmp_path, name + '_t.nii.gz'),
+                               '-ot', 'complex', '-a', '50', '-cm', 'red',
+                               op.join(tmp_path, name + '_d.nii.gz'),
+                               '-ot', 'complex', '-a', '50', '-cm', 'blue'])
 
         img_t = read_nifti_mrs(op.join(tmp_path, name + '_t.nii.gz'))
         img_d = read_nifti_mrs(op.join(tmp_path, name + '_d.nii.gz'))

tests/test_twix.py

@@ -200,3 +200,55 @@ def test_XA_HERCULES(tmp_path):
     assert hdr_ext['dim_5'] == 'DIM_COIL'
     assert hdr_ext['dim_6'] == 'DIM_EDIT'
     assert hdr_ext['dim_7'] == 'DIM_DYN'
+
+
+def test_twix_mrsi_orientation(tmp_path):
+    '''Test that the (empty) mrsi has information matching the twix equivalent.'''
+
+    def run_comparison(file_twix, file_dcm):
+        subprocess.run(
+            ['spec2nii', 'dicom',
+             '-o', tmp_path,
+             '-f', 'dicom',
+             file_dcm])
+
+        subprocess.run(
+            ['spec2nii', 'twix',
+             '-e', 'image',
+             '-o', tmp_path,
+             '-f', 'twix',
+             file_twix])
+
+        img_d = read_nifti_mrs(tmp_path / 'dicom.nii.gz')
+        img_t = read_nifti_mrs(tmp_path / 'twix_empty.nii.gz')
+
+        assert np.allclose(img_t.affine, img_d.affine, atol=1E-3)
+        assert img_t.shape == img_d.shape
+
+    vb_pairs = [
+        ['meas_MID145_csi_se_Tra_sat_FID108735.dat', 'csi_se_Tra_sat_9_1'],
+        ['meas_MID143_csi_se_3D_C_S23_5_T20_3_10_FID108733.dat', 'csi_se_3D_C>S23.5>T20.3_10_8_1'],
+        ['meas_MID141_csi_se_3D_S_T23_5_C20_3_10_FID108731.dat', 'csi_se_3D_S>T23.5>C20.3_10_7_1'],
+        ['meas_MID139_csi_se_3D_T_C23_5_S20_3_10_FID108729.dat', 'csi_se_3D_T>C23.5>S20.3_10_6_1'],
+        ['meas_MID137_csi_se_C_S23_5_T20_3_10_FID108727.dat', 'csi_se_C>S23.5>T20.3_10_5_1'],
+        ['meas_MID135_csi_se_S_T23_5_C20_3_10_FID108725.dat', 'csi_se_S>T23.5>C20.3_10_4_1'],
+        ['meas_MID133_csi_se_T_C23_5_S20_3_10_FID108723.dat', 'csi_se_T>C23.5>S20.3_10_3_1']]
+
+    for pair in vb_pairs:
+        twix = siemens_path / 'VBData' / 'Twix' / pair[0]
+        dcm = siemens_path / 'VBData' / 'DICOM' / pair[1]
+        run_comparison(twix, dcm)
+
+    ve_pairs = [
+        ['meas_MID00223_FID62728_csi_se_t_c23_5_s20_3_R10.dat', 'csi_se_t>c23.5>s20.3_R10_4_1'],
+        ['meas_MID00225_FID62730_csi_se_s_t23_5_c20_3_R10.dat', 'csi_se_s>t23.5>c20.3_R10_5_1'],
+        ['meas_MID00227_FID62732_csi_se_c_s23_5_t20_3_R10.dat', 'csi_se_c>s23.5>t20.3_R10_6_1'],
+        ['meas_MID00229_FID62734_csi_se_3D_t_c23_5_s20_3_R10.dat', 'csi_se_3D_t>c23.5>s20.3_R10_7_1'],
+        ['meas_MID00231_FID62736_csi_se_3D_s_t23_5_c20_3_R10.dat', 'csi_se_3D_s>t23.5>c20.3_R10_8_1'],
+        ['meas_MID00233_FID62738_csi_se_3D_c_s23_5_t20_3_R10.dat', 'csi_se_3D_c>s23.5>t20.3_R10_9_1'],
+        ['meas_MID00244_FID62749_csi_se_iso_tra_sat.dat', 'csi_se_iso_tra_sat_14_1']]
+
+    for pair in ve_pairs:
+        twix = siemens_path / 'VEData' / 'Twix' / pair[0]
+        dcm = siemens_path / 'VEData' / 'DICOM' / pair[1]
+        run_comparison(twix, dcm)

tests/test_uih_orientation.py

@@ -61,7 +61,8 @@ def test_uih_svs(tmp_path):
                            '-of', tmp_path / 'svs.png',
                            '-vl', '245', '296', '10',
                            '-hc', data_base / 'mrs_data/2D_tra.nii.gz',
-                           tmp_path / 'svs.nii.gz', '-a', '50', '-cm', 'blue'])
+                           tmp_path / 'svs.nii.gz', '-ot', 'complex',
+                           '-a', '50', '-cm', 'blues'])
 
     fsl_render = flip_first_third(Image.open(tmp_path / 'svs.png'))
     screenshot = Image.open(screenshots['svs'])
@@ -90,7 +91,8 @@ def test_uih_2d_csi(tmp_path):
                            '-of', tmp_path / 'csi_2d.png',
                            '-vl', '245', '296', '9',
                            '-hc', data_base / 'mrs_data/2D_tra.nii.gz',
-                           tmp_path / 'csi_2d.nii.gz', '-a', '50', '-cm', 'blue', '-dr', '-0.9', '9.7'])
+                           tmp_path / 'csi_2d.nii.gz', '-ot', 'complex',
+                           '-a', '50', '-cm', 'blue', '-dr', '-0.9', '9.7'])
 
     fsl_render = flip_first_third(Image.open(tmp_path / 'csi_2d.png'))
     screenshot = Image.open(screenshots['csi_2d'])
@@ -119,7 +121,8 @@ def test_uih_3d_csi(tmp_path):
                            '-of', tmp_path / 'csi_3d.png',
                            '-vl', '245', '296', '9',
                            '-hc', data_base / 'mrs_3d/3d_tra.nii.gz',
-                           tmp_path / 'csi_3d.nii.gz', '-a', '50', '-cm', 'blue', '-dr', '-0.9', '9.7'])
+                           tmp_path / 'csi_3d.nii.gz', '-ot', 'complex',
+                           '-a', '50', '-cm', 'blue', '-dr', '-0.9', '9.7'])
 
     fsl_render = flip_first_third(Image.open(tmp_path / 'csi_3d.png'))
     screenshot = Image.open(screenshots['csi_3d'])