From a12538cc95c7dcbb4a3601b18fd2b48599ca0bb6 Mon Sep 17 00:00:00 2001 From: William T Clarke Date: Thu, 9 Mar 2023 22:20:38 +0000 Subject: [PATCH] Twix MRSI pathway - Generate empty NIfTI-MRS (#69) * Enable empty file generation for twix MRSI. * Lint * Update fsleyes orientation render commands. --- CHANGELOG.md | 4 + README.md | 8 +- build_process.md | 7 +- spec2nii/Siemens/twixfunctions.py | 237 +++++++++++++++++--- tests/test_ge_pfile_orientation.py | 6 +- tests/test_philips_dicom_orientation.py | 3 +- tests/test_philips_sdat_spar_orientation.py | 3 +- tests/test_siemens_mrsi_orientation.py | 10 +- tests/test_siemens_svs_orientation.py | 16 +- tests/test_twix.py | 52 +++++ tests/test_uih_orientation.py | 9 +- 11 files changed, 298 insertions(+), 57 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15145a4..9a91454 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/README.md b/README.md index 1e19f0c..c1e8313 100644 --- a/README.md +++ b/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` diff --git a/build_process.md b/build_process.md index 7a90ac5..77db110 100644 --- a/build_process.md +++ b/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 diff --git a/spec2nii/Siemens/twixfunctions.py b/spec2nii/Siemens/twixfunctions.py index f71f0bb..471f621 100644 --- a/spec2nii/Siemens/twixfunctions.py +++ b/spec2nii/Siemens/twixfunctions.py @@ -83,17 +83,16 @@ def is_product(): return xa_or_vx(hdr) == 'xa' and is_product() -def process_twix(twixObj, base_name_out, name_in, dataKey, dim_overrides, quiet=False, verbose=False, remove_os=False): - """Process a twix file. Identify type of MRS and then pass to the relevant function.""" - if seq_name(twixObj.hdr) == 'fid': +def id_sequence_type(twixhdr): + if seq_name(twixhdr) == 'fid': n_voxels = 0 - elif twixObj.hdr.Meas.lFinalMatrixSizePhase \ - and twixObj.hdr.Meas.lFinalMatrixSizeRead: - n_voxels = twixObj.hdr.Meas.lFinalMatrixSizeSlice \ - * twixObj.hdr.Meas.lFinalMatrixSizePhase \ - * twixObj.hdr.Meas.lFinalMatrixSizeRead - elif twixObj.hdr.Meas.lFinalMatrixSizeSlice: - n_voxels = twixObj.hdr.Meas.lFinalMatrixSizeSlice + elif twixhdr.Meas.lFinalMatrixSizePhase \ + and twixhdr.Meas.lFinalMatrixSizeRead: + n_voxels = twixhdr.Meas.lFinalMatrixSizeSlice \ + * twixhdr.Meas.lFinalMatrixSizePhase \ + * twixhdr.Meas.lFinalMatrixSizeRead + elif twixhdr.Meas.lFinalMatrixSizeSlice: + n_voxels = twixhdr.Meas.lFinalMatrixSizeSlice else: # If lFinalMatrixSize{Slice,Phase,Read} are all empty # Either unlocalised or unusually filled in headers. @@ -102,8 +101,19 @@ def process_twix(twixObj, base_name_out, name_in, dataKey, dim_overrides, quiet= n_voxels = 1 if n_voxels > 1: - return process_mrsi(twixObj, base_name_out, name_in, dataKey, quiet=quiet, verbose=verbose) + return 'mrsi' elif n_voxels == 0: + return 'fid' + elif n_voxels == 1: + return 'svs' + + +def process_twix(twixObj, base_name_out, name_in, dataKey, dim_overrides, quiet=False, verbose=False, remove_os=False): + """Process a twix file. Identify type of MRS and then pass to the relevant function.""" + + if id_sequence_type(twixObj.hdr) == 'mrsi': + return process_mrsi(twixObj, base_name_out, name_in, dataKey, quiet=quiet, verbose=verbose) + elif id_sequence_type(twixObj.hdr) == 'fid': return process_fid( twixObj, base_name_out, @@ -113,7 +123,7 @@ def process_twix(twixObj, base_name_out, name_in, dataKey, dim_overrides, quiet= remove_os, quiet=quiet, verbose=verbose) - else: + elif id_sequence_type(twixObj.hdr) == 'svs': return process_svs( twixObj, base_name_out, @@ -125,9 +135,77 @@ def process_twix(twixObj, base_name_out, name_in, dataKey, dim_overrides, quiet= verbose=verbose) -def process_mrsi(twixObj, base_name_out, name_in, dataKey, quiet=False, verbose=False): +def process_mrsi(twixObj, base_name_out, name_in, dataKey, remove_os=True, quiet=False, verbose=False): """Identify correct MRSI pathway, either simple internal reconstruction or to ismrmrd""" - raise NotImplementedError('MRSI pathway not yet implemented.') + + if not quiet: + print('Generating empty NIfTI-MRS file for MRSI input.') + if verbose: + print('Spec2nii is not a reconstruction program.') + print( + 'Insufficient information is stored in the twix/.dat header to carry out ' + 'generalised reconstruction of arbitrary k,t-space data.') + print( + 'This routine generates an empty file, ' + 'so that data reconstructed through another ' + 'pathway can be inserted into it.') + + # Data size + data_size = [] + if twixObj.hdr.Meas.lFinalMatrixSizePhase: + data_size.append(int(twixObj.hdr.Meas.lFinalMatrixSizePhase)) + else: + data_size.append(int(1)) + if twixObj.hdr.Meas.lFinalMatrixSizeRead: + data_size.append(int(twixObj.hdr.Meas.lFinalMatrixSizeRead)) + else: + data_size.append(int(1)) + if twixObj.hdr.Meas.lFinalMatrixSizeSlice: + data_size.append(int(twixObj.hdr.Meas.lFinalMatrixSizeSlice)) + else: + data_size.append(int(1)) + + data_size.append(int(twixObj.hdr.Meas.lVectorSize)) + + # Perform Orientation calculations + # 1) Calculate dicom like imageOrientationPatient, imagePositionPatient, pixelSpacing and slicethickness + orient = twix2DCMOrientation(twixObj['hdr'], verbose=verbose) + imageOrientationPatient, imagePositionPatient, pixelSpacing, slicethickness, dim_swapped = orient + + # 2) In the style of dcm2niix calculate the affine matrix + orientation = dcm_to_nifti_orientation(imageOrientationPatient, + imagePositionPatient, + np.append(pixelSpacing, slicethickness), + data_size[:3], + half_shift=True, + verbose=verbose) + + # Account for the first two dimensions being swapped if main orientation is Transverse + # See CSIOrientations function for logic + if dim_swapped: + data_size = [data_size[1], data_size[0]] + data_size[2:] + data = np.zeros(data_size, dtype=complex) + + # Extract dwellTime + dwellTime = twixObj['hdr']['MeasYaps'][('sRXSPEC', 'alDwellTime', '0')] / 1E9 + if remove_os: + dwellTime *= 2 + + # Extract metadata + meta_obj = extractTwixMetadata(twixObj['hdr'], basename(twixObj[dataKey].filename)) + + # Identify what those indices are + # If cha is one: loop over 3rd and higher dims and make 2D images + # If cha isn't present one: loop over 2nd and higher dims and make 1D images + # Don't write here, just fill up class property lists for later writing + if base_name_out: + out_str = base_name_out + else: + out_str = name_in.split('.')[0] + + out_str += '_empty' + + return [gen_nifti_mrs_hdr_ext(data, dwellTime, meta_obj, orientation.Q44), ], [out_str, ] def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_os, quiet=False, verbose=False): @@ -155,7 +233,7 @@ def process_svs(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # Perform Orientation calculations # 1) Calculate dicom like imageOrientationPatient,imagePositionPatient,pixelSpacing and slicethickness orient = twix2DCMOrientation(twixObj['hdr'], force_svs=True, verbose=verbose) - imageOrientationPatient, imagePositionPatient, pixelSpacing, slicethickness = orient + imageOrientationPatient, imagePositionPatient, pixelSpacing, slicethickness, _ = orient # 2) In the style of dcm2niix calculate the affine matrix orientation = dcm_to_nifti_orientation(imageOrientationPatient, @@ -357,7 +435,7 @@ def process_fid(twixObj, base_name_out, name_in, dataKey, dim_overrides, remove_ # Perform Orientation calculations # 1) Calculate dicom like imageOrientationPatient,imagePositionPatient,pixelSpacing and slicethickness orient = twix2DCMOrientation(twixObj['hdr'], force_svs=True, verbose=verbose) - imageOrientationPatient, imagePositionPatient, pixelSpacing, slicethickness = orient + imageOrientationPatient, imagePositionPatient, pixelSpacing, slicethickness, _ = orient # 2) In the style of dcm2niix calculate the affine matrix orientation = dcm_to_nifti_orientation(imageOrientationPatient, @@ -560,32 +638,35 @@ def twix2DCMOrientation(mapVBVDHdr, force_svs=False, verbose=False): # Orientation information # Added the force_svs because in some sequences there are slice objects initialised # and recorded but this seems sporadic behaviour. - if ('sSpecPara', 'sVoI', 'sNormal', 'dSag') in mapVBVDHdr['MeasYaps']: - NormaldSag = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dSag')] - elif ('sSliceArray', 'asSlice', '0', 'sNormal', 'dSag') in mapVBVDHdr['MeasYaps']\ + if ('sSliceArray', 'asSlice', '0', 'sNormal', 'dSag') in mapVBVDHdr['MeasYaps']\ and not force_svs: # This is for slice-selective spectroscopy NormaldSag = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sNormal', 'dSag')] + elif ('sSpecPara', 'sVoI', 'sNormal', 'dSag') in mapVBVDHdr['MeasYaps']: + NormaldSag = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dSag')] else: NormaldSag = 0.0 - if ('sSpecPara', 'sVoI', 'sNormal', 'dCor') in mapVBVDHdr['MeasYaps']: - NormaldCor = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dCor')] - elif ('sSliceArray', 'asSlice', '0', 'sNormal', 'dCor') in mapVBVDHdr['MeasYaps']\ + if ('sSliceArray', 'asSlice', '0', 'sNormal', 'dCor') in mapVBVDHdr['MeasYaps']\ and not force_svs: NormaldCor = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sNormal', 'dCor')] + elif ('sSpecPara', 'sVoI', 'sNormal', 'dCor') in mapVBVDHdr['MeasYaps']: + NormaldCor = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dCor')] else: NormaldCor = 0.0 - if ('sSpecPara', 'sVoI', 'sNormal', 'dTra') in mapVBVDHdr['MeasYaps']: - NormaldTra = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dTra')] - elif ('sSliceArray', 'asSlice', '0', 'sNormal', 'dTra') in mapVBVDHdr['MeasYaps']\ + if ('sSliceArray', 'asSlice', '0', 'sNormal', 'dTra') in mapVBVDHdr['MeasYaps']\ and not force_svs: NormaldTra = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sNormal', 'dTra')] + elif ('sSpecPara', 'sVoI', 'sNormal', 'dTra') in mapVBVDHdr['MeasYaps']: + NormaldTra = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sNormal', 'dTra')] else: NormaldTra = 0.0 - if ('sSpecPara', 'sVoI', 'dInPlaneRot') in mapVBVDHdr['MeasYaps']: + if ('sSliceArray', 'asSlice', '0', 'dInPlaneRot') in mapVBVDHdr['MeasYaps']\ + and not force_svs: + inplaneRotation = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'dInPlaneRot')] + elif ('sSpecPara', 'sVoI', 'dInPlaneRot') in mapVBVDHdr['MeasYaps']: inplaneRotation = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'dInPlaneRot')] else: inplaneRotation = 0.0 @@ -606,11 +687,6 @@ def twix2DCMOrientation(mapVBVDHdr, force_svs=False, verbose=False): RoFoV = 10000.0 PeFoV = 10000.0 - dColVec_vector, dRowVec_vector = GSL.calc_prs(TwixSliceNormal, inplaneRotation, verbose) - - imageOrientationPatient = np.stack((dRowVec_vector, dColVec_vector), axis=0) - - pixelSpacing = np.array([PeFoV, RoFoV]) # [RoFoV PeFoV]; if ('sSliceArray', 'asSlice', '0', 'dThickness') in mapVBVDHdr['MeasYaps']\ and not force_svs: sliceThickness = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'dThickness')] @@ -620,17 +696,26 @@ def twix2DCMOrientation(mapVBVDHdr, force_svs=False, verbose=False): sliceThickness = 10000.0 # Position info (including table position) - if ('sSpecPara', 'sVoI', 'sPosition', 'dSag') in mapVBVDHdr['MeasYaps']: + if ('sSliceArray', 'asSlice', '0', 'sPosition', 'dSag') in mapVBVDHdr['MeasYaps']\ + and not force_svs: + PosdSag = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sPosition', 'dSag')] + elif ('sSpecPara', 'sVoI', 'sPosition', 'dSag') in mapVBVDHdr['MeasYaps']: PosdSag = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sPosition', 'dSag')] else: PosdSag = 0.0 - if ('sSpecPara', 'sVoI', 'sPosition', 'dCor') in mapVBVDHdr['MeasYaps']: + if ('sSliceArray', 'asSlice', '0', 'sPosition', 'dCor') in mapVBVDHdr['MeasYaps']\ + and not force_svs: + PosdCor = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sPosition', 'dCor')] + elif ('sSpecPara', 'sVoI', 'sPosition', 'dCor') in mapVBVDHdr['MeasYaps']: PosdCor = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sPosition', 'dCor')] else: PosdCor = 0.0 - if ('sSpecPara', 'sVoI', 'sPosition', 'dTra') in mapVBVDHdr['MeasYaps']: + if ('sSliceArray', 'asSlice', '0', 'sPosition', 'dTra') in mapVBVDHdr['MeasYaps']\ + and not force_svs: + PosdTra = mapVBVDHdr['MeasYaps'][('sSliceArray', 'asSlice', '0', 'sPosition', 'dTra')] + elif ('sSpecPara', 'sVoI', 'sPosition', 'dTra') in mapVBVDHdr['MeasYaps']: PosdTra = mapVBVDHdr['MeasYaps'][('sSpecPara', 'sVoI', 'sPosition', 'dTra')] else: PosdTra = 0.0 @@ -642,15 +727,93 @@ def twix2DCMOrientation(mapVBVDHdr, force_svs=False, verbose=False): if ('lScanRegionPosTra',) in mapVBVDHdr['MeasYaps']: PosdTra += mapVBVDHdr['MeasYaps'][('lScanRegionPosTra',)] - basePosition = np.array([PosdSag, PosdCor, PosdTra], dtype=float) - imagePositionPatient = basePosition + if id_sequence_type(mapVBVDHdr) == 'mrsi': + data_size_pe = 1 + if mapVBVDHdr.Meas.lFinalMatrixSizePhase: + data_size_pe = mapVBVDHdr.Meas.lFinalMatrixSizePhase + + data_size_ro = 1 + if mapVBVDHdr.Meas.lFinalMatrixSizeRead: + data_size_ro = mapVBVDHdr.Meas.lFinalMatrixSizeRead + + data_size_sl = 1 + if mapVBVDHdr.Meas.lFinalMatrixSizeSlice: + data_size_sl = mapVBVDHdr.Meas.lFinalMatrixSizeSlice + + base_pos_info = np.array([PosdSag, PosdCor, PosdTra], dtype=float) + + imageOrientationPatient, imagePositionPatient, pixelSpacing, sliceThickness, dim_swapped = CSIOrientations( + TwixSliceNormal, + inplaneRotation, + PeFoV, + RoFoV, + sliceThickness, + data_size_pe, + data_size_ro, + data_size_sl, + base_pos_info, + verbose) + + else: + dColVec_vector, dRowVec_vector = GSL.calc_prs(TwixSliceNormal, inplaneRotation, verbose) + imageOrientationPatient = np.stack((dRowVec_vector, dColVec_vector), axis=0) + + pixelSpacing = np.array([PeFoV, RoFoV]) # [RoFoV PeFoV]; + + imagePositionPatient = np.array([PosdSag, PosdCor, PosdTra], dtype=float) + + dim_swapped = False + if verbose: print(f'imagePositionPatient is {imagePositionPatient.ravel()}') print(f'imageOrientationPatient is \n{imageOrientationPatient}') print(f'{imageOrientationPatient.ravel()}') print(f'pixelSpacing is {pixelSpacing}') - return imageOrientationPatient, imagePositionPatient, pixelSpacing, sliceThickness + return imageOrientationPatient, imagePositionPatient, pixelSpacing, sliceThickness, dim_swapped + + +def CSIOrientations(slice_normal, ip_rot, fov_pe, fov_ro, fov_sl, n_pe, n_ro, n_sl, base_pos, verbose=False): + SAGITTAL = 0 + CORONAL = 1 + TRANSVERSE = 2 + mo_case = GSL.class_ori(slice_normal[SAGITTAL], slice_normal[CORONAL], slice_normal[TRANSVERSE], False) + dColVec_vector, dRowVec_vector = GSL.calc_prs(slice_normal, ip_rot, verbose) + print('To achieve matching orientation to equivalent DICOM output:') + + if n_sl > 1: + # 3D MRSI + base_pos -= slice_normal * (fov_sl / 2 - fov_sl / n_sl / 2) + fov_sl /= n_sl + + # Sagital + if mo_case == 0: + print('Mirror along ROW/readout direction: VB = LIN, VE+ = SEG') + dRowVec_vector *= -1.0 + pixelSpacing = np.array([fov_ro / n_ro, fov_pe / n_pe]) + dColVec_vector, dRowVec_vector = dRowVec_vector, dColVec_vector + imagePositionPatient =\ + base_pos - (dRowVec_vector * fov_pe / 2) - (dColVec_vector * fov_ro / 2) + dim_swapped = False + # Coronal + elif mo_case == 1: + pixelSpacing = np.array([fov_ro / n_ro, fov_pe / n_pe]) + dColVec_vector, dRowVec_vector = dRowVec_vector, dColVec_vector + imagePositionPatient =\ + base_pos - (dRowVec_vector * fov_pe / 2) - (dColVec_vector * fov_ro / 2) + dim_swapped = False + # Transverse + elif mo_case == 2: + print('Mirror along ROW/readout direction: VB = LIN, VE+ = SEG') + dRowVec_vector *= -1.0 + pixelSpacing = np.array([fov_pe / n_pe, fov_ro / n_ro]) + print('Swap 1st and 2nd diemensions. For VB swap LIN and PHS, for VE+ swap SEG and LIN dimensions') + dim_swapped = True + imagePositionPatient =\ + base_pos - (dRowVec_vector * fov_ro / 2) - (dColVec_vector * fov_pe / 2) + imageOrientationPatient = np.stack((dRowVec_vector, dColVec_vector), axis=0) + + return imageOrientationPatient, imagePositionPatient, pixelSpacing, fov_sl, dim_swapped def examineTwix(twixObj, fileName, mraid): diff --git a/tests/test_ge_pfile_orientation.py b/tests/test_ge_pfile_orientation.py index ed38cf1..f2b5167 100644 --- a/tests/test_ge_pfile_orientation.py +++ b/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) diff --git a/tests/test_philips_dicom_orientation.py b/tests/test_philips_dicom_orientation.py index c0524f0..9c769cb 100644 --- a/tests/test_philips_dicom_orientation.py +++ b/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 diff --git a/tests/test_philips_sdat_spar_orientation.py b/tests/test_philips_sdat_spar_orientation.py index 7824e3a..34caa85 100644 --- a/tests/test_philips_sdat_spar_orientation.py +++ b/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 diff --git a/tests/test_siemens_mrsi_orientation.py b/tests/test_siemens_mrsi_orientation.py index fc2db13..eb80911 100644 --- a/tests/test_siemens_mrsi_orientation.py +++ b/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 diff --git a/tests/test_siemens_svs_orientation.py b/tests/test_siemens_svs_orientation.py index 462d4ec..e2ad66f 100644 --- a/tests/test_siemens_svs_orientation.py +++ b/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')) diff --git a/tests/test_twix.py b/tests/test_twix.py index 12332b7..57b7743 100644 --- a/tests/test_twix.py +++ b/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) diff --git a/tests/test_uih_orientation.py b/tests/test_uih_orientation.py index 63beb2f..a77352a 100644 --- a/tests/test_uih_orientation.py +++ b/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'])