Commit 75fc2ef4 authored by William Clarke's avatar William Clarke
Browse files

Merge branch 'master' into 'master'

Master

See merge request !14
parents 69b1442e 086ab398
Pipeline #6972 passed with stages
in 5 minutes and 8 seconds
###########################################################################
# This file defines the build process for fsl-mrs, as hosted at:
#
# https://git.fmrib.ox.ac.uk/fsl/fsl_mrs
#
# The build pipeline currently comprises the following stages:
#
# 1. style: Check coding style - allowed to fail
#
# 2. test: Unit tests
#
# 3. doc: Building user documentation which appears at:
# https://open.win.ox.ac.uk/pages/fsl/fsl_mrs/
#
# 4. build
# & deploy: Build in three stages fsl-ci-pre, fsl-ci-build,
# sl-ci-deploy
#
# A custom docker image is used for the test job - images are
# available at:
#
# https://hub.docker.com/u/wtclarke/
#
# Stage run conditions:
# Style is run in all cases, but allowed to fail.
# Test is run in all cases.
# Doc is only run on master branches.
# Build stages are run according to the rules associated
# with https://git.fmrib.ox.ac.uk/fsl/fsl-ci-rules
#
###########################################################################
include:
- project: fsl/fsl-ci-rules
file: .gitlab-ci.yml
stages:
- Static Analysis
- Test
- style
- test
- doc
- fsl-ci-pre
- fsl-ci-build
- fsl-ci-deploy
####################################
# These anchors are used to restrict
# when and where jobs are executed.
####################################
.only_upstream: &only_upstream
only:
- branches@fsl/fsl_mrs
.only_master: &only_master
only:
- master
.only_releases: &only_releases
only:
- tags@fsl/fsl_mrs
.except_releases: &except_releases
except:
- tags
# ############
# # 1. style
# ############
flake8:
image: python:3.7-slim-buster
stage: Static Analysis
stage: style
before_script:
- python --version
- pip install flake8
script:
- flake8 --max-line-length=120 fsl_mrs
- flake8 fsl_mrs
allow_failure: true
############
# 2. test
############
pytest:
image: wtclarke/fsl_mrs_tests:1.0
stage: Test
stage: test
variables:
GIT_SUBMODULE_STRATEGY: normal
before_script:
......@@ -36,14 +100,18 @@ pytest:
script:
- pytest fsl_mrs/tests
############
# 3. doc
############
pages:
<<: *only_master
image: python:3.7
stage: doc
script:
- pip install -U sphinx sphinx_rtd_theme
- pip install .
- sphinx-build -b html ./docs/user_docs public
artifacts:
paths:
- public
rules:
- if: '$CI_COMMIT_BRANCH == "master"'
This document contains the FSL-MRS release history in reverse chronological order.
1.0.6 (Tuesday 12th January 2021)
---------------------------------
- Internal changes to core MRS class.
- New plotting functions added, utility functions for plotting added to MRS class.
- fsl_mrs/aux folder renamed for Windows compatibility.
- Moved online documentation to open.win.ox.ac.uk/pages/fsl/fsl_mrs/.
- Fixed small bugs in preprocessing display.
- Synthetic spectra now use fitting model directly.
- Bug fixes in the fsl_Mrs commandline interface. Thanks to Alex Craig-Craven.
- WIP: Dynamic fitting model and dynamic experiment simulation.
- spec2nii requirement pinned to 0.2.11 during NIFTI-MRS development.
1.0.5 (Friday 9th October 2020)
-------------------------------
- Extended documentation of hardcoded constants, including MCMC priors.
......
......@@ -9,7 +9,7 @@ FSL-MRS is a collection of python modules and wrapper scripts for pre-processing
### Installation
#### Conda package
The primary installation method is via _conda_. After installing conda and creating or activating a suitable [enviroment](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html) you can install FSL-MRS from the FSL conda channel.
The primary installation method is via _conda_. After installing conda and creating or activating a suitable [environment](https://docs.conda.io/projects/conda/en/latest/user-guide/tasks/manage-environments.html) you can install FSL-MRS from the FSL conda channel.
conda install -c conda-forge \
-c https://fsl.fmrib.ox.ac.uk/fsldownloads/fslconda/channel/ \
......@@ -32,7 +32,7 @@ or
pip install spec2nii
After installation see the [quick start guide](https://users.fmrib.ox.ac.uk/~saad/fsl_mrs/html/quick_start.html).
After installation see the [quick start guide](https://open.win.ox.ac.uk/pages/fsl/fsl_mrs/quick_start.html).
---
......@@ -59,7 +59,7 @@ After installation see the [quick start guide](https://users.fmrib.ox.ac.uk/~saa
### Documentation
Documentation can be found online [here](https://users.fmrib.ox.ac.uk/~saad/fsl_mrs/html/index.html).
Documentation can be found online on the [WIN open science website](https://open.win.ox.ac.uk/pages/fsl/fsl_mrs/).
For each of the wrapper scripts above, simply type `<name_of_script> --help` to get the usage.
......
import datetime
import fsl_mrs
date = datetime.date.today()
# Configuration file for the Sphinx documentation builder.
......@@ -25,7 +26,7 @@ copyright = f'{date.year}, Will Clarke & Saad Jbabdi, University of Oxford, Oxfo
author = 'William Clarke'
# The full version, including alpha/beta/rc tags
version = '1.0.5'
version = fsl_mrs.__version__
release = version
# From PM's fsleyes doc
......
......@@ -30,7 +30,7 @@ By default this option will add the following basis spectra (in separate metabol
Additional peaks may be added int he interactive environment by calling :code:`add_MM_peaks` with optional arguments to override the defaults.
References
==========
~~~~~~~~~~
.. [CUDA12] Cudalbu C, Mlynárik V, Gruetter R. Handling Macromolecule Signals in the Quantification of the Neurochemical Profile. Journal of Alzheimer’s Disease 2012;31:S101–S115 doi: 10.3233/JAD-2012-120100.
......
......@@ -5,134 +5,104 @@
# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
# Will Clarke <william.clarke@ndcn.ox.ac.uk>
#
# Copyright (C) 2019 University of Oxford
# Copyright (C) 2019 University of Oxford
# SHBASECOPYRIGHT
import warnings
from os.path import isfile
from fsl_mrs.utils import mrs_io as io
from fsl_mrs.utils import misc
from fsl_mrs.utils.constants import *
from fsl_mrs.utils.constants import GYRO_MAG_RATIO, PPM_SHIFT, PPM_RANGE
import numpy as np
#------------------------------------------------
#
#
#------------------------------------------------
class MRS(object):
"""
MRS Class - container for FID, Basis, and sequence info
"""
def __init__(self,FID=None,header=None,basis=None,names=None,basis_hdr=None,H2O=None,cf=None,bw=None):
# If FID and basis are files then read data from file
#if FID is not None and basis is not None:
# if isfile(FID) and isfile(basis):
# self.from_files(FID,basis)
# return
def __init__(self, FID=None, header=None, basis=None, names=None,
basis_hdr=None, H2O=None, cf=None, bw=None, nucleus=None):
# Read in class data input
# (now copying the data - looks ugly but better than referencing.
# now I can run multiple times with different setups)
# If there is no FID data then return empty class object.
# Data is copied.
if FID is not None:
self.set_FID(FID)
else:
return
if H2O is not None:
self.H2O = H2O.copy()
else:
self.H2O = None
self.set_H2O(H2O)
# Set FID class attributes
if header is not None:
self.set_acquisition_params(centralFrequency=header['centralFrequency'],bandwidth=header['bandwidth'])
self.set_acquisition_params(
centralFrequency=header['centralFrequency'],
bandwidth=header['bandwidth'])
self.set_nucleus(header=header, nucleus=nucleus)
self.calculate_axes()
elif (cf is not None) and (bw is not None):
self.set_acquisition_params(centralFrequency=cf,bandwidth=bw)
self.set_acquisition_params(
centralFrequency=cf,
bandwidth=bw)
self.set_nucleus(nucleus=nucleus)
self.calculate_axes()
else:
raise ValueError('You must pass a header or bandwidth and central frequency.')
raise ValueError('You must pass a header'
' or bandwidth, nucleus, and central frequency.')
# Set Basis info
if basis is not None:
self.basis = basis.copy()
# Handle single basis spectra
if self.basis.ndim==1:
self.basis = self.basis[:,np.newaxis]
# Assume that there will always be more timepoints than basis spectra.
if self.basis.shape[0] < self.basis.shape[1]:
self.basis = self.basis.T
self.numBasis = self.basis.shape[1]
self.numBasisPoints = self.basis.shape[0]
if (names is not None) and (basis_hdr is not None):
self.names = names.copy()
self.set_acquisition_params_basis(1/basis_hdr['bandwidth'])
else:
raise ValueError('Pass basis names and header with basis.')
self.set_basis(basis, names, basis_hdr)
# Now interpolate the basis to the same time axis.
self.resample_basis()
# Other properties that need defaults
self.metab_groups = None
self.scaling = {'FID': 1.0, 'basis': 1.0}
else:
self.basis = None
self.names = None
self.numBasis = None
self.basis_dwellTime = None
self.basis_bandwidth = None
# Other properties
self.metab_groups = None
self.scaling = {'FID':1.0,'basis':1.0}
def from_files(self,FID_file,Basis_file,H2O_file=None):
FID,FIDheader = io.read_FID(FID_file)
basis,names,Bheader = io.read_basis(Basis_file)
def from_files(self, FID_file, Basis_file, H2O_file=None):
'''Load data from files into empty MRS class object'''
FID, FIDheader = io.read_FID(FID_file)
basis, names, Bheader = io.read_basis(Basis_file)
if H2O_file is not None:
H2O,_ = io.read_FID(H2O_file)
else:
H2O, _ = io.read_FID(H2O_file)
else:
H2O = None
cf = FIDheader['centralFrequency']
bw = FIDheader['bandwidth']
MRSArgs = {'header': FIDheader,
'basis': basis, 'basis_hdr': Bheader[0],
'names': names}
MRSArgs = {'bw':bw,'cf':cf,
'basis':basis,'basis_hdr':Bheader[0],
'names':names}
self.__init__(FID=FID, H2O=H2O, **MRSArgs)
self.__init__(FID=FID,H2O=H2O,**MRSArgs)
def __str__(self):
cf_MHz = self.centralFrequency / 1e6
cf_T = self.centralFrequency / self.gyromagnetic_ratio / 1e6
return
out = '------- MRS Object ---------\n'
out += f' FID.shape = {self.FID.shape}\n'
out += f' FID.centralFreq (MHz) = {cf_MHz:0.3f}\n'
out += f' FID.nucleus = {self.nucleus}\n'
out += f' FID.centralFreq (T) = {cf_T:0.3f}\n'
out += f' FID.bandwidth (Hz) = {self.bandwidth:0.1f}\n'
out += f' FID.dwelltime (s) = {self.dwellTime:0.5e}\n'
def __str__(self):
out = '------- MRS Object ---------\n'
out += ' FID.shape = {}\n'.format(self.FID.shape)
out += ' FID.centralFreq (MHz) = {}\n'.format(self.centralFrequency/1e6)
out += ' FID.centralFreq (T) = {}\n'.format(self.centralFrequency/H1_gamma/1e6)
out += ' FID.bandwidth (Hz) = {}\n'.format(self.bandwidth)
out += ' FID.dwelltime (s) = {}\n'.format(self.dwellTime)
if self.basis is not None:
out += ' basis.shape = {}\n'.format(self.basis.shape)
out += ' Metabolites = {}\n'.format(self.names)
out += ' numBasis = {}\n'.format(self.numBasis)
out += ' timeAxis = {}\n'.format(self.timeAxis.shape)
out += ' freqAxis = {}\n'.format(self.frequencyAxis.shape)
out += f' basis.shape = {self.basis.shape}\n'
out += f' Metabolites = {self.names}\n'
out += f' numBasis = {self.numBasis}\n'
out += f' timeAxis = {self.timeAxis.shape}\n'
out += f' freqAxis = {self.frequencyAxis.shape}\n'
return out
def __repr__(self) -> str:
return str(self)
# Acquisition parameters
def set_acquisition_params(self,centralFrequency,bandwidth):
def set_acquisition_params(self, centralFrequency, bandwidth):
"""
Set useful params for fitting
......@@ -140,66 +110,122 @@ class MRS(object):
----------
centralFrequency : float (unit=Hz)
bandwidth : float (unit=Hz)
echotime : float (unit=sec)
"""
# Store CF in Hz
self.centralFrequency = misc.checkCFUnits(centralFrequency)
self.centralFrequency = misc.checkCFUnits(centralFrequency)
self.bandwidth = bandwidth
self.dwellTime = 1 / self.bandwidth
self.bandwidth = bandwidth
self.dwellTime = 1/self.bandwidth
def set_acquisition_params_basis(self, dwelltime):
"""
sets basis-specific timing params
"""
# Basis has different dwelltime
self.basis_dwellTime = dwelltime
self.basis_bandwidth = 1 / dwelltime
axes = misc.calculateAxes(self.basis_bandwidth,
self.centralFrequency,
self.numBasisPoints,
self.default_ppm_shift)
self.basis_frequencyAxis = axes['freq']
self.basis_timeAxis = axes['time']
def set_nucleus(self, header=None, nucleus=None):
if nucleus is not None:
self.nucleus = nucleus
elif header is not None and 'ResonantNucleus' in header:
# Look for nucleus string in header
self.nucleus = header['ResonantNucleus']
elif header is not None \
and 'json' in header \
and 'ResonantNucleus' in header['json']:
# Look for nucleus string in header
self.nucleus = header['json']['ResonantNucleus']
else:
# Else try to infer from central frequency range
self.nucleus = self.infer_nucleus(self.centralFrequency)
# Set associated parameters
self.gyromagnetic_ratio = GYRO_MAG_RATIO[self.nucleus]
self.default_ppm_shift = PPM_SHIFT[self.nucleus]
self.default_ppm_range = PPM_RANGE[self.nucleus]
@staticmethod
def infer_nucleus(cf):
cf_MHz = cf / 1e6
for key in GYRO_MAG_RATIO:
onefivet_range = GYRO_MAG_RATIO[key] * np.asarray([1.445, 1.505])
siemensthreet_range = GYRO_MAG_RATIO[key] * np.asarray([2.890, 2.895])
threet_range = GYRO_MAG_RATIO[key] * np.asarray([2.995, 3.005])
sevent_range = GYRO_MAG_RATIO[key] * np.asarray([6.975, 7.005])
ninefourt_range = GYRO_MAG_RATIO[key] * np.asarray([9.35, 9.45])
elevensevent_range = GYRO_MAG_RATIO[key] * np.asarray([11.74, 11.8])
if (cf_MHz > onefivet_range[0] and cf_MHz < onefivet_range[1]) or \
(cf_MHz > siemensthreet_range[0] and cf_MHz < siemensthreet_range[1]) or \
(cf_MHz > threet_range[0] and cf_MHz < threet_range[1]) or \
(cf_MHz > sevent_range[0] and cf_MHz < sevent_range[1]) or \
(cf_MHz > ninefourt_range[0] and cf_MHz < ninefourt_range[1]) or \
(cf_MHz > elevensevent_range[0] and cf_MHz < elevensevent_range[1]):
#print(f'Identified as {key} nucleus data.'
# f' Esitmated field: {cf_MHz/GYRO_MAG_RATIO[key]} T.')
return key
raise ValueError(f'Unidentified nucleus,'
f' central frequency is {cf_MHz} MHz.'
'Pass nucleus parameter explicitly.')
def calculate_axes(self):
''' Calculate axes'''
axes = misc.calculateAxes(self.bandwidth,
self.centralFrequency,
self.numPoints)
self.numPoints,
self.default_ppm_shift)
self.timeAxis = axes['time']
self.frequencyAxis = axes['freq']
self.ppmAxis = axes['ppm']
self.ppmAxisShift = axes['ppmshift']
self.ppmAxisFlip = np.flipud(self.ppmAxisShift)
# turn into column vectors
self.timeAxis = self.timeAxis[:,None]
self.frequencyAxis = self.frequencyAxis[:,None]
self.ppmAxisShift = self.ppmAxisShift[:,None]
self.timeAxis = axes['time']
self.frequencyAxis = axes['freq']
self.ppmAxis = axes['ppm']
self.ppmAxisShift = axes['ppmshift']
# turn into column vectors
self.timeAxis = self.timeAxis[:, None]
self.frequencyAxis = self.frequencyAxis[:, None]
self.ppmAxisShift = self.ppmAxisShift[:, None]
def set_acquisition_params_basis(self,dwelltime):
"""
sets basis-specific timing params
"""
# Basis has different dwelltime
self.basis_dwellTime = dwelltime
self.basis_bandwidth = 1/dwelltime
self.basis_frequencyAxis = np.linspace(-self.basis_bandwidth/2,
self.basis_bandwidth/2,
self.numBasisPoints)
self.basis_timeAxis = np.linspace(self.basis_dwellTime,
self.basis_dwellTime*self.numBasisPoints,
self.numBasisPoints)
def getSpectrum(self,ppmlim=None,shift=True):
def get_spec(self, ppmlim=None, shift=True):
'''Returns spectrum over defined ppm limits
Parameters:
-----------
ppmlim : tuple
shift: bool
'''
spectrum = misc.FIDToSpec(self.FID)
f,l = self.ppmlim_to_range(ppmlim,shift=shift)
return spectrum[f:l]
def getAxes(self,axis='ppmshift',ppmlim=None):
first, last = self.ppmlim_to_range(ppmlim, shift=shift)
return spectrum[first:last]
def getAxes(self, axis='ppmshift', ppmlim=None):
''' Return x axis over defined limits
Axis must be one of ppmshift, ppm, freq,
or time.
'''
if axis.lower() == 'ppmshift':
f,l = self.ppmlim_to_range(ppmlim,shift=True)
return np.squeeze(self.ppmAxisShift[f:l])
first, last = self.ppmlim_to_range(ppmlim, shift=True)
return np.squeeze(self.ppmAxisShift[first:last])
elif axis.lower() == 'ppm':
f,l = self.ppmlim_to_range(ppmlim,shift=False)
return np.squeeze(self.ppmAxis[f:l])
first, last = self.ppmlim_to_range(ppmlim, shift=False)
return np.squeeze(self.ppmAxis[first:last])
elif axis.lower() == 'freq':
f,l = self.ppmlim_to_range(ppmlim,shift=False)
return np.squeeze(self.frequencyAxis[f:l])
first, last = self.ppmlim_to_range(ppmlim, shift=False)
return np.squeeze(self.frequencyAxis[first:last])
elif axis.lower() == 'time':
return np.squeeze(self.timeAxis)
else:
raise ValueError('axis must be one of ppmshift, ppm, freq or time.')
raise ValueError('axis must be one of ppmshift, '
'ppm, freq or time.')
def ppmlim_to_range(self,ppmlim=None,shift=True):
def ppmlim_to_range(self, ppmlim=None, shift=True):
"""
turns ppmlim into data range
......@@ -215,76 +241,79 @@ class MRS(object):
int : last position
"""
if ppmlim is not None:
if shift:
ppm2range = lambda x: np.argmin(np.abs(self.ppmAxisShift-x))
else:
ppm2range = lambda x: np.argmin(np.abs(self.ppmAxis-x))
first = ppm2range(ppmlim[0])
last = ppm2range(ppmlim[1])
if first>last:
first,last = last,first
def ppm2range(x, shift):
if shift:
return np.argmin(np.abs(self.ppmAxisShift - x))
else:
return np.argmin(np.abs(self.ppmAxis - x))
first = ppm2range(ppmlim[0], shift)
last = ppm2range(ppmlim[1], shift)
if first > last:
first, last = last, first
else:
first,last = 0,self.numPoints
return int(first),int(last)
first, last = 0, self.numPoints
return int(first), int(last)
def resample_basis(self):
"""
Usually the basis is simulated using different timings/number of points
Usually the basis is simulated using different
timings and/or number of points.
This interpolates the basis to match the FID
"""
# RESAMPLE BASIS FUNCTION
# bdt = self.basis_dwellTime
# bbw = self.basis_bandwidth
# bn = self.numBasisPoints
# bt = np.linspace(bdt,bdt*bn,bn)-bdt
# fidt = self.timeAxis.flatten()-self.dwellTime
# f = interp1d(bt,self.basis,axis=0)
# newiFB = f(fidt)
self.basis = misc.ts_to_ts(self.basis,self.basis_dwellTime,self.dwellTime,self.numPoints)
self.basis = misc.ts_to_ts(self.basis,
self.basis_dwellTime,
self.dwellTime,
self.numPoints)
self.basis_dwellTime = self.dwellTime
self.basis_bandwidth = 1/self.dwellTime
self.basis_bandwidth = 1 / self.dwellTime
self.numBasisPoints = self.numPoints
# Helper functions