Commit 96bc72d4 authored by William Clarke's avatar William Clarke
Browse files

fsl_mrs_preproc modified and test passes.

parent 796eae2c
......@@ -85,8 +85,12 @@ class NIFTI_MRS(Image):
super().__init__(*args, **kwargs)
# Check that file meets minimum requirements
if float(self.mrs_nifti_version) < 0.2:
raise NotNIFTI_MRS('NIFTI-MRS > V0.2 required.')
try:
nmrs_version = self.mrs_nifti_version
if float(nmrs_version) < 0.2:
raise NotNIFTI_MRS('NIFTI-MRS > V0.2 required.')
except IndexError:
raise NotNIFTI_MRS('NIFTI-MRS intent code not set.')
if 44 not in self.header.extensions.get_codes():
raise NotNIFTI_MRS('NIFTI-MRS must have a header extension.')
......
This diff is collapsed.
# Test the loading functions
# Test io functions
import fsl_mrs.utils.mrs_io as mrsio
import fsl_mrs.utils.mrs_io.fsl_io as fslio
from fsl_mrs.utils import plotting
......@@ -6,11 +6,12 @@ import numpy as np
import os.path as op
testsPath = op.dirname(__file__)
SVSTestData = {'nifti':op.join(testsPath,'testdata/mrs_io/metab.nii'),
'raw':op.join(testsPath,'testdata/mrs_io/metab.RAW'),
'txt':op.join(testsPath,'testdata/mrs_io/metab.txt')}
SVSTestData = {'nifti': op.join(testsPath, 'testdata/mrs_io/metab.nii'),
'raw': op.join(testsPath, 'testdata/mrs_io/metab.RAW'),
'txt': op.join(testsPath, 'testdata/mrs_io/metab.txt')}
headerReqFields = ['centralFrequency', 'bandwidth', 'dwelltime']
headerReqFields = ['centralFrequency','bandwidth','dwelltime']
def test_read_FID_SVS():
# Test the loading of the three types of data we handle for SVS data
......@@ -18,35 +19,40 @@ def test_read_FID_SVS():
# .raw
# .txt
data_nifti,header_nifti = mrsio.read_FID(SVSTestData['nifti'])
data_raw,header_raw = mrsio.read_FID(SVSTestData['raw'])
data_txt,header_txt = mrsio.read_FID(SVSTestData['txt'])
data_nifti = mrsio.read_FID(SVSTestData['nifti'])
data_raw = mrsio.read_FID(SVSTestData['raw'])
data_txt = mrsio.read_FID(SVSTestData['txt'])
# Check that the data from each of these matches - it should they are all the same bit of data.
datamean = np.mean([data_nifti,data_raw,data_txt],axis=0)
datamean = np.mean([data_nifti.data,
data_raw.data,
data_txt.data], axis=0)
assert np.isclose(data_nifti,datamean).all()
assert np.isclose(data_raw,datamean).all()
assert np.isclose(data_txt,datamean).all()
assert np.isclose(data_nifti.data, datamean).all()
assert np.isclose(data_raw.data, datamean).all()
assert np.isclose(data_txt.data, datamean).all()
# Check that the headers each contain the required fields
for r in headerReqFields:
assert r in header_nifti
assert r in header_raw
assert r in header_txt
headerMean = np.mean([header_nifti[r],header_raw[r],header_txt[r]])
assert np.isclose(header_nifti[r],headerMean)
assert np.isclose(header_raw[r],headerMean)
assert np.isclose(header_txt[r],headerMean)
# # Check that the headers each contain the required fields
# for r in headerReqFields:
# assert r in header_nifti
# assert r in header_raw
# assert r in header_txt
# headerMean = np.mean([header_nifti[r], header_raw[r], header_txt[r]])
# assert np.isclose(header_nifti[r], headerMean)
# assert np.isclose(header_raw[r], headerMean)
# assert np.isclose(header_txt[r], headerMean)
# TODO: Make MRSI test function (and find data)
# def test_read_FID_MRSI()
BasisTestData = {'fsl':op.join(testsPath,'testdata/mrs_io/basisset_FSL'),
'raw':op.join(testsPath,'testdata/mrs_io/basisset_LCModel_raw'),
'txt':op.join(testsPath,'testdata/mrs_io/basisset_JMRUI'),
'lcm':op.join(testsPath,'testdata/mrs_io/basisset_LCModel.BASIS')}
BasisTestData = {'fsl': op.join(testsPath, 'testdata/mrs_io/basisset_FSL'),
'raw': op.join(testsPath, 'testdata/mrs_io/basisset_LCModel_raw'),
'txt': op.join(testsPath, 'testdata/mrs_io/basisset_JMRUI'),
'lcm': op.join(testsPath, 'testdata/mrs_io/basisset_LCModel.BASIS')}
def test_read_Basis():
# Test the loading of the four types of data we handle for basis specta
# fsl_mrs - folder of json
......@@ -54,13 +60,13 @@ def test_read_Basis():
# lcmodel - folder of .raw
# jmrui - folder of .txt
basis_fsl,names_fsl,headers_fsl = mrsio.read_basis(BasisTestData['fsl'])
basis_raw,names_raw,headers_raw = mrsio.read_basis(BasisTestData['raw'])
basis_txt,names_txt,headers_txt = mrsio.read_basis(BasisTestData['txt'])
basis_lcm,names_lcm,headers_lcm = mrsio.read_basis(BasisTestData['lcm'])
basis_fsl, names_fsl, headers_fsl = mrsio.read_basis(BasisTestData['fsl'])
basis_raw, names_raw, headers_raw = mrsio.read_basis(BasisTestData['raw'])
basis_txt, names_txt, headers_txt = mrsio.read_basis(BasisTestData['txt'])
basis_lcm, names_lcm, headers_lcm = mrsio.read_basis(BasisTestData['lcm'])
#lcm basis file is zeropadded by a factor of 2 remove
basis_lcm = basis_lcm[:2048,:]
# lcm basis file is zeropadded by a factor of 2, remove
basis_lcm = basis_lcm[:2048, :]
# Test that all contain the same amount of data.
expectedDataSize = (2048, 21)
......@@ -68,44 +74,51 @@ def test_read_Basis():
assert basis_raw.shape == expectedDataSize
assert basis_txt.shape == expectedDataSize
assert basis_lcm.shape == expectedDataSize
# Test that the number of names match the amount of data
numNames = 21
assert len(names_fsl)==numNames
assert len(names_raw)==numNames
assert len(names_txt)==numNames
assert len(names_lcm)==numNames
assert len(names_fsl) == numNames
assert len(names_raw) == numNames
assert len(names_txt) == numNames
assert len(names_lcm) == numNames
# Check that the headers each contain the required fields
# Exclude raw, we know it doesn't contain everything
for r in headerReqFields:
assert r in headers_fsl[0]
assert r in headers_fsl[0]
assert r in headers_txt[0]
assert r in headers_lcm[0]
headerMean = np.mean([headers_fsl[0][r],headers_txt[0][r],headers_lcm[0][r]])
headerMean = np.mean([headers_fsl[0][r],
headers_txt[0][r],
headers_lcm[0][r]])
if r == 'centralFrequency':
assert np.isclose(headers_fsl[0][r],headerMean,rtol=2e-01, atol=1e05)
assert np.isclose(headers_txt[0][r],headerMean,rtol=2e-01, atol=1e05)
assert np.isclose(headers_lcm[0][r],headerMean,rtol=2e-01, atol=1e05)
assert np.isclose(headers_fsl[0][r], headerMean, rtol=2e-01, atol=1e05)
assert np.isclose(headers_txt[0][r], headerMean, rtol=2e-01, atol=1e05)
assert np.isclose(headers_lcm[0][r], headerMean, rtol=2e-01, atol=1e05)
else:
assert np.isclose(headers_fsl[0][r],headerMean)
assert np.isclose(headers_txt[0][r],headerMean)
assert np.isclose(headers_lcm[0][r],headerMean)
assert np.isclose(headers_fsl[0][r], headerMean)
assert np.isclose(headers_txt[0][r], headerMean)
assert np.isclose(headers_lcm[0][r], headerMean)
# Conjugate fsl and jMRUI
basis_fsl = basis_fsl.conj()
basis_txt = basis_txt.conj()
# Test that all contain roughly the same data when scaled.
metabToCheck = 'Cr'
checkIdx = names_raw.index('Cr')
normAbsSpec = lambda spec:np.abs(spec)/np.max(np.abs(spec))
convertToLimitedSpec = lambda fid: normAbsSpec(plotting.FID2Spec(fid)[900:1000])
meanSpec = np.mean([convertToLimitedSpec(basis_fsl[:,checkIdx]),
convertToLimitedSpec(basis_raw[:,checkIdx]),
convertToLimitedSpec(basis_txt[:,checkIdx]),
convertToLimitedSpec(basis_lcm[:,checkIdx])],axis=0)
checkIdx = names_raw.index(metabToCheck)
def normAbsSpec(spec):
return np.abs(spec) / np.max(np.abs(spec))
def convertToLimitedSpec(fid):
return normAbsSpec(plotting.FID2Spec(fid)[900:1000])
meanSpec = np.mean([convertToLimitedSpec(basis_fsl[:, checkIdx]),
convertToLimitedSpec(basis_raw[:, checkIdx]),
convertToLimitedSpec(basis_txt[:, checkIdx]),
convertToLimitedSpec(basis_lcm[:, checkIdx])], axis=0)
# breakpoint()
# import matplotlib.pyplot as plt
# plt.plot(convertToLimitedSpec(basis_fsl[:,checkIdx]))
......@@ -113,17 +126,21 @@ def test_read_Basis():
# plt.plot(convertToLimitedSpec(basis_txt[:,checkIdx]),'-.')
# plt.plot(convertToLimitedSpec(basis_lcm[:,checkIdx]),':')
# plt.show()
assert np.allclose(convertToLimitedSpec(basis_fsl[:,checkIdx]),meanSpec,rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_raw[:,checkIdx]),meanSpec,rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_txt[:,checkIdx]),meanSpec,rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_lcm[:,checkIdx]),meanSpec,rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_fsl[:, checkIdx]), meanSpec, rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_raw[:, checkIdx]), meanSpec, rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_txt[:, checkIdx]), meanSpec, rtol=2e-01, atol=1e-03)
assert np.allclose(convertToLimitedSpec(basis_lcm[:, checkIdx]), meanSpec, rtol=2e-01, atol=1e-03)
def test_fslBasisRegen():
pointsToGen = 10
basis_fsl,names_fsl,headers_fsl = mrsio.read_basis(BasisTestData['fsl'])
basis_fsl2,names_fsl2,headers_fsl2 = fslio.readFSLBasisFiles(BasisTestData['fsl'],readoutShift=4.65,bandwidth=4000,points=pointsToGen)
basis_fsl, names_fsl, headers_fsl = mrsio.read_basis(BasisTestData['fsl'])
basis_fsl2, names_fsl2, headers_fsl2 = fslio.readFSLBasisFiles(BasisTestData['fsl'],
readoutShift=4.65,
bandwidth=4000,
points=pointsToGen)
assert np.isclose(basis_fsl2,basis_fsl[:10,:]).all()
assert np.allclose(basis_fsl2, basis_fsl[:10, :])
assert names_fsl2 == names_fsl
for r in headerReqFields:
assert headers_fsl[0][r] == headers_fsl2[0][r]
assert headers_fsl[0][r] == headers_fsl2[0][r]
......@@ -3,27 +3,30 @@ from pathlib import Path
testsPath = Path(__file__).parent
data = testsPath / 'testdata/fsl_mrs_preproc'
t1 = str(testsPath / 'testdata/svs_segment/T1.anat/T1_biascorr.nii.gz')
def test_preproc(tmp_path):
allfiles_metab = list(data.glob('steam_metab_raw*.nii.gz'))
allfiles_wrefc = list(data.glob('steam_wref_comb_raw*.nii.gz'))
allfiles_wrefq = list(data.glob('steam_wref_quant_raw*.nii.gz'))
allfiles_ecc = list(data.glob('steam_ecc_raw*.nii.gz'))
metab = str(data / 'metab_raw.nii.gz')
wrefc = str(data / 'wref_raw.nii.gz')
wrefq = str(data / 'quant_raw.nii.gz')
ecc = str(data / 'ecc.nii.gz')
retcode = subprocess.check_call(
['fsl_mrs_preproc',
'--output', str(tmp_path),
'--data'] + allfiles_metab +
['--reference'] + allfiles_wrefc +
['--quant'] + allfiles_wrefq +
['--ecc'] + allfiles_ecc +
['--hlsvd',
'--leftshift', '1',
'--overwrite',
'--report',
'--verbose'])
['fsl_mrs_preproc',
'--output', str(tmp_path),
'--data', metab,
'--reference', wrefc,
'--quant', wrefq,
'--ecc', ecc,
'--t1', t1,
'--hlsvd',
'--leftshift', '1',
'--overwrite',
'--report',
'--verbose'])
assert retcode == 0
assert (tmp_path / 'mergedReports.html').exists()
assert (tmp_path / 'voxel_location.png').exists()
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
......@@ -11,7 +11,43 @@ import os.path as op
from fsl_mrs.core.NIFTI_MRS import gen_new_nifti_mrs
# Read jMRUI style text files
def readjMRUItxt_fid(txtfile):
'''Read FID files in jMRUI txt format
:param txtfile: .txt format file path
:return: NIFTI_MRS
'''
data, header = readjMRUItxt(txtfile, unpack_header=True)
if 'TypeOfNucleus' in header['jmrui']:
nucleus = header['jmrui']['TypeOfNucleus']
else:
nucleus = '1H'
data = data.reshape((1, 1, 1) + data.shape)
return gen_new_nifti_mrs(data, header['dwelltime'], header['centralFrequency'], nucleus=nucleus)
# Read jMRUI .txt files containing basis
def read_txtBasis_files(txtfiles):
basis = []
names = []
header = []
for file in txtfiles:
b, h = readjMRUItxt(file)
basis.append(b)
prefix, _ = op.splitext(op.basename(file))
names.append(prefix)
header.append(h)
basis = np.array(basis).conj().T
return basis, names, header
# generically read jMRUI style text files
def readjMRUItxt(filename, unpack_header=True):
"""
Read .txt format file
......@@ -57,12 +93,7 @@ def readjMRUItxt(filename, unpack_header=True):
# Clean up header
header = translateHeader(header)
if 'TypeOfNucleus' in header['jmrui']:
nucleus = header['jmrui']['TypeOfNucleus']
else:
nucleus = '1H'
return gen_new_nifti_mrs(data, header['dwelltime'], header['centralFrequency'], nucleus=nucleus)
return data, header
# Translate jMRUI header to mandatory fields
......@@ -84,24 +115,6 @@ def num(s):
return s
# Read jMRUI .txt files containing basis
def read_txtBasis_files(txtfiles):
basis = []
names = []
header = []
for file in txtfiles:
b, h = readjMRUItxt(file)
basis.append(b)
prefix, _ = op.splitext(op.basename(file))
names.append(prefix)
header.append(h)
basis = np.array(basis).conj().T
return basis, names, header
# Write functions
def writejMRUItxt(fileout, FID, paramDict):
......
......@@ -64,6 +64,7 @@ def read_lcm_raw_h2o(filename):
:return: NIFTI_MRS
"""
data, header = readLCModelRaw(filename)
data = data.reshape((1, 1, 1) + data.shape)
return gen_new_nifti_mrs(data, header['dwelltime'], header['centralFrequency'])
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment