Commit 34462d12 authored by William Clarke's avatar William Clarke
Browse files

Merge branch 'merge_split' into 'master'

mrs_tools. Merge/split NIfTI-MRS handling.

See merge request wclarke/fsl_mrs!4
parents 7d7a1bce dea06f07
This document contains the FSL-MRS release history in reverse chronological order.
1.1.3 (TBC)
------------------------------
- Added mrs_tools script. Replaces mrs_vis and mrs_info. Adds split/merge/reorder functionality.
1.1.2 (Friday 16th April 2021)
------------------------------
- Added 2H information
......
......@@ -51,10 +51,8 @@ After installation see the [quick start guide](https://open.win.ox.ac.uk/pages/f
: Pre-packaged processing for non-edited SVS.
- **fsl\_mrs\_sim**
: simulate basis spectra
- **mrs_vis**
: quick visualisation of the spectra or basis spectra
- **mrs_info**
: quick information on a NIfTI-MRS file
- **mrs_tools**
: Collection of tools for NIfTI-MRS. Includes quick visualisation and information.
- **svs_segment & mrsi_segment**
: Run tissue segmentation for SVS/MRSI from T1 image.
......
Manipulating MRS Data
=====================
Handling NIfTI-MRS
------------------
MRS data stored in NIfTI-MRS format can contain multiple higher dimensions. For example it might contain dimensions encoding multiple receive coils, multiple temporal averages, or even a spectral editing dimension.
Data might need to be manipulated within the NIfTI-MRS storage framework before, after, or during preprocessing. For this, FLS-MRS provides the :code:`mrs_tools` command line script. :code:`mrs_tools` has the ability to merge and split NIfTI-MRS files along the higher encoding dimensions. It can also reorder the higher dimensions, or create a new singleton dimension for further manipulation.
:code:`mrs_tools` also contains the :code:`mrs_tools vis` and :code:`mrs_tools info` options to provide quick visualisation and information on the command line. See the :ref:`Visualisation <visualisation>` page for more information on :code:`mrs_tools vis/info`.
:code:`mrs_tools split` takes a single file and splits it along a specified dimension e.g. :code:`--dim DIM_DYN`, at a single point (:code:`--index 8`) or extracting multiple elements into a second file (:code:`--indices 8 9 10`).
:code:`mrs_tools merge` takes two or more files and merges them along a specified dimension e.g. :code:`--dim DIM_EDIT`.
:code:`mrs_tools reorder` permutes the dimensions of an existing NIfTI-MRS file. For example, the 5th through 7th dimensions can be changed from :code:`DIM_COIL, DIM_DYN, DIM_EDIT` to :code:`DIM_DYN, DIM_EDIT, DIM_COIL` using :code:`--dim_order DIM_DYN DIM_EDIT DIM_COIL`.
\ No newline at end of file
......@@ -19,6 +19,7 @@ If this is your first experience with FSL-MRS, get started with the :ref:`Quick
install
quick_start
data_conversion
data_handling
processing
fitting
quantitation
......
......@@ -34,9 +34,9 @@ But note that there are frequently multiple calibration scans for e.g. shimming
1.1 Take a look at your data
^^^^^^^^^^^^^^^^^^^^^^^^^^^^
You can use :code:`mrs_info` and :code:`mrs_vis` on the command line to view your data at any stage of the process. First mrs_info to see the dimensionality of the data::
You can use :code:`mrs_tools info` and :code:`mrs_tools vis` on the command line to view your data at any stage of the process. First :code:`mrs_tools info` to see the dimensionality of the data::
mrs_info my_metab_file.nii.gz
mrs_tools info my_metab_file.nii.gz
Read file my_metab_file.nii.gz (/path_to_file).
NIfTI-MRS version 0.2
......@@ -47,25 +47,25 @@ You can use :code:`mrs_info` and :code:`mrs_vis` on the command line to view you
Nucleus: 1H
Field Strength: 6.98 T
Then mrs_vis to visualise the data::
Then :code:`mrs_tools vis` to visualise the data::
mrs_vis my_metab_file.nii.gz
mrs_tools vis my_metab_file.nii.gz
:code:`mrs_vis` will automatically perform coil combination and averaging down to a single spectrum for display purposes only.
:code:`mrs_tools vis` will automatically perform coil combination and averaging down to a single spectrum for display purposes only.
.. image:: data/raw_conv.png
:width: 600
You can also quickly view data across one of the NIfTI-MRS higher dimensions (those containing uncombined coils, or averages etc.) In this case we plot all the transients stored in the dimension tagged *DIM_DYN* (i.e. averages)::
mrs_vis my_metab_file.nii.gz --display_dim DIM_DYN
mrs_tools vis my_metab_file.nii.gz --display_dim DIM_DYN
.. image:: data/mrs_vis_dir.png
:width: 600
If you see a significantly different picture (no data, just noise, etc.) stop and investigate. See :ref:`Troubleshooting <TS_4>`.
Have a look at the :ref:`Visualisation <visualisation>` page for more information on :code:`mrs_vis`.
Have a look at the :ref:`Visualisation <visualisation>` page for more information on :code:`mrs_tools vis`.
2. Process your raw data
~~~~~~~~~~~~~~~~~~~~~~~~
......@@ -101,9 +101,9 @@ The fitting in FSL-MRS requires the user to provide basis spectra. Basis spectra
`my_sequence_description.json` contains a description of the sequence broken down into blocks of RF pulses and gradients. This must be created for each sequence manually once. `metabs.txt` contains a list of metabolites to simulate. Much more information on constructing a suitable sequence description JSON file can be found on the :ref:`Basis Spectra Simulation <simulation>` page.
Have a quick check of your basis set using mrs_vis::
Have a quick check of your basis set using :code:`mrs_tools vis`::
mrs_vis my_basis_spectra/
mrs_tools vis my_basis_spectra/
4. Tissue Segmentation
~~~~~~~~~~~~~~~~~~~~~~
......
......@@ -31,7 +31,7 @@ Troubleshooting hints
.. _TS_4:
4. Data looks 'wrong' after conversion
If when using :code:`mrs_vis` you see no signal and just noise try conjugating the data using :code:`fsl_mrs_proc conj` or try expanding the ppm range plotted :code:`--ppmlim -10 10`. If you see a flat line, then conversion failed. The data might be corrupted - did the acquisition complete successfully?
If when using :code:`mrs_tools vis` you see no signal and just noise try conjugating the data using :code:`fsl_mrs_proc conj` or try expanding the ppm range plotted :code:`--ppmlim -10 10`. If you see a flat line, then conversion failed. The data might be corrupted - did the acquisition complete successfully?
.. image:: data/bad_data.png
:width: 600
\ No newline at end of file
......@@ -15,32 +15,32 @@ There are 4 ways of visualising/interacting with MRS data in FSL-MRS:
1. Quick glance
---------------
The first thing one might want to do when given a FID file or simulated spectra is to have a quick look at the spectra to see if they look like one would expect. To get a sense of the dimensionality and basic status of the data run :code:`mrs_info` for a quick text summary. FSL-MRS then provides a light-weight script (:code:`mrs_vis`) to quickly visualise the MRS or MRSI data. For example, running :code:`mrs_vis` on the provided example SVS data:
The first thing one might want to do when given a FID file or simulated spectra is to have a quick look at the spectra to see if they look like one would expect. To get a sense of the dimensionality and basic status of the data run :code:`mrs_tools info` for a quick text summary. FSL-MRS then provides a light-weight script (:code:`mrs_tools vis`) to quickly visualise the MRS or MRSI data. For example, running :code:`mrs_tools vis` on the provided example SVS data:
::
mrs_vis example_usage/example_data/metab.nii
mrs_tools vis example_usage/example_data/metab.nii
gives the following basic plot:
.. image:: data/mrs_vis_svs.png
:width: 400
Note that the reason :code:`mrs_vis` "knows" how to scale the x-axis is that the relevant information is stored in the NIfTI-MRS MRS header extension (namely the *dwell time* and the *central frequency*).
Note that the reason :code:`mrs_tools vis` "knows" how to scale the x-axis is that the relevant information is stored in the NIfTI-MRS MRS header extension (namely the *dwell time* and the *central frequency*).
:code:`mrs_vis` can also visualise a folder of mrs data::
:code:`mrs_tools vis` can also visualise a folder of mrs data::
mrs_vis ./converted_data_dir/
mrs_tools vis ./converted_data_dir/
.. image:: data/mrs_vis_dir.png
:width: 600
We can also run :code:`mrs_vis` to visualise the metabolite basis. Again below we do so for the simulated basis provided with the example data:
We can also run :code:`mrs_tools vis` to visualise the metabolite basis. Again below we do so for the simulated basis provided with the example data:
::
mrs_vis example_usage/example_data/steam_11ms/
mrs_tools vis example_usage/example_data/steam_11ms/
.. image:: data/mrs_vis_basis.png
......
......@@ -138,7 +138,7 @@
"outputs": [],
"source": [
"%%capture\n",
"%sx mrs_vis --save basis_fig.png steam_basis"
"%sx mrs_tools vis --save basis_fig.png steam_basis"
]
},
{
......@@ -170,4 +170,4 @@
},
"nbformat": 4,
"nbformat_minor": 4
}
}
\ No newline at end of file
%% Cell type:markdown id: tags:
# Example basis spectra creation
This notebook demos the process of creating basis spectra for fitting in FSL-MRS.
### Contents:
- [1. Sequence JSON file](#1.-Sequence-JSON-file)
- [1.1. Defining MM](#1.1-Defining-MM)
- [2. Simulation](#2.-Simulation)
- [3. Visualisation](#3.-Visualisation)
Will Clarke
June 2020
University of Oxford
%% Cell type:markdown id: tags:
## 1. Sequence JSON file
%% Cell type:markdown id: tags:
The sequence to simulate and the system parameters are descibed in a json file. The file breaks the sequence into a series of RF pulses (+ slice selcect gradients) and delays with optional rephasing gradients. A coherence filter is used to crush unwanted coherences.
%% Cell type:markdown id: tags:
## 1.1 Defining MM
An experimentally measured (or otherwise derived) macromollecues basis spectrum can be included in the basis spectrum by defining an additional JSON file.
%% Cell type:markdown id: tags:
## 2. Simulation
%% Cell type:code id: tags:
``` python
%sx fsl_mrs_sim -b example_data/metabs.txt \
-o steam_basis \
-p -2.66 \
--MM example_data/macSeed.json \
--overwrite \
-e 11.0 \
example_data/steam11.json
```
%%%% Output: execute_result
["Identified spinsystems: ['Ala', 'Asc', 'Asp', 'GPC', 'PCh', 'Cr', 'PCr', 'GABA', 'Glc', 'Gln', 'Glu', 'GSH', 'Ins', 'Lac', 'NAA', 'NAAG', 'PE', 'Scyllo', 'Tau']",
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Auto-phase adjustment. Phasing peak position = 1.99 ppm',
'Rx_Phase: 2.3408',
'Additional phase: 0.059',
'Final Rx_Phase: 2.400',
'Running simulation on Asc.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Glu.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on NAA.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on PCr.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on GSH.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on GABA.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Lac.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Scyllo.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Tau.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Cr.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Ins.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Ala.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Gln.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on NAAG.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Asp.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on Glc.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on PCh.',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Simulation running using mode 1d. Axis order = [0 1 2].',
'Running simulation on PE.',
'Simulation running using mode 1d. Axis order = [0 1 2].']
%% Cell type:markdown id: tags:
## 3. Visualisation
%% Cell type:code id: tags:
``` python
%%capture
%sx mrs_vis --save basis_fig.png steam_basis
%sx mrs_tools vis --save basis_fig.png steam_basis
```
%% Cell type:markdown id: tags:
![basis spectra](basis_fig.png)
......
......@@ -8,6 +8,7 @@
import json
from pathlib import Path
import re
import numpy as np
from nibabel.nifti1 import Nifti1Extension
......@@ -106,12 +107,11 @@ class NIFTI_MRS(Image):
std_tags = ['DIM_COIL', 'DIM_DYN', 'DIM_INDIRECT_0']
for idx in range(3):
curr_dim = idx + 5
if self.ndim >= curr_dim:
curr_tag = f'dim_{curr_dim}'
if curr_tag in self.hdr_ext:
self.dim_tags[idx] = self.hdr_ext[curr_tag]
else:
self.dim_tags[idx] = std_tags[idx]
curr_tag = f'dim_{curr_dim}'
if curr_tag in self.hdr_ext:
self.dim_tags[idx] = self.hdr_ext[curr_tag]
elif curr_dim < self.ndim:
self.dim_tags[idx] = std_tags[idx]
def __getitem__(self, sliceobj):
'''Apply conjugation at use. This swaps from the
......@@ -171,6 +171,7 @@ class NIFTI_MRS(Image):
extension = Nifti1Extension(44, json_s.encode('UTF-8'))
self.header.extensions.clear()
self.header.extensions.append(extension)
self._set_dim_tags()
def dim_position(self, dim_tag):
'''Return position of dim if it exists.'''
......@@ -187,6 +188,83 @@ class NIFTI_MRS(Image):
dim += 4
return dim
def add_hdr_field(self, key, value):
"""Add a field to the header extension
To do: validate type (standard or user)
:param key: Field key
:type key: str
:param value: Value of field to add
"""
dim_n = re.compile(r'dim_[567].*')
if dim_n.match(key):
raise ValueError('Modify dimension headers through dedicated methods.')
current_hdr_ext = self.hdr_ext
current_hdr_ext.update({key: value})
self.hdr_ext = current_hdr_ext
def remove_hdr_field(self, key):
"""Remove a field from the header extension
:param key: Key to remove
:type key: str
"""
if key == 'SpectrometerFrequency' or key == 'ResonantNucleus':
raise ValueError('You cannot remove the required metadata.')
dim_n = re.compile(r'dim_[567].*')
if dim_n.match(key):
raise ValueError('Modify dimension headers through dedicated methods.')
current_hdr_ext = self.hdr_ext
current_hdr_ext.pop(key, None)
self.hdr_ext = current_hdr_ext
def set_dim_info(self, dim, info_str):
"""Set or update the 'dim_N_info' field
:param dim: The dim tag or python dimension index (i.e. N-1)
:type dim: str or int
:param info_str: New info string
:type info_str: str
"""
dim = self._dim_tag_to_index(dim)
current_hdr_ext = self.hdr_ext
current_hdr_ext[f'dim_{dim + 1}_info'] = info_str
self.hdr_ext = current_hdr_ext
def set_dim_header(self, dim, hdr_obj):
"""Set or update the 'dim_N_header' field
hdr_obj replaces the current value.
:param dim: The dim tag or python dimension index (i.e. N-1)
:type dim: str or int
:param hdr_obj: dict containing the dimension headers
:type hdr_obj: dict
"""
dim = self._dim_tag_to_index(dim)
# Check size
def size_chk(obj):
# Allow for expansion along the next dimension
if dim == self.ndim:
dim_len = 1
else:
dim_len = self.shape[dim]
if len(obj) != dim_len:
raise ValueError(f'New dim header length must be {self.shape[dim]}')
for key in hdr_obj:
if isinstance(hdr_obj[key], list):
size_chk(hdr_obj[key])
elif isinstance(hdr_obj[key], dict)\
and 'Value' in hdr_obj[key]:
size_chk(hdr_obj[key]['Value'])
current_hdr_ext = self.hdr_ext
current_hdr_ext[f'dim_{dim + 1}_header'] = hdr_obj
self.hdr_ext = current_hdr_ext
def copy(self, remove_dim=None):
'''Return a copy of this image, optionally with a dimension removed.
Args:
......@@ -204,10 +282,13 @@ class NIFTI_MRS(Image):
if dd > new_obj.ndim:
hdr_ext.pop(f'dim_{dd}', None)
hdr_ext.pop(f'dim_{dd}_header', None)
hdr_ext.pop(f'dim_{dd}_info', None)
elif dd >= dim:
hdr_ext[f'dim_{dd}'] = hdr_ext[f'dim_{dd + 1}']
if f'dim_{dd + 1}_header' in hdr_ext:
hdr_ext[f'dim_{dd}_header'] = hdr_ext[f'dim_{dd + 1}_header']
if f'dim_{dd + 1}_info' in hdr_ext:
hdr_ext[f'dim_{dd}_info'] = hdr_ext[f'dim_{dd + 1}_info']
new_obj.hdr_ext = hdr_ext
new_obj._set_dim_tags()
......
#!/usr/bin/env python
# mrs_info - quick NIfTI MRS information
#
# Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
#
# Copyright (C) 2021 University of Oxford
# SHBASECOPYRIGHT
import argparse
from pathlib import Path
def main():
p = argparse.ArgumentParser(description='FSL MRS Tools: NIfTI-MRS information')
p.add_argument('file', type=Path, metavar='FILE or list of FILEs',
help='NIfTI MRS file(s)', nargs='+')
args = p.parse_args()
from fsl_mrs.utils.mrs_io import read_FID
from fsl_mrs.utils.constants import GYRO_MAG_RATIO
for file in args.file:
data = read_FID(str(file))
print(f'\nRead file {file.name} ({file.parent.resolve()}).')
print(f'NIfTI-MRS version {data.mrs_nifti_version}')
print(f'Data shape {data.shape}')
print(f'Dimension tags: {data.dim_tags}')
print(f'Spectrometer Frequency: {data.spectrometer_frequency[0]} MHz')
print(f'Dwelltime (Bandwidth): {data.dwelltime:0.3E}s ({data.bandwidth:0.0f} Hz)')
print(f'Nucleus: {data.nucleus[0]}')
if data.nucleus[0] in GYRO_MAG_RATIO:
field_strength = data.spectrometer_frequency[0] / GYRO_MAG_RATIO[data.nucleus[0]]
print(f'Field Strength: {field_strength:0.2f} T')
print()
if __name__ == '__main__':
main()
#!/usr/bin/env python
""" mrs_tools - top level script for calling general mrs handling tools
Author: William Clarke <william.clarke@ndcn.ox.ac.uk>
Copyright (C) 2021 University of Oxford
SHBASECOPYRIGHT
"""
import argparse
from pathlib import Path
from fsl_mrs import __version__
def main():
# Parse command-line arguments
p = argparse.ArgumentParser(description="FSL Magnetic Resonance Spectroscopy - Tools")
p.add_argument('-v', '--version', action='version', version=__version__)
sp = p.add_subparsers(title='subcommands',
description='Availible tools',
required=True,
dest='subcommand')
# Info tool
infoparser = sp.add_parser(
'info',
help='Information about the NIfTI-MRS file.')
infoparser.add_argument(
'file',
type=Path,
metavar='FILE or list of FILEs',
help='NIfTI MRS file(s)', nargs='+')
infoparser.set_defaults(func=info)
# Vis tool
visparser = sp.add_parser(
'vis',
help='Quick visualisation of a NIfTI-MRS file or FSL-MRS basis set.')
visparser.add_argument('file', type=Path, metavar='FILE or DIR',
help='NIfTI file or directory of basis sets')
visparser.add_argument('--ppmlim', default=(.2, 4.2), type=float,
nargs=2, metavar=('LOW', 'HIGH'),
help='limit the fit to a freq range (default=(.2,4.2))')
visparser.add_argument('--mask', default=None, type=str, help='Mask for MRSI')
visparser.add_argument('--save', default=None, type=str, help='Save fig to path')
visparser.add_argument('--display_dim', default=None, type=str,
help='NIFTI-MRS tag. Do not average across this dimension.')
visparser.set_defaults(func=vis)
# Merge tool - Merge NIfTI MRS along higher dimensions
mergeparser = sp.add_parser(
'merge',
help='Merge NIfTI-MRS along higher dimensions.')
mergeparser.add_argument('--files', type=Path, required=True, nargs='+',
help='List of files to merge')
mergeparser.add_argument('--dim', type=str, required=True,
help='NIFTI-MRS dimension tag to merge across.')
mergeparser.add_argument('--output',
required=True, type=Path, default=Path('.'),
help='output folder (defaults to current directory)')
mergeparser.add_argument('--filename', type=str,
help='Override output file name.')
mergeparser.set_defaults(func=merge)
# Split tool
splitparser = sp.add_parser(
'split',
help='Split NIfTI-MRS along higher dimensions.')
splitparser.add_argument('--file', type=Path, required=True,
help='File to split')
splitparser.add_argument('--dim', type=str, required=True,
help='NIFTI-MRS dimension tag to split across.')
group = splitparser.add_mutually_exclusive_group(required=True)
group.add_argument('--indices', type=int, nargs='+',
help='List of indices to extract into second file.'
'All indices are zero-indexed.')
group.add_argument('--index', type=int, nargs='+',
help='Index to split at (split after index, zero-indexed).')
splitparser.add_argument('--output',
required=True, type=Path, default=Path('.'),
help='output folder (defaults to current directory)')
splitparser.add_argument('--filename', type=str,
help='Override output file names.')
splitparser.set_defaults(func=split)
# Reorder tool
reorderparser = sp.add_parser(
'reorder',
help='Reorider higher dimensions of NIfTI-MRS.')
reorderparser.add_argument('--file', type=Path, required=True,
help='File to reorder')
reorderparser.add_argument('--dim_order', type=str, nargs='+', required=True,
help='NIFTI-MRS dimension tags in desired order. '
'Enter as strings (min:1, max:3). '
'Can create singleton dimension at end.')
reorderparser.add_argument('--output',
required=True, type=Path, default=Path('.'),
help='output folder (defaults to current directory)')
reorderparser.add_argument('--filename', type=str,
help='Override output file names.')
reorderparser.set_defaults(func=reorder)
# Parse command-line arguments
args = p.parse_args()
# Call function
args.func(args)
def info(args):
"""Prints basic information about NIfTI-MRS files
:param args: Argparse interpreted arguments
:type args: Namespace
"""
from fsl_mrs.utils.mrs_io import read_FID
from fsl_mrs.utils.constants import GYRO_MAG_RATIO
for file in args.file:
data = read_FID(str(file))
print(f'\nRead file {file.name} ({file.parent.resolve()}).')
print(f'NIfTI-MRS version {data.mrs_nifti_version}')
print(f'Data shape {data.shape}')
print(f'Dimension tags: {data.dim_tags}')
print(f'Spectrometer Frequency: {data.spectrometer_frequency[0]} MHz')
print(f'Dwelltime (Bandwidth): {data.dwelltime:0.3E}s ({data.bandwidth:0.0f} Hz)')
print(f'Nucleus: {data.nucleus[0]}')
if data.nucleus[0] in GYRO_MAG_RATIO:
field_strength = data.spectrometer_frequency[0] / GYRO_MAG_RATIO[data.nucleus[0]]
print(f'Field Strength: {field_strength:0.2f} T')
print()
def vis(args):
"""Visualiser for NIfTI-MRS files
:param args: Argparse interpreted arguments
:type args: Namespace
"""
from fsl_mrs.utils.plotting import plot_spectrum, FID2Spec, plot_spectra
from fsl_mrs.utils.mrs_io import read_FID, read_basis
import matplotlib.pyplot as plt
from fsl_mrs.core import MRS
import numpy as np
from fsl_mrs.utils.preproc import nifti_mrs_proc
import nibabel as nib
# Some logic to figure out what we are dealing with
p = args.file
nifti_files = list(p.glob('*.nii*'))
# Identify BASIS
if p.is_dir() and len(nifti_files) == 0 or \
p.suffix.upper() == '.BASIS':
basis, names, basishdr = read_basis(args.file)
fid = np.zeros(basis.shape[0])
mrs = MRS(FID=fid,
header=basishdr[0],
basis=basis,
names=names,
basis_hdr=basishdr[0])
mrs.check_Basis(repair=True)
first, last = mrs.ppmlim_to_range(ppmlim=args.ppmlim)
plt.figure(figsize=(8, 8))
for idx, n in enumerate(names):
plt.plot(mrs.getAxes(ppmlim=args.ppmlim),
np.real(FID2Spec(mrs.basis[:, idx]))[first:last],
label=n)
plt.gca().invert_xaxis()
plt.xlabel('Chemical shift (ppm)')
plt.legend()
if args.save is not None:
plt.savefig(args.save)
else:
plt.show()
# Identify directory of nifti files
elif p.is_dir() and len(nifti_files) > 0:
raise ValueError('mrs_tools vis should be called on a single'
' NIFTI-MRS file, not a directory (unless'
' it contains basis files).')
# Single nifti file
elif p.is_file():
data = read_FID(args.file)
if data.ndim > 4 and 'DIM_COIL' in data.dim_tags:
print('Performing coil combination')
data = nifti_mrs_proc.coilcombine(data)
if np.prod(data.shape[:3]) == 1:
# SVS
if args.display_dim:
for idx in range(data.ndim - 4):
if data.dim_tags[idx] != args.display_dim:
print(f'Averaging {data.dim_tags[idx]}')
data = nifti_mrs_proc.average(data, data.dim_tags[idx])
fig = plot_spectra(data.mrs(), ppmlim=args.ppmlim)
else:
while data.ndim > 4:
print(f'Averaging {data.dim_tags[0]}')
data = nifti_mrs_proc.average(data, data.dim_tags[0])
fig = plot_spectrum(data.mrs(), ppmlim=args.ppmlim)
if args.save is not None:
fig.savefig(args.save)
else:
plt.show()
else:
while data.ndim > 4:
print(f'Averaging {data.dim_tags[0]}')
data = nifti_mrs_proc.average(data, data.dim_tags[0])
mrsi = data.mrs()
if args.mask is not None:
mask_hdr = nib.load(args.mask)
mask = np.asanyarray(mask_hdr.dataobj)
if mask.ndim == 2:
mask = np.expand_dims(mask, 2)
mrsi.set_mask(mask)
mrsi.plot()
def merge(<