From 9b868689e1baa157b1b567b6f3d4b51a7f8fdffb Mon Sep 17 00:00:00 2001 From: spinicist Date: Thu, 25 Jun 2020 08:42:41 +0100 Subject: [PATCH] ENH qipype is born --- Python/QUIT/__init__.py | 0 Python/QUIT/interfaces/relax.py | 828 ------------------ Python/T1_mapping.ipynb | 190 ++-- Python/Tests/test_despot_sc.py | 36 +- Python/Tests/test_mt.py | 59 +- Python/Tests/test_perfusion.py | 31 +- Python/Tests/test_relax.py | 18 +- Python/Tests/test_ssfp.py | 40 +- Python/Tests/test_utils.py | 13 +- Python/example_notebook.ipynb | 2 +- Python/qipype/__init__.py | 4 + Python/{QUIT => qipype}/base.py | 165 ++-- Python/qipype/interfaces/b1.py | 121 +++ Python/{QUIT => qipype}/interfaces/core.py | 2 +- Python/{QUIT => qipype}/interfaces/fsl.py | 0 Python/{QUIT => qipype}/interfaces/mt.py | 232 ++--- .../{QUIT => qipype}/interfaces/perfusion.py | 67 +- Python/qipype/interfaces/relax.py | 449 ++++++++++ Python/{QUIT => qipype}/interfaces/rufis.py | 0 Python/{QUIT => qipype}/interfaces/stats.py | 0 .../interfaces/susceptibility.py | 0 Python/{QUIT => qipype}/interfaces/utils.py | 8 +- Python/{QUIT => qipype}/sims.py | 4 + Python/{QUIT => qipype}/workflows/cest.py | 6 +- Python/{QUIT => qipype}/workflows/mpm.py | 4 +- Python/{QUIT => qipype}/workflows/prep.py | 4 +- Python/setup.py | 18 +- README.md | 62 +- Source/Core/Args.h | 2 +- Source/Core/SimulateModel.h | 2 +- Source/MT/qi_qmt.cpp | 7 +- Source/Perfusion/qi_ase_oef.cpp | 5 - Source/Relaxometry/qidespot2.cpp | 6 +- Source/Relaxometry/qidespot2fm.cpp | 6 +- 34 files changed, 1070 insertions(+), 1321 deletions(-) delete mode 100644 Python/QUIT/__init__.py delete mode 100644 Python/QUIT/interfaces/relax.py create mode 100644 Python/qipype/__init__.py rename Python/{QUIT => qipype}/base.py (59%) create mode 100644 Python/qipype/interfaces/b1.py rename Python/{QUIT => qipype}/interfaces/core.py (98%) rename Python/{QUIT => qipype}/interfaces/fsl.py (100%) rename Python/{QUIT => qipype}/interfaces/mt.py (56%) rename Python/{QUIT => qipype}/interfaces/perfusion.py (73%) create mode 100644 Python/qipype/interfaces/relax.py rename Python/{QUIT => qipype}/interfaces/rufis.py (100%) rename Python/{QUIT => qipype}/interfaces/stats.py (100%) rename Python/{QUIT => qipype}/interfaces/susceptibility.py (100%) rename Python/{QUIT => qipype}/interfaces/utils.py (98%) rename Python/{QUIT => qipype}/sims.py (95%) rename Python/{QUIT => qipype}/workflows/cest.py (98%) rename Python/{QUIT => qipype}/workflows/mpm.py (95%) rename Python/{QUIT => qipype}/workflows/prep.py (98%) diff --git a/Python/QUIT/__init__.py b/Python/QUIT/__init__.py deleted file mode 100644 index e69de29b..00000000 diff --git a/Python/QUIT/interfaces/relax.py b/Python/QUIT/interfaces/relax.py deleted file mode 100644 index a5f9af2c..00000000 --- a/Python/QUIT/interfaces/relax.py +++ /dev/null @@ -1,828 +0,0 @@ -#! /usr/bin/env python -# -*- coding: utf-8 -*- - -""" -Implementation of nipype interfaces for QUIT relaxometry utilities -Requires that the QUIT tools are in your your system path -""" - -from os import path, getcwd -from nipype.interfaces.base import TraitedSpec, File, traits -from .. import base as QI - - -############################ qidespot1 ############################ - -class DESPOT1InputSpec(QI.FitInputSpec): - # Additional Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - algo = traits.String(desc="Choose algorithm (l/w/n)", argstr="--algo=%s") - iterations = traits.Int( - desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') - clamp_PD = traits.Float( - desc='Clamp PD between 0 and value', argstr='-f %f') - clamp_T1 = traits.Float( - desc='Clamp T1 between 0 and value', argstr='--clampT1=%f') - - -class DESPOT1OutputSpec(TraitedSpec): - t1_map = File('D1_T1.nii.gz', desc="Path to T1 map", usedefault=True) - pd_map = File('D1_PD.nii.gz', desc="Path to PD map", usedefault=True) - rmse_map = File('D1_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class DESPOT1(QI.FitCommand): - """ - Run DESPOT1 analysis with qidespot1 - - Example with parameter file - ------- - >>> from QUIT.nipype.relaxometry import QIDespot1 - >>> d1 = QIDespot1(prefix='nipype_', param_file='spgr_params.json') - >>> d1.inputs.in_file = 'SPGR.nii.gz' - >>> d1_res = d1.run() - >>> print(d1_res.outputs) - - Example with parameter dictionary - ------- - >>> from QUIT.nipype.relaxometry import QIDespot1 - >>> params = {'SPGR': {'TR':5E-3, 'FA':[5,10]} } - >>> d1 = QIDespot1(prefix='nipype_', param_dict=params) - >>> d1.inputs.in_file = 'SPGR.nii.gz' - >>> d1_res = d1.run() - >>> print(d1_res.outputs) - - """ - - _cmd = 'qi despot1' - input_spec = DESPOT1InputSpec - output_spec = DESPOT1OutputSpec - -############################ qidespot1sim ############################ - - -class DESPOT1SimInputSpec(QI.SimInputSpec): - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - - -class DESPOT1Sim(QI.SimCommand): - """ - Run DESPOT1 simulation with qidespot1 - - Example with parameter dictionary - ------- - >>> from QUIT.nipype.relaxometry import QIDespot1Sim - >>> params = {'SPGR': {'TR':5E-3, 'FA':[5,10]}, - 'T1File': 'D1_T1.nii.gz', - 'PDFile': 'D1_PD.nii.gz'} - >>> d1sim = QIDespot1Sim(prefix='nipype_', param_dict=params) - >>> d1sim.inputs.spgr_file = 'SPGR.nii.gz' - >>> d1sim_res = d1.run() - >>> print(d1sim_res.outputs) - - """ - - _cmd = 'qi despot1' - _param_files = ['PD', 'T1'] - input_spec = DESPOT1SimInputSpec - output_spec = QI.SimOutputSpec - -############################ qidespot1hifi ############################ - - -class HIFIInputSpec(QI.InputSpec): - # Inputs - spgr_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to SPGR data') - - mprage_file = File(exists=True, argstr='%s', mandatory=True, - position=-1, desc='Path to MPRAGE data') - - # Options - residuals = traits.Bool( - desc='Write out residuals for each data-point', argstr='--resids') - clamp_T1 = traits.Float( - desc='Clamp T1 between 0 and value', argstr='--clampT1=%f') - - -class HIFIOutputSpec(TraitedSpec): - t1_map = File('HIFI_T1.nii.gz', desc="Path to T1 map", usedefault=True) - pd_map = File('HIFI_PD.nii.gz', desc="Path to PD map", usedefault=True) - b1_map = File('HIFI_B1.nii.gz', desc="Path to B1 map", usedefault=True) - rmse_map = File('HIFI_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class HIFI(QI.FitCommand): - """ - Calculate T1 & B1 map with the DESPOT1-HIFI method - - Example - ------- - >>> from QUIT.nipype.relaxometry import QIDespot1Hifi - >>> params = {'SPGR': {'TR':5E-3, 'FA':[5,10]}, - 'MPRAGE': { 'FA': 5, 'TR': 5E-3, 'TI': 0.45, 'TD': 0, 'eta': 1, 'ETL': 64, 'k0': 0 }} - >>> hifi = QIDespot1Hifi(prefix='nipype_', param_dict=params) - >>> hifi.inputs.spgr_file = 'SPGR.nii.gz' - >>> hifi.inputs.mprage_file = 'MPRAGE.nii.gz' - >>> hifi_res = hifi.run() - >>> print(hifi_res.outputs) - - """ - - _cmd = 'qi despot1hifi' - input_spec = HIFIInputSpec - output_spec = HIFIOutputSpec - -############################ qidespot1hifisim ############################ - - -class HIFISimInputSpec(QI.SimInputBaseSpec): - # Inputs - spgr_file = File(exists=False, argstr='%s', mandatory=True, - position=-2, desc='Path for output SPGR image') - - mprage_file = File(exists=False, argstr='%s', mandatory=True, - position=-1, desc='Path for output MPRAGE image') - - -class HIFISimOutputSpec(TraitedSpec): - spgr_file = File(desc="Output SPGR image") - mprage_file = File(desc="Output MPRAGE image") - - -class HIFISim(QI.SimCommand): - """ - Simulate SPGR/FLASH and MPRAGE images using DESPOT1-HIFI model - - Example - ------- - >>> from QUIT.nipype.relaxometry import QIDespot1HifiSim - >>> params = {'SPGR': {'TR':5E-3, 'FA':[5,10]}, - 'MPRAGE': { 'FA': 5, 'TR': 5E-3, 'TI': 0.45, 'TD': 0, 'eta': 1, 'ETL': 64, 'k0': 0 }, - 'PDFile': 'PD.nii.gz', - 'T1File': 'T1.nii.gz'} - >>> hifisim = QIDespot1HifiSim(prefix='nipype_', param_dict=params) - >>> hifisim.inputs.spgr_file = 'SPGR.nii.gz' - >>> hifisim.inputs.mprage_file = 'MPRAGE.nii.gz' - >>> hifisim_res = hifi.run() - >>> print(hifisim_res.outputs) - - """ - - _cmd = 'qi despot1hifi' - _param_files = ['PD', 'T1', 'B1'] - input_spec = HIFISimInputSpec - output_spec = HIFISimOutputSpec - - def _list_outputs(self): - outputs = self.output_spec().get() - outputs['spgr_file'] = self.inputs.spgr_file - outputs['mprage_file'] = self.inputs.mprage_file - return outputs - -############################ qidespot2 ############################ - - -class DESPOT2InputSpec(QI.FitInputSpec): - # Inputs - t1_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to T1 map') - - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - algo = traits.Enum("LLS", "WLS", "NLS", - desc="Choose algorithm", argstr="--algo=%d") - ellipse = traits.Bool( - desc="Data is ellipse geometric solution", argstr='--gs') - iterations = traits.Int( - desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') - clamp_PD = traits.Float( - desc='Clamp PD between 0 and value', argstr='-f %f') - clamp_T2 = traits.Float( - desc='Clamp T2 between 0 and value', argstr='--clampT1=%f') - - -class DESPOT2OutputSpec(TraitedSpec): - t2_map = File('D2_T2.nii.gz', desc="Path to T2 map", usedefault=True) - pd_map = File('D2_PD.nii.gz', desc="Path to PD map", usedefault=True) - rmse_map = File('D2_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class DESPOT2(QI.FitCommand): - """ - Run DESPOT2 analysis with qidespot2 - - Example with parameter file - ------- - >>> from QUIT.nipype.relaxometry import QIDespot2 - >>> d1 = QIDespot2(prefix='nipype_', param_file='ssfp_params.json') - >>> d2.inputs.in_file = 'SSFP.nii.gz' - >>> d2.inputs.t1_file = 'D1_T1.nii.gz' - >>> d2_res = d2.run() - >>> print(d2_res.outputs) - - Example with parameter dictionary - ------- - >>> from QUIT.nipype.relaxometry import QIDespot2 - >>> params = {'SSFP': {'TR':5E-3, 'FA':[10,50]} } - >>> d2 = QIDespot2(prefix='nipype_', param_dict=params) - >>> d2.inputs.in_file = 'SSFP.nii.gz' - >>> d2.inputs.t1_file = 'D1_T1.nii.gz' - >>> d2_res = d2.run() - >>> print(d2_res.outputs) - - """ - - _cmd = 'qi despot2' - input_spec = DESPOT2InputSpec - output_spec = DESPOT2OutputSpec - -############################ qidespot2sim ############################ - - -class DESPOT2SimInputSpec(QI.SimInputSpec): - # Inputs - t1_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path for input T1 map') - - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - ellipse = traits.Bool( - desc="Data is ellipse geometric solution", argstr='--gs') - - -class DESPOT2Sim(QI.SimCommand): - """ - Run DESPOT2 simulation - - Example with parameter dictionary - ------- - >>> from QUIT.nipype.relaxometry import QIDespot2Sim - >>> params = {'SSFP': {'TR':5E-3, 'FA':[10,50]}, - 'PDFile': 'PD.nii.gz', - 'T1File': 'T1.nii.gz' } - >>> d2sim = QIDespot2Sim(prefix='nipype_', param_dict=params) - >>> d2sim.inputs.in_file = 'SSFP.nii.gz' - >>> d2sim.inputs.t1_file = 'D1_T1.nii.gz' - >>> d2sim_res = d2.run() - >>> print(d2sim_res.outputs) - - """ - - _cmd = 'qi despot2' - _param_files = ['PD', 'T2'] - input_spec = DESPOT2SimInputSpec - output_spec = QI.SimOutputSpec - -############################ qidespot2fm ############################ - - -class FMInputSpec(QI.FitInputSpec): - # Inputs - t1_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to T1 map') - - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - asym = traits.Bool( - desc="Fit asymmetric (+/-) off-resonance frequency", argstr='--asym') - algo = traits.Enum("LLS", "WLS", "NLS", - desc="Choose algorithm", argstr="--algo=%d") - - -class FMOutputSpec(TraitedSpec): - t2_map = File('FM_T2.nii.gz', desc="Path to T2 map", usedefault=True) - pd_map = File('FM_PD.nii.gz', desc="Path to PD map", usedefault=True) - f0_map = File('FM_f0.nii.gz', - desc="Path to f0 (off-resonance) map", usedefault=True) - rmse_map = File('FM_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class FM(QI.FitCommand): - """ - Run DESPOT2-FM analysis - - Example - ------- - >>> from QUIT.nipype.relaxometry import QIDespot2FM - >>> params = {'SSFP': {'TR':5E-3, 'FA':[10,10,50,50], 'PhaseInc':[180,180,0,0] } - >>> fm = QIDespot2FM(prefix='nipype_', param_dict=params) - >>> fm.inputs.in_file = 'SSFP.nii.gz' - >>> fm.inputs.t1_file = 'D1_T1.nii.gz' - >>> fm_res = fm.run() - >>> print(fm_res.outputs) - - """ - - _cmd = 'qi despot2fm' - input_spec = FMInputSpec - output_spec = FMOutputSpec - - -class FMSimInputSpec(QI.SimInputSpec): - # Inputs - t1_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path for input T1 map') - - # Options - b1_map = File(desc='B1 map (ratio)', argstr='--B1=%s') - asym = traits.Bool( - desc="Fit asymmetric (+/-) off-resonance frequency", argstr='--asym') - - -class FMSim(QI.SimCommand): - """ - Run DESPOT2-FM simulation - - """ - - _cmd = 'qi despot2fm' - _param_files = ['PD', 'T2', 'f0'] - input_spec = FMSimInputSpec - output_spec = QI.SimOutputSpec - -############################ qi_jsr ############################ - - -class JSRInputSpec(QI.InputSpec): - # Inputs - spgr_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to SPGR data') - - ssfp_file = File(exists=True, argstr='%s', mandatory=True, - position=-1, desc='Path to SSFP data') - - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - residuals = traits.Bool( - desc='Write out residuals for each data-point', argstr='--resids') - npsi = traits.Int( - desc='Number of psi/off-resonance starts', argstr='--npsi=%d') - - -class JSROutputSpec(TraitedSpec): - pd_map = File('JSR_PD.nii.gz', desc='Path to PD map', usedefault=True) - t1_map = File('JSR_T1.nii.gz', desc='Path to T1 map', usedefault=True) - t2_map = File('JSR_T2.nii.gz', desc='Path to T2 map', usedefault=True) - f0_map = File('JSR_df0.nii.gz', - desc='Path to off-resonance map', usedefault=True) - rmse_map = File('JSR_rmse.nii.gz', - desc='Path to residual map', usedefault=True) - - -class JSR(QI.FitCommand): - """ - Calculate T1 &T2 map with Joint System Relaxometry - - """ - - _cmd = 'qi jsr' - input_spec = JSRInputSpec - output_spec = JSROutputSpec - -############################ qimcdespot ############################ -# Status: Everything is there but not tested - - -class QIMCDespotInputSpec(QI.InputSpec): - # Inputs - spgr_file = File(exists=True, argstr='%s', mandatory=True, - position=0, desc='SPGR file') - ssfp_file = File(exists=True, argstr='%s', mandatory=True, - position=1, desc='SSFP file') - param_dict = traits.Dict(desc='Parameters as dictionary', position=2, - argstr='', mandatory=True, xor=['param_file']) - param_file = File(desc='Parameters as .json file', position=2, argstr='--json=%s', - xor=['param_dict'], mandatory=True, exists=True) - - # Options - f0_map = File(desc='f0 map (Hertz)', argstr='--f0=%s', exists=True) - B1_map = File(desc='B1 map (ratio)', argstr='--B1=%s', exists=True) - mask = File(desc='Only process voxels within the mask', - argstr='--mask=%s', exists=True) - residuals = traits.Bool( - desc='Write out residuals for each data-point', argstr='--resid') - model = traits.Enum("1", "2", "2nex", "3", "3_f0", "3nex", desc="Select model to fit - 1/2/2nex/3/3_f0/3nex, default 3", - argstr="--model=%d") - scale = traits.Bool( - desc='Normalize signals to mean (a good idea)', argstr='--scale') - algo = traits.Enum( - "S", "G", desc="Select (S)tochastic or (G)aussian Region Contraction", argstr="--algo=%d") - iterations = traits.Int( - desc='Max iterations, default 4', argstr='---its=%d') - field_strength = traits.Float( - desc='Specify field-strength for fitting regions - 3/7/u for user input', argstr='--tesla=%f') - threads = traits.Int( - desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') - prefix = traits.String( - desc='Add a prefix to output filenames', argstr='--out=%s') - - -class QIMCDespotOutputSpec(TraitedSpec): - # Specify which outputs there are - output_3C_T1_m = File(desc="T1 of myelin water") - output_3C_T2_m = File(desc="T2 of myelin water") - output_3C_T1_ie = File(desc="T1 of intra/extra-cellular water") - output_3C_T2_ie = File(desc="T2 of intra/extra-cellular water") - output_3C_T1_csf = File(desc="T1 of CSF") - output_3C_T2_csf = File(desc="T2 of CSF") - output_3C_tau_m = File( - desc="The residence time of myelin water (reciprocal of forward exchange rate)") - output_3C_f_m = File(desc="The Myelin Water Fraction (MWF)") - output_3C_f_csf = File(desc="The CSF Fraction") - output_3C_f0 = File( - desc="The off-resonance frequency. If this was specified on the command line, it will be a copy of that file") - output_3C_B1 = File( - desc="The relative flip-angle map. If this was specified on the command line, it will be a copy of that file") - - -class QIMCDespot(QI.FitCommand): - """ - Interace for qimcdespot - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QiMcDespot - >>> interface = QiMcDespot(prefix='nipype_', param_file='mcdespot_params.json') - """ - - _cmd = 'qi mcdespot' - input_spec = QIMCDespotInputSpec - output_spec = QIMCDespotOutputSpec - - # If the command requires a json input file - def _format_arg(self, name, spec, value): - return self._process_params(name, spec, value) - - def _list_outputs(self): - outputs = self.output_spec().get() - - # Specify output files - outputs = ['3C_T1_m', '3C_T2_m', '3C_T1_ie', '3C_T2_ie', '3C_T1_csf', '3C_T2_csf', - '3C_tau_m', '3C_f_m', '3C_f_csf', '3C_f0', '3C_B1'] - for op in outputs: - outputs['output_{}'.format(op)] = os.path.abspath( - self._add_prefix('{}.nii.gz'.format(op))) - return outputs - -############################ qimp2rage ############################ -# Implemented but not tested # - - -class MP2RAGEInputSpec(QI.InputSpec): - # Inputs - in_file = File(exists=True, argstr='%s', mandatory=True, - position=0, desc='Path to complex MP-RAGE data') - - # Commonly used options - threads = traits.Int( - desc='Use N threads (default=hardware)', argstr='--threads=%d') - prefix = traits.String( - desc='Add a prefix to output filenames', argstr='--out=%s') - beta = traits.Float(desc='Regularisation paramter', argstr='--beta=%f') - - -class MP2RAGEOutputSpec(TraitedSpec): - # Specify which outputs there are - uni_file = File('MP2_UNI.nii.gz', - desc='The Uniform MP2 contrast image', usedefault=True) - t1_map = File('MP2_T1.nii.gz', desc='T1 Map', usedefault=True) - - -class MP2RAGE(QI.FitCommand): - """ - Interface for qimp2rage - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QIMP2RAGE - >>> interface = QIMP2RAGE(prefix='nipype_', param_file='spgr_params.json') - - """ - - _cmd = 'qi mp2rage' - input_spec = MP2RAGEInputSpec - output_spec = MP2RAGEOutputSpec - -############################ qimultiecho ############################ - - -class MultiechoInputSpec(QI.FitInputSpec): - # Options - algo = traits.String(desc="Choose algorithm (l/a/n)", argstr="--algo=%s") - iterations = traits.Int( - desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') - thresh_PD = traits.Float( - desc='Only output maps when PD exceeds threshold value', argstr='-t=%f') - clamp_T2 = traits.Float( - desc='Clamp T2 between 0 and value', argstr='-p=%f') - - -class MultiechoOutputSpec(TraitedSpec): - t2_map = File( - 'ME_T2.nii.gz', desc='The T2 map. Units are the same as TE1 and ESP', usedefault=True) - pd_map = File( - 'ME_PD.nii.gz', desc='The apparent proton-density map (intercept of the decay curve at TE=0)', usedefault=True) - rmse_map = File('ME_rmse.nii.gz', - desc='Path to residual map', usedefault=True) - - -class Multiecho(QI.FitCommand): - """ - help for myInterface - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QiMultiecho - >>> interface = QiMultiecho(prefix='nipype_', param_file='me_params.json') - - """ - - _cmd = 'qi multiecho' - input_spec = MultiechoInputSpec - output_spec = MultiechoOutputSpec - - -class MultiechoSim(QI.SimCommand): - """ - help for myInterface - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QiMultiecho - >>> interface = QiMultiecho(prefix='nipype_', param_file='me_params.json') - - """ - - _cmd = 'qi multiecho' - _param_files = ['PD', 'T2'] - input_spec = QI.SimInputSpec - output_spec = QI.SimOutputSpec - -############################ qi_mpm_r2s ############################ - - -class MPMR2sInputSpec(QI.InputSpec): - # Inputs - pdw_file = File(exists=True, argstr='%s', mandatory=True, - position=-3, desc='Path to PD-weighted data') - - t1w_file = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to T1-weighted data') - - mtw_file = File(exists=True, argstr='%s', mandatory=True, - position=-1, desc='Path to MT-weighted data') - - # Options - residuals = traits.Bool( - desc='Write out residuals for each data-point', argstr='--resids') - - -class MPMR2sOutputSpec(TraitedSpec): - r2s_map = File( - 'MPM_R2s.nii.gz', desc='The R2* map. Units are the same as TE1 and ESP', usedefault=True) - s0_pdw = File('MPM_S0_PDw.nii.gz', - desc='Intercept of the decay curve at TE=0 for PDw', usedefault=True) - s0_t1w = File('MPM_S0_T1w.nii.gz', - desc='Intercept of the decay curve at TE=0 for T1w', usedefault=True) - s0_mtw = File('MPM_S0_MTw.nii.gz', - desc='Intercept of the decay curve at TE=0 for MTw', usedefault=True) - - -class MPMR2s(QI.FitCommand): - """ - Runs qi_mpm_r2s - - """ - - _cmd = 'qi mpm_r2s' - input_spec = MPMR2sInputSpec - output_spec = MPMR2sOutputSpec - -############################ qidream ############################ -# Implemented but not tested # - - -class QIDreamInputSpec(QI.InputSpec): - # Inputs - dream_file = File(exists=True, argstr='%s', mandatory=True, - position=0, desc='Input file. Must have 2 volumes (FID and STE)') - - # Options - threads = traits.Int( - desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') - prefix = traits.String( - desc='Add a prefix to output filenames', argstr='--out=%s') - order = traits.String( - desc='Volume order - f/s/v - fid/ste/vst first', argstr='--order=%s') - mask_file = File( - desc='Only process voxels within the mask', argstr='--mask=%s') - alpha = traits.Float( - desc="Nominal flip-angle (default 55)", argstr="--alpha=%f") - - -class QIDreamOutputSpec(TraitedSpec): - b1_rel_map = File(desc="The relative flip-angle map.") - b1_act_map = File(desc="The actual achieved angle in each voxel.") - - -class QIDream(QI.BaseCommand): - """ - Interface for qidream - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QiDream - >>> interface = QiDream(prefix='nipype_', param_file='spgr_params.json') - - """ - - _cmd = 'qi dream' - input_spec = QIDreamInputSpec - output_spec = QIDreamOutputSpec - - # If the command requires a json input file - def _format_arg(self, name, spec, value): - return self._process_params(name, spec, value) - - def _list_outputs(self): - outputs = self.output_spec().get() - outputs['b1_rel_map'] = os.path.abspath( - self._add_prefix('DREAM_B1.nii.gz')) - outputs['b1_act_map'] = os.path.abspath( - self._add_prefix('DREAM_angle.nii.gz')) - - return outputs - -############################ qiafi ############################ -# Implemented but not tested # - - -class QIAFIInputSpec(QI.InputSpec): - # Inputs - afi_file = File(exists=True, argstr='%s', mandatory=True, - position=0, desc='Input file') - - # Options - threads = traits.Int( - desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') - prefix = traits.String( - desc='Add a prefix to output filenames', argstr='--out=%s') - alpha = traits.Float( - desc="Specify nominal flip-angle, default 55", argstr="--flip=%f") - tr_ratio = traits.Float( - desc="Specify TR2:TR1 ratio, default 5", argstr="--ratio=%f") - save_act_b1 = traits.Bool( - desc="Write out the actual flip-angle as well as B1", argstr="--save") - - -class QIAFIOutputSpec(TraitedSpec): - # Specify which outputs there are - b1_rel_map = File(desc="The relative flip-angle map.") - b1_act_map = File(desc="The actual flip-angle map.") - - -class QIAFI(QI.BaseCommand): - """ - help for myInterface - - Example 1 - ------- - >>> from QUIT.nipype.relaxometry import QiAfi - >>> interface = QiAfi(prefix='nipype_', param_file='spgr_params.json') - - """ - - _cmd = 'qi afi' - input_spec = QIAFIInputSpec - output_spec = QIAFIOutputSpec - - # If the command requires a json input file - def _format_arg(self, name, spec, value): - return self._process_params(name, spec, value) - - def _list_outputs(self): - outputs = self.output_spec().get() - outputs['b1_rel_map'] = os.path.abspath( - self._add_prefix('AFI_B1.nii.gz')) - outputs['b1_act_map'] = os.path.abspath( - self._add_prefix('AFI_angle.nii.gz')) - return outputs - -############################ qi_ssfp_elipse ############################ - - -class EllipseInputSpec(QI.FitInputSpec): - # Additional Options - algo = traits.String(desc='Choose algorithm (h/d)', argstr='--algo=%s') - - -class EllipseOutputSpec(TraitedSpec): - G_map = File('ES_G.nii.gz', desc='Path to G map', usedefault=True) - a_map = File('ES_a.nii.gz', desc='Path to a map', usedefault=True) - b_map = File('ES_b.nii.gz', desc='Path to b map', usedefault=True) - theta0_map = File('ES_theta_0.nii.gz', - desc='Path to theta 0 (off-resonance phase) map', usedefault=True) - phi_rf_map = File('ES_phi_rf.nii.gz', - desc='Path to RF phase map', usedefault=True) - rmse_map = File('ES_rmse.nii.gz', - desc='Path to residual map', usedefault=True) - - -class Ellipse(QI.FitCommand): - """ - Fit an ellipse to SSFP data - - """ - - _cmd = 'qi ssfp_ellipse' - input_spec = EllipseInputSpec - output_spec = EllipseOutputSpec - - -class EllipseSim(QI.SimCommand): - """ - Simulate SSFP data from ellipse parameters - - """ - - _cmd = 'qi ssfp_ellipse' - _param_files = ['G', 'a', 'b', 'theta_0', 'phi_rf'] - input_spec = QI.SimInputSpec - output_spec = QI.SimOutputSpec - -############################ qi_planet ############################ - - -class PLANETInputSpec(QI.InputSpec): - # Inputs - G_map = File(exists=True, argstr='%s', mandatory=True, - position=-3, desc='Path to G parameter map') - a_map = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to a parameter map') - b_map = File(exists=True, argstr='%s', mandatory=True, - position=-1, desc='Path to b parameter map') - - # Options - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - residuals = traits.Bool( - desc='Write out residuals for each data-point', argstr='--resids') - - -class PLANETOutputSpec(TraitedSpec): - PD_map = File('PLANET_PD.nii.gz', desc="Path to PD map", usedefault=True) - T1_map = File('PLANET_T1.nii.gz', desc="Path to T1 map", usedefault=True) - T2_map = File('PLANET_T2.nii.gz', desc="Path to T2 map", usedefault=True) - rmse_map = File('PLANET_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class PLANET(QI.FitCommand): - """ - Calculate T1/T2 from ellipse parameters - - """ - - _cmd = 'qi planet' - input_spec = PLANETInputSpec - output_spec = PLANETOutputSpec - - -class PLANETSimInputSpec(QI.SimInputBaseSpec): - G_file = File(argstr='%s', mandatory=True, - position=-3, desc='Output G file') - a_file = File(argstr='%s', mandatory=True, - position=-2, desc='Output a file') - b_file = File(argstr='%s', mandatory=True, - position=-1, desc='Output b file') - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - - -class PLANETSimOutputSpec(TraitedSpec): - G_file = File(desc='Output G file') - a_file = File(desc='Output a file') - b_file = File(desc='Output b file') - - -class PLANETSim(QI.SimCommand): - """ - Simulate ellipse parameters from T1/T2 - - """ - - _cmd = 'qi planet' - _param_files = ['PD', 'T1', 'T2'] - input_spec = PLANETSimInputSpec - output_spec = PLANETSimOutputSpec - - def _list_outputs(self): - outputs = self.output_spec().get() - outputs['G_file'] = path.abspath(self.inputs.G_file) - outputs['a_file'] = path.abspath(self.inputs.a_file) - outputs['b_file'] = path.abspath(self.inputs.b_file) - return outputs diff --git a/Python/T1_mapping.ipynb b/Python/T1_mapping.ipynb index 0f73b72b..f0fc7e3d 100644 --- a/Python/T1_mapping.ipynb +++ b/Python/T1_mapping.ipynb @@ -10,7 +10,7 @@ "\n", "This pipeline (http://github.com/spinicist/nanslice) downloads the BrainWeb brain & B1 phantoms from http://brainweb.bic.mni.mcgill.ca. It then replaces the tissue classification labels with values of Proton Density and T1, simulates an SPGR/FLASH image with some added noise, and finally uses that simulated data to fit for T1 and PD again.\n", "\n", - "In order to display the images, and to do some " + "In order to display the images, please install the `nanslice` package from `pip`." ] }, { @@ -22,25 +22,22 @@ }, { "cell_type": "code", - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ "%matplotlib inline\n", - "from QUIT.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2\n", - "from nanslice import Layer\n", + "from qipype.interfaces.relax import DESPOT1, DESPOT1Sim\n", + "from qipype.sims import init_brainweb, make_phantom\n", "import nanslice.jupyter as ns\n", - "import nibabel as nib\n", - "import numpy as np\n", - "import requests\n", - "import gzip\n", - "import os.path" - ] + "import nibabel as nib" + ], + "execution_count": null }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ "## Create Phantom Data\n", "\n", @@ -53,64 +50,25 @@ "source": [ "# Data Fetching\n", "\n", - "We now download the Brainweb brain & T1 phantoms, if they are not already present in the current directory. These are stored in MINC format, so we use " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Display the phantom data as a check" - ] - }, - { - "cell_type": "code", - "metadata": {}, - "outputs": [], - "source": [ - "if not isfile('classes.mnc'):\n", - " print('Downloading classes')\n", - " params = {'download_for_real':'[Start Download!]',\n", - " 'do_download_alias':'phantom_1.0mm_normal_crisp',\n", - " 'format_value':'minc',\n", - " 'who_name': 'Tobias Wood',\n", - " 'who_institution': 'KCL',\n", - " 'who_email': 'tobias.wood@kcl.ac.uk'}\n", - " response = requests.get(url='http://brainweb.bic.mni.mcgill.ca/cgi/brainweb1', params=params)\n", - " minc_file = open('classes.mnc', 'wb')\n", - " minc_file.write(response.content)\n", - "if not isfile('rf20_C.mnc'):\n", - " print('Downloading B1')\n", - " params = {'download_for_real':'[Start Download!]',\n", - " 'do_download_alias':'rf20_C.mnc.gz',\n", - " 'format_value':'minc',\n", - " 'zip_value':'none',\n", - " 'who_name': 'Tobias Wood',\n", - " 'who_institution': 'KCL',\n", - " 'who_email': 'tobias.wood@kcl.ac.uk'}\n", - " response = requests.get(url='https://brainweb.bic.mni.mcgill.ca/cgi/brainweb1', params=params)\n", - " b1_file = open('rf20_C.mnc', 'wb')\n", - " b1_file.write(response.content)\n", - "classes = nanslice.Layer('classes.mnc')\n", - "b0_minc = nanslice.Layer('rf20_B.mnc')\n", - "b1_minc = nanslice.Layer('rf20_C.mnc')" + "We now download the BrainWeb brain & T1 phantoms, if they are not already present in the current directory. QUIT provides a utility function for this." ] }, { "cell_type": "code", - "metadata": {}, + "execution_count": null, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ - "ns.three_plane(classes)" + "init_brainweb()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now we define a function that allows us to convert the tissue classes to T1 and PD/M0 values, and use it to create our reference T1, PD and B1 images. These will be saved into the current directory. This function also allows us to subsample the images, which is useful when investigating methods that take longer to fit.\n", - "\n", - "The BrainWeb phantom uses the following mapping:\n", + "Now we use another QUIT function to convert the BrainWeb phantoms (which are tissue class maps) into parameter maps, in this case Proton Density and T1. There are 10 tissue classes:\n", "\n", "0. Background\n", "1. CSF\n", @@ -123,32 +81,23 @@ "8. Glial Matter\n", "9. Connective\n", "\n", - "To keep things simple, we only specify values for CSF, GM & WM, and set the other tissue types to zero. The B1 data is used as is from the phantom." + "The `make_phantom` function requires a dictionary, where every parameter map we want to make is a key, and the value is a 10 element array with the values for each class. To keep things simple here, we only specify values for CSF, GM & WM, and set the other tissue types to zero. We also specify that we want a B1 map for later use, and sub-sample the images by a factor of 2 to make it faster." ] }, { "cell_type": "code", - "metadata": {}, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ - "def classes_to_params(M0Vals, T1Vals, subsamp=1):\n", - " class_data = classes.image.get_data().astype('int32')\n", - " M0data = np.choose(class_data[::subsamp,::subsamp,::subsamp], M0vals).astype('float32')\n", - " T1data = np.choose(class_data[::subsamp,::subsamp,::subsamp], T1vals).astype('float32')\n", - " T2data = np.choose(class_data[::subsamp,::subsamp,::subsamp], T2vals).astype('float32')\n", - " B1data = b1_minc.image.get_data().astype('float32')[::subsamp,::subsamp,::subsamp]\n", - " # PDdata = np.array(list(map(PDFunc, classes.image.get_data())))\n", - " M0image = nib.nifti1.Nifti1Image(M0data, affine=classes.image.affine)\n", - " T1image = nib.nifti1.Nifti1Image(T1data, affine=classes.image.affine)\n", - " B1image = nib.nifti1.Nifti1Image(B1data, affine=classes.image.affine)\n", - " nib.save(M0image, 'M0.nii.gz')\n", - " nib.save(T1image, 'T1.nii.gz')\n", - " nib.save(B1image, 'B1.nii.gz')\n", - "\n", - "M0vals = np.array([0, 1, 0.8, 0.7, 0, 0, 0, 0, 0, 0])\n", - "T1vals = np.array([0, 3.0, 1.3, 0.9, 0, 0, 0, 0, 0, 0])\n", - "classes_to_params(M0vals, T1vals, 1)" - ] + "param_dict = {'PD': [0, 1, 0.8, 0.7, 0, 0, 0, 0, 0, 0],\n", + " 'T1': [0, 3.0, 1.3, 0.9, 0, 0, 0, 0, 0, 0]}\n", + "make_phantom(param_dict, B1=True, subsamp=2)\n", + "display(ns.three_plane('T1.nii.gz'))\n", + "display(ns.three_plane('B1.nii.gz'))" + ], + "execution_count": null }, { "cell_type": "markdown", @@ -169,13 +118,12 @@ "outputs": [], "source": [ "sequence_parameters = {'SPGR': {'TR': 10e-3, 'FA': [3, 18]}}\n", - "simulator_results = DESPOT1Sim(sequence=d1seq, in_file='sim_spgr.nii.gz', noise=0.001, PD='PD.nii.gz', T1='T1.nii.gz').run()\n", + "simulator_results = DESPOT1Sim(sequence=sequence_parameters, out_file='sim_spgr.nii.gz', noise=0.001, PD='PD.nii.gz', T1='T1.nii.gz').run()\n", "\n", - "simulated_spgr = Layer(simulator_results.outputs.out_file, volume=0)\n", - "display(ns.three_plane(spgr))\n", - "spgr.volume = 1\n", - "display(ns.three_plane(spgr))" - ] + "display(ns.three_plane(simulator_results.outputs.out_file, volume=0))\n", + "display(ns.three_plane(simulator_results.outputs.out_file, volume=1))" + ], + "execution_count": null }, { "cell_type": "markdown", @@ -186,25 +134,85 @@ "Now we have simulated data we can run the DESPOT1/VFA fitting code and compare the results to our phantom reference data." ] }, + { + "cell_type": "code", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fit = DESPOT1(sequence=sequence_parameters, in_file=simulator_results.outputs.out_file)\n", + "print(fit)\n", + "print(dir(fit))\n", + "fitting_results = fit.run()\n", + "\n", + "display(ns.compare('T1.nii.gz', fitting_results.outputs.t1_map, title='T1 Comparison'))\n", + "display(ns.compare('PD.nii.gz', fitting_results.outputs.pd_map, title='PD Comparison'))" + ], + "execution_count": null + }, { "cell_type": "markdown", "metadata": {}, "source": [ "## B1 Mapping\n", "\n", - "SPGR images are affected by B1 inhomogeneity. In the above simulation, we did not specify a B1 map and assumed that B1 was flat." + "SPGR images are affected by B1 inhomogeneity. In the above simulation, we did not specify a B1 map and assumed that B1 was flat, but this is not true in real images. We will now add B1 into both our simulation and fitting.\n", + "\n", + "For DESPOT1, and most other methods, B1 is a \"fixed\" parameter, in contrast to T1 & PD which are \"varying\" parameters. This means that DESPOT1 does not fit a value for B1, it must be measured with another scan and is then fixed within each voxel. There are two methods in QUIT for calculating B1 maps (see https://quit.readthedocs.io/en/latest/Docs/B1.html).\n", + "\n", + "First we simulate our input images again, supplying the B1 map, but fit without supplying the B1 map. We add the `prefix` parameter to the fit to avoid overwriting the output of the previous fitting run." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "b1_sim_results = DESPOT1Sim(sequence=sequence_parameters, out_file='sim_b1.nii.gz', b1_map='B1.nii.gz', noise=0.001, PD='PD.nii.gz', T1='T1.nii.gz').run()\n", + "no_b1_fit_results = DESPOT1(prefix='no_b1', sequence=sequence_parameters, in_file=b1_sim_results.outputs.out_file).run()\n", + "display(ns.compare('T1.nii.gz', no_b1_fit_results.outputs.t1_map, clim=(0.8,1.5), title='Phantom Vs No-B1 in Fitting'))" ] }, { "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "print(no_b1_fit_results.outputs.t1_map)" + ] + }, + { + "cell_type": "markdown", "metadata": {}, + "source": [ + "Finally, we fit again, this time supplying the B1 map, and compare the results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [] + }, "outputs": [], "source": [ - "fitting_results = DESPOT1(sequence=d1seq, in_file='sim_spgr.nii.gz').run()\n", - "\n", - "display(ns.compare('T1.nii.gz', fitting_results.outputs.t1_map, title='T1 Comparison'))\n", - "display(ns.compare('PD.nii.gz', fitting_results.outputs.pd_map, title='PD Comparison'))" + "b1_fit_results = DESPOT1(prefix='b1', sequence=sequence_parameters, in_file=b1_sim_results.outputs.out_file, b1_map='B1.nii.gz').run()\n", + "display(ns.compare(b1_fit_results.outputs.t1_map, no_b1_fit_results.outputs.t1_map, title='B1 in fitting versus no B1'))" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -226,7 +234,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.4" + "version": "3.7.6-final" }, "nteract": { "version": "0.8.4" diff --git a/Python/Tests/test_despot_sc.py b/Python/Tests/test_despot_sc.py index 49d23283..32e2daac 100644 --- a/Python/Tests/test_despot_sc.py +++ b/Python/Tests/test_despot_sc.py @@ -1,13 +1,22 @@ +from pathlib import Path +from os import chdir import unittest from nipype.interfaces.base import CommandLine -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2, DESPOT2Sim, HIFI, HIFISim, FM, FMSim +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2, DESPOT2Sim, HIFI, HIFISim, FM, FMSim vb = True CommandLine.terminal_output = 'allatonce' class DESPOT_SC(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_despot1(self): seq = {'SPGR': {'TR': 10e-3, 'FA': [3, 18]}} spgr_file = 'sim_spgr.nii.gz' @@ -19,9 +28,9 @@ def test_despot1(self): NewImage(img_size=img_sz, grad_dim=1, grad_vals=( 0.8, 1.3), out_file='T1.nii.gz', verbose=vb).run() - DESPOT1Sim(sequence=seq, in_file=spgr_file, + DESPOT1Sim(sequence=seq, out_file=spgr_file, noise=noise, verbose=vb, - PD='PD.nii.gz', T1='T1.nii.gz').run() + PD_map='PD.nii.gz', T1_map='T1.nii.gz').run() DESPOT1(sequence=seq, in_file=spgr_file, verbose=vb, residuals=True).run() @@ -50,7 +59,7 @@ def test_hifi(self): HIFISim(sequence=seqs, spgr_file=spgr_file, mprage_file=mprage_file, noise=noise, verbose=vb, - PD='PD.nii.gz', T1='T1.nii.gz', B1='B1.nii.gz').run() + PD_map='PD.nii.gz', T1_map='T1.nii.gz', B1_map='B1.nii.gz').run() HIFI(sequence=seqs, spgr_file=spgr_file, mprage_file=mprage_file, verbose=vb, residuals=True).run() @@ -79,11 +88,11 @@ def test_despot2(self, gs=False, tol=20): NewImage(img_size=img_sz, grad_dim=2, grad_vals=(0.04, 0.1), out_file='T2.nii.gz', verbose=vb).run() - DESPOT2Sim(sequence=seq, in_file=ssfp_file, - t1_file='T1.nii.gz', ellipse=gs, noise=noise, verbose=vb, - PD='PD.nii.gz', T2='T2.nii.gz').run() + DESPOT2Sim(sequence=seq, out_file=ssfp_file, + ellipse=gs, noise=noise, verbose=vb, + PD_map='PD.nii.gz', T2_map='T2.nii.gz', T1_map='T1.nii.gz').run() DESPOT2(sequence=seq, in_file=ssfp_file, - t1_file='T1.nii.gz', ellipse=gs, verbose=vb, residuals=True).run() + T1_map='T1.nii.gz', ellipse=gs, verbose=vb, residuals=True).run() diff_T2 = Diff(in_file='D2_T2.nii.gz', baseline='T2.nii.gz', noise=noise, verbose=vb).run() @@ -113,11 +122,12 @@ def test_fm(self): NewImage(img_size=img_sz, grad_dim=2, grad_vals=(-100, 100), out_file='f0.nii.gz', verbose=vb).run() - FMSim(sequence=seq, in_file=ssfp_file, asym=False, - t1_file='T1.nii.gz', noise=noise, verbose=vb, - PD='PD.nii.gz', T2='T2.nii.gz', f0='f0.nii.gz').run() + sim = FMSim(sequence=seq, out_file=ssfp_file, asym=False, + T1_map='T1.nii.gz', noise=noise, verbose=vb, + PD_map='PD.nii.gz', T2_map='T2.nii.gz', f0_map='f0.nii.gz') + sim.run() FM(sequence=seq, in_file=ssfp_file, asym=False, - t1_file='T1.nii.gz', verbose=vb, residuals=True).run() + T1_map='T1.nii.gz', verbose=vb, residuals=True).run() diff_T2 = Diff(in_file='FM_T2.nii.gz', baseline='T2.nii.gz', noise=noise, verbose=vb).run() diff --git a/Python/Tests/test_mt.py b/Python/Tests/test_mt.py index c11a45c1..a270fa2d 100644 --- a/Python/Tests/test_mt.py +++ b/Python/Tests/test_mt.py @@ -1,14 +1,23 @@ +from pathlib import Path +from os import chdir import unittest import numpy as np from nipype.interfaces.base import CommandLine -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.mt import Lorentzian, LorentzianSim, Lineshape, qMT, qMTSim, ZSpec +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.mt import Lorentzian, LorentzianSim, Lineshape, qMT, qMTSim, ZSpec vb = True CommandLine.terminal_output = 'allatonce' class MT(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_lorentzian(self): sequence = {'MTSat': {'pulse': {'p1': 0.4, 'p2': 0.3, @@ -34,13 +43,14 @@ def test_lorentzian(self): NewImage(out_file='A.nii.gz', verbose=vb, img_size=img_sz, grad_dim=2, grad_vals=(0.5, 1)).run() - LorentzianSim(sequence=sequence, pools=pools, in_file=lorentz_file, - noise=noise, verbose=vb, - DS_f0='f0.nii.gz', - DS_fwhm='fwhm.nii.gz', - DS_A='A.nii.gz').run() - Lorentzian(sequence=sequence, pools=pools, - in_file=lorentz_file, verbose=vb).run() + L1Sim = LorentzianSim(pools) + L1Sim(sequence=sequence, out_file=lorentz_file, + noise=noise, verbose=vb, + DS_f0_map='f0.nii.gz', + DS_fwhm_map='fwhm.nii.gz', + DS_A_map='A.nii.gz').run() + L1Fit = Lorentzian(pools) + L1Fit(sequence=sequence, in_file=lorentz_file, verbose=vb).run() diff_f0 = Diff(in_file='LTZ_DS_f0.nii.gz', baseline='f0.nii.gz', noise=noise, verbose=vb).run() @@ -94,16 +104,17 @@ def test_lorentzian2(self): NewImage(out_file='MTA.nii.gz', verbose=vb, img_size=img_sz, fill=0.4).run() - LorentzianSim(sequence=sequence, pools=pools, in_file=lorentz_file, - noise=noise, verbose=vb, - DS_f0='f0.nii.gz', - DS_fwhm='fwhm.nii.gz', - DS_A='A.nii.gz', - MT_fwhm='MTfwhm.nii.gz', - MT_f0='f0.nii.gz', - MT_A='MTA.nii.gz').run() - Lorentzian(sequence=sequence, pools=pools, - in_file=lorentz_file, verbose=vb).run() + L2Sim = LorentzianSim(pools) + L2Sim(sequence=sequence, out_file=lorentz_file, + noise=noise, verbose=vb, + DS_f0_map='f0.nii.gz', + DS_fwhm_map='fwhm.nii.gz', + DS_A_map='A.nii.gz', + MT_fwhm_map='MTfwhm.nii.gz', + MT_f0_map='f0.nii.gz', + MT_A_map='MTA.nii.gz').run() + L2Fit = Lorentzian(pools) + L2Fit(sequence=sequence, in_file=lorentz_file, verbose=vb).run() diff_fwhm = Diff(in_file='LTZ_DS_fwhm.nii.gz', baseline='fwhm.nii.gz', noise=noise, verbose=vb).run() @@ -153,10 +164,14 @@ def test_qMT(self): NewImage(out_file=t1app, verbose=vb, img_size=img_sz, fill=1.0).run() - qMTSim(sequence=qmt, in_file=qmt_file, t1_map=t1app, lineshape=lineshape_file, + qMTSim(sequence=qmt, out_file=qmt_file, T1_map=t1app, lineshape=lineshape_file, noise=noise, verbose=vb, - M0_f='M0_f.nii.gz', F_over_R1_f='F_over_R1_f.nii.gz', T2_b='T2_b.nii.gz', T1_f_over_T2_f='T1_f_over_T2_f.nii.gz', k='k.nii.gz').run() - qMT(sequence=qmt, in_file=qmt_file, t1_map=t1app, + M0_f_map='M0_f.nii.gz', + F_over_R1_f_map='F_over_R1_f.nii.gz', + T2_b_map='T2_b.nii.gz', + T1_f_over_T2_f_map='T1_f_over_T2_f.nii.gz', + k_map='k.nii.gz').run() + qMT(sequence=qmt, in_file=qmt_file, T1_map=t1app, lineshape=lineshape_file, verbose=vb).run() diff_PD = Diff(in_file='QMT_M0_f.nii.gz', baseline='M0_f.nii.gz', diff --git a/Python/Tests/test_perfusion.py b/Python/Tests/test_perfusion.py index e47ce29b..6c06b6cc 100644 --- a/Python/Tests/test_perfusion.py +++ b/Python/Tests/test_perfusion.py @@ -1,14 +1,23 @@ +from pathlib import Path +from os import chdir import unittest from math import sqrt from nipype.interfaces.base import CommandLine -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.perfusion import ASL, ASE, ASESim, ZShim +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.perfusion import ASL, ASE, ASESim, ASEDBVSim, ZShim vb = True CommandLine.terminal_output = 'allatonce' class Perfusion(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_asl(self): seq = {'CASL': {'TR': 4.0, 'label_time': 3.0, 'post_label_delay': [0.3]}} @@ -42,11 +51,11 @@ def test_oef(self): NewImage(img_size=img_sz, grad_dim=2, grad_vals=(0.005, 0.025), out_file='DBV.nii.gz', verbose=vb).run() - ASESim(sequence=seq, in_file=ase_file, noise=noise, verbose=vb, - S0='S0.nii.gz', - dT='dT.nii.gz', - R2p='R2p.nii.gz', - DBV='DBV.nii.gz').run() + ASEDBVSim(sequence=seq, out_file=ase_file, noise=noise, verbose=vb, + S0_map='S0.nii.gz', + dT_map='dT.nii.gz', + R2p_map='R2p.nii.gz', + DBV_map='DBV.nii.gz').run() ASE(sequence=seq, in_file=ase_file, verbose=vb).run() diff_R2p = Diff(in_file='ASE_R2p.nii.gz', baseline='R2p.nii.gz', @@ -72,11 +81,11 @@ def test_oef_fixed_dbv(self): NewImage(img_size=img_sz, grad_dim=1, grad_vals=(1.0, 3.0), out_file='R2p.nii.gz', verbose=vb).run() - ASESim(sequence=seq, in_file=ase_file, + ASESim(sequence=seq, out_file=ase_file, fix_DBV=DBV, noise=noise, verbose=vb, - S0='S0.nii.gz', - dT='dT.nii.gz', - R2p='R2p.nii.gz').run() + S0_map='S0.nii.gz', + dT_map='dT.nii.gz', + R2p_map='R2p.nii.gz').run() ASE(sequence=seq, in_file=ase_file, fix_DBV=DBV, verbose=vb).run() diff_R2p = Diff(in_file='ASE_R2p.nii.gz', baseline='R2p.nii.gz', diff --git a/Python/Tests/test_relax.py b/Python/Tests/test_relax.py index 86be8e2e..4e132b6b 100644 --- a/Python/Tests/test_relax.py +++ b/Python/Tests/test_relax.py @@ -1,14 +1,22 @@ +from pathlib import Path +from os import chdir import unittest from nipype.interfaces.base import CommandLine -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.relax import Multiecho, MultiechoSim -from QUIT.interfaces.mt import Lineshape +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.relax import Multiecho, MultiechoSim vb = True CommandLine.terminal_output = 'allatonce' class Relax(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_multiecho(self): me = {'MultiEcho': {'TR': 10, 'TE1': 0.01, 'ESP': 0.01, 'ETL': 5}} @@ -21,8 +29,8 @@ def test_multiecho(self): NewImage(img_size=img_sz, grad_dim=2, grad_vals=(0.04, 0.1), out_file='T2.nii.gz', verbose=vb).run() - MultiechoSim(sequence=me, in_file=me_file, - PD='PD.nii.gz', T2='T2.nii.gz', + MultiechoSim(sequence=me, out_file=me_file, + PD_map='PD.nii.gz', T2_map='T2.nii.gz', noise=noise, verbose=vb).run() Multiecho(sequence=me, in_file=me_file, verbose=vb).run() diff --git a/Python/Tests/test_ssfp.py b/Python/Tests/test_ssfp.py index e6effd7f..f716628f 100644 --- a/Python/Tests/test_ssfp.py +++ b/Python/Tests/test_ssfp.py @@ -1,15 +1,24 @@ import unittest import numpy as np +from pathlib import Path +from os import chdir from nipype.interfaces.base import CommandLine, DynamicTraitedSpec -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.relax import Ellipse, EllipseSim, PLANET, PLANETSim -from QUIT.interfaces.mt import eMT, eMTSim +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.relax import Ellipse, EllipseSim, PLANET, PLANETSim +from qipype.interfaces.mt import eMT, eMTSim vb = True CommandLine.terminal_output = 'allatonce' class SSFP(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_planet(self): ellipse_seq = {"SSFP": { "FA": [15, 15, 15, 15, 15, 15], @@ -42,12 +51,13 @@ def test_planet(self): PLANETSim(sequence=planet_seq, G_file=planet_G, a_file=planet_a, b_file=planet_b, noise=0, verbose=vb, - PD='PD.nii.gz', - T1='T1.nii.gz', - T2='T2.nii.gz').run() - EllipseSim(sequence=ellipse_seq, in_file=ellipse_file, + PD_map='PD.nii.gz', + T1_map='T1.nii.gz', + T2_map='T2.nii.gz').run() + EllipseSim(sequence=ellipse_seq, out_file=ellipse_file, noise=noise, verbose=vb, - G=planet_G, a=planet_a, b=planet_b, theta_0='zero.nii.gz', phi_rf='zero.nii.gz').run() + G_map=planet_G, a_map=planet_a, b_map=planet_b, + theta_0_map='zero.nii.gz', phi_rf_map='zero.nii.gz').run() Ellipse(sequence=ellipse_seq, in_file=ellipse_file, verbose=vb).run() PLANET(sequence=planet_seq, G_map=planet_G, a_map=planet_a, b_map=planet_b, verbose=vb).run() @@ -113,14 +123,14 @@ def test_emt(self): eMTSim(sequence=emt_seq, G_file=emt_G, a_file=emt_a, b_file=emt_b, noise=0, verbose=vb, - PD='PD.nii.gz', - T1_f='T1_f.nii.gz', - T2_f='T2_f.nii.gz', - f_b='f_b.nii.gz', - k_bf='k_bf.nii.gz').run() - EllipseSim(sequence=ellipse_sim, in_file=ellipse_file, + PD_map='PD.nii.gz', + T1_f_map='T1_f.nii.gz', + T2_f_map='T2_f.nii.gz', + f_b_map='f_b.nii.gz', + k_bf_map='k_bf.nii.gz').run() + EllipseSim(sequence=ellipse_sim, out_file=ellipse_file, noise=noise, verbose=vb, - G=emt_G, a=emt_a, b=emt_b, theta_0='zero.nii.gz', phi_rf='zero.nii.gz').run() + G_map=emt_G, a_map=emt_a, b_map=emt_b, theta_0_map='zero.nii.gz', phi_rf_map='zero.nii.gz').run() Ellipse(sequence=ellipse_fit, in_file=ellipse_file, verbose=vb).run() eMT(sequence=emt_seq, G_map=emt_G, a_map=emt_a, b_map=emt_b, verbose=vb).run() diff --git a/Python/Tests/test_utils.py b/Python/Tests/test_utils.py index f065f36d..fcc0b0f6 100644 --- a/Python/Tests/test_utils.py +++ b/Python/Tests/test_utils.py @@ -1,14 +1,23 @@ +from pathlib import Path +from os import chdir import unittest from math import sqrt from nipype.interfaces.base import CommandLine -from QUIT.interfaces.core import NewImage, Diff -from QUIT.interfaces.utils import PolyImage, PolyFit, Filter, RFProfile +from qipype.interfaces.core import NewImage, Diff +from qipype.interfaces.utils import PolyImage, PolyFit, Filter, RFProfile vb = True CommandLine.terminal_output = 'allatonce' class Utils(unittest.TestCase): + def setUp(self): + Path('testdata').mkdir(exist_ok=True) + chdir('testdata') + + def tearDown(self): + chdir('../') + def test_poly(self): sz = 16 scale = sqrt(3 * (sz/2)**2) diff --git a/Python/example_notebook.ipynb b/Python/example_notebook.ipynb index f813a41d..63ddfa95 100644 --- a/Python/example_notebook.ipynb +++ b/Python/example_notebook.ipynb @@ -26,7 +26,7 @@ "outputs": [], "source": [ "%matplotlib inline\n", - "from QUIT.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2\n", + "from quit.interfaces.relax import DESPOT1, DESPOT1Sim, DESPOT2\n", "from nanslice import Layer\n", "import nanslice.jupyter as ns\n", "import nibabel as nib\n", diff --git a/Python/qipype/__init__.py b/Python/qipype/__init__.py new file mode 100644 index 00000000..b78f4581 --- /dev/null +++ b/Python/qipype/__init__.py @@ -0,0 +1,4 @@ +#!/usr/bin/env python +""" +nipype interfaces and workflows for QUIT. +""" diff --git a/Python/QUIT/base.py b/Python/qipype/base.py similarity index 59% rename from Python/QUIT/base.py rename to Python/qipype/base.py index 167838f0..c36e8417 100644 --- a/Python/QUIT/base.py +++ b/Python/qipype/base.py @@ -30,9 +30,7 @@ class InputSpec(InputBaseSpec): """ Input Specification for most QUIT tools """ - # Most QUIT programs take similar arguments json = traits.File(exists=True, desc='JSON Input file', argstr='--json=%s') - subregion = traits.String( desc='Only process a subregion of the image. Argument should be a string "start_x,start_y,start_z,size_x,size_y,size_z"', argstr='--subregion=%s') threads = traits.Int( @@ -56,60 +54,69 @@ class FitInputSpec(InputSpec): desc='Write out residuals for each data-point', argstr='--resids') -class SimInputBaseSpec(DynamicTraitedSpec): +def SimInputSpec(name, varying, fixed=[], out_files=None, extras=None): """ Input specification for tools in simulation mode """ - ignore_exception = traits.Bool( - False, - usedefault=True, - nohash=True, - deprecated='1.0.0', - desc='Print an error message instead of throwing an exception ' - 'in case the interface fails to run') - args = traits.Str(argstr='%s', desc='Additional parameters to the command') - environ = traits.DictStrStr( - desc='Environment variables', usedefault=True, nohash=True) - # This input does not have a "usedefault=True" so the set_default_terminal_output() - # method would work - terminal_output = traits.Enum( - 'stream', - 'allatonce', - 'file', - 'none', - deprecated='1.0.0', - desc=('Control terminal output: `stream` - ' - 'displays to terminal immediately (default), ' - '`allatonce` - waits till command is ' - 'finished to display output, `file` - ' - 'writes output to file, `none` - output' - ' is ignored'), - nohash=True) - - verbose = traits.Bool(desc='Print more information', argstr='--verbose') - environ = {'QUIT_EXT': 'NIFTI_GZ'} - json = traits.File(exists=True, desc='JSON Input file', argstr='--json=%s') - noise = traits.Float(desc='Noise level to add to simulation', - argstr='--simulate=%f', default_value=0.0, usedefault=True) - subregion = traits.String( - desc='Only process a subregion of the image. Argument should be a string "start_x,start_y,start_z,size_x,size_y,size_z"', argstr='--subregion=%s') - threads = traits.Int( - desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') - mask_file = File( - desc='Only process voxels within the mask', argstr='--mask=%s') - + if out_files is None: + out_files = ['out'] + + attrs = {'noise': traits.Float(desc='Noise level to add to simulation', + argstr='--simulate=%f', default_value=0.0, usedefault=True)} + varying_names = [] + for v in varying: + aname = '{}_map'.format(v) + desc = 'Path to input {} map'.format(v) + attrs[aname] = File(exists=True, desc='Path to {} map'.format(v)) + varying_names.append(aname) + attrs['varying_names'] = varying_names + + for f in fixed: + aname = '{}_map'.format(f) + desc = 'Path to fixed {} map'.format(f) + args = '--{}=%s'.format(f) + attrs[aname] = File(argstr=args, exists=True, desc=desc) + + for idx, o in enumerate(out_files): + aname = '{}_file'.format(o) + desc = 'Output {} file'.format(o) + attrs[aname] = File(argstr='%s', mandatory=True, + position=idx, desc=desc) + if extras: + for k, v in extras.items(): + attrs[k] = v + T = type(name + 'SimInputSpec', (InputSpec,), attrs) + return T -class SimInputSpec(SimInputBaseSpec): - """ - Input specification for tools in simulation mode - """ - in_file = File(argstr='%s', mandatory=True, - position=-1, desc='Output simulated file') ################################# Output Specs ################################# - - -class SimOutputSpec(TraitedSpec): - out_file = File(desc="Path to output image") +# +# QUIT output specifications follow a consistent pattern, so these are defined +# as functions that return a type. Closest thing to C++ templates I could find + + +def FitOutputSpec(prefix, parameters): + attrs = {} + for p in parameters: + pname = '{}_map'.format(p) + fname = '{}_{}.nii.gz'.format(prefix, p) + desc = 'Path to {}'.format(p) + attrs[pname] = File(fname, desc=desc, usedefault=True) + attrs['rmse_map'] = File('{}_rmse.nii.gz'.format( + prefix), desc='Path to RMSE', usedefault=True) + T = type(prefix + 'FitOutputSpec', (TraitedSpec,), attrs) + return T + + +def SimOutputSpec(name, files=None): + if files is None: + files = ['out'] + attrs = {} + for f in files: + aname = '{}_file'.format(f) + desc = 'Path to {} image file'.format(f) + attrs[aname] = File(desc=desc) + T = type(name + 'SimOutputSpec', (TraitedSpec,), attrs) + return T ################################### Commands ################################### @@ -183,7 +190,7 @@ def _gen_fname(self, def _parse_inputs(self, skip=None): # Make sequence dictionary into a .json file for input to interface if hasattr(self, '_json') and not isdefined(self.inputs.json): - fname = '_input.json' + fname = '_' + self.__class__.__name__ + '_input.json' with open(fname, 'w') as outfile: json.dump(self._json, outfile, indent=2) self.inputs.json = fname @@ -195,37 +202,15 @@ class FitCommand(BaseCommand): Support for standard fitting tools. """ - def _add_prefixes(self, outputs): - """ - Add prefix to everything in .outputs - """ - for key, value in outputs.items(): - if isinstance(value, (list,)): - outputs[key] = [self._gen_fname( - x, prefix=self.inputs.prefix) for x in value] - else: - outputs[key] = self._gen_fname( - value, prefix=self.inputs.prefix) - return outputs - def _list_outputs(self): - outputs = {} - if hasattr(self, '_param_files'): - for p in self._param_files: - pname = '{}_file'.format(p) - fname = '{}_{}.nii.gz'.format(self._prefix, p) - outputs[pname] = fname - return self._add_prefixes(outputs) + outputs = self.output_spec().get() + prefixed = {} + for key, trait in outputs.items(): + prefixed[key] = self._gen_fname(trait, prefix=self.inputs.prefix) + return prefixed def __init__(self, sequence={}, **kwargs): self._json = deepcopy(sequence) - if hasattr(self, '_param_files'): - for p in self._param_files: - pname = '{}_file'.format(p) - fname = '{}_{}.nii.gz'.format(self._prefix, p) - desc = 'Path to {}'.format(p) - setattr(self.output_spec, pname, File( - fname, desc=desc, usedefault=True)) super().__init__(**kwargs) @@ -236,21 +221,23 @@ class SimCommand(BaseCommand): def __init__(self, sequence={}, **kwargs): self._json = deepcopy(sequence) - for pf in self._param_files: - setattr(self.input_spec, pf, File( - exists=True, desc='Path to %s map' % pf)) super().__init__(**kwargs) - def _parse_inputs(self, skip=[]): - for pf in self._param_files: - if isdefined(getattr(self.inputs, pf)): - self._json[pf] = getattr(self.inputs, pf) - skip.append(pf) - return super()._parse_inputs(skip) + def _parse_inputs(self, skip=None): + if skip is None: + skip = [] + for v in self.input_spec.varying_names: + if isdefined(getattr(self.inputs, v)): + self._json[v] = getattr(self.inputs, v) + skip.append(v) + parsed = super()._parse_inputs(skip) + return parsed def _list_outputs(self): + inputs = self.inputs.get() outputs = self.output_spec().get() - outputs['out_file'] = path.abspath(self.inputs.in_file) + for k, v in outputs.items(): + outputs[k] = path.abspath(inputs[k]) return outputs diff --git a/Python/qipype/interfaces/b1.py b/Python/qipype/interfaces/b1.py new file mode 100644 index 00000000..43b50144 --- /dev/null +++ b/Python/qipype/interfaces/b1.py @@ -0,0 +1,121 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Implementation of nipype interfaces for QUIT B1 utilities +Requires that the QUIT tools are in your your system path +""" + +from os import path, getcwd +from nipype.interfaces.base import TraitedSpec, File, traits +from .. import base + +############################ qidream ############################ +# Implemented but not tested # + + +class QIDreamInputSpec(base.InputSpec): + # Inputs + dream_file = File(exists=True, argstr='%s', mandatory=True, + position=0, desc='Input file. Must have 2 volumes (FID and STE)') + + # Options + threads = traits.Int( + desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') + prefix = traits.String( + desc='Add a prefix to output filenames', argstr='--out=%s') + order = traits.String( + desc='Volume order - f/s/v - fid/ste/vst first', argstr='--order=%s') + mask_file = File( + desc='Only process voxels within the mask', argstr='--mask=%s') + alpha = traits.Float( + desc="Nominal flip-angle (default 55)", argstr="--alpha=%f") + + +class QIDreamOutputSpec(TraitedSpec): + b1_rel_map = File(desc="The relative flip-angle map.") + b1_act_map = File(desc="The actual achieved angle in each voxel.") + + +class QIDream(base.BaseCommand): + """ + Interface for qidream + + Example 1 + ------- + >>> from qipype.interfaces.relax import QiDream + >>> interface = QiDream(prefix='nipype_', param_file='spgr_params.json') + + """ + + _cmd = 'qi dream' + input_spec = QIDreamInputSpec + output_spec = QIDreamOutputSpec + + # If the command requires a json input file + def _format_arg(self, name, spec, value): + return self._process_params(name, spec, value) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['b1_rel_map'] = os.path.abspath( + self._add_prefix('DREAM_B1.nii.gz')) + outputs['b1_act_map'] = os.path.abspath( + self._add_prefix('DREAM_angle.nii.gz')) + + return outputs + +############################ qiafi ############################ +# Implemented but not tested # + + +class QIAFIInputSpec(base.InputSpec): + # Inputs + afi_file = File(exists=True, argstr='%s', mandatory=True, + position=0, desc='Input file') + + # Options + threads = traits.Int( + desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') + prefix = traits.String( + desc='Add a prefix to output filenames', argstr='--out=%s') + alpha = traits.Float( + desc="Specify nominal flip-angle, default 55", argstr="--flip=%f") + tr_ratio = traits.Float( + desc="Specify TR2:TR1 ratio, default 5", argstr="--ratio=%f") + save_act_b1 = traits.Bool( + desc="Write out the actual flip-angle as well as B1", argstr="--save") + + +class QIAFIOutputSpec(TraitedSpec): + # Specify which outputs there are + b1_rel_map = File(desc="The relative flip-angle map.") + b1_act_map = File(desc="The actual flip-angle map.") + + +class QIAFI(base.BaseCommand): + """ + help for myInterface + + Example 1 + ------- + >>> from qipype.interfaces.relax import QiAfi + >>> interface = QiAfi(prefix='nipype_', param_file='spgr_params.json') + + """ + + _cmd = 'qi afi' + input_spec = QIAFIInputSpec + output_spec = QIAFIOutputSpec + + # If the command requires a json input file + def _format_arg(self, name, spec, value): + return self._process_params(name, spec, value) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['b1_rel_map'] = os.path.abspath( + self._add_prefix('AFI_B1.nii.gz')) + outputs['b1_act_map'] = os.path.abspath( + self._add_prefix('AFI_angle.nii.gz')) + return outputs diff --git a/Python/QUIT/interfaces/core.py b/Python/qipype/interfaces/core.py similarity index 98% rename from Python/QUIT/interfaces/core.py rename to Python/qipype/interfaces/core.py index 54fd10bc..8474b107 100644 --- a/Python/QUIT/interfaces/core.py +++ b/Python/qipype/interfaces/core.py @@ -50,7 +50,7 @@ class NewImage(CommandLine): Example usage ------- - >>> from QUIT.nipype.CoreProgs import QINewImage + >>> from quit.nipype.CoreProgs import QINewImage >>> qinewimage = QINewImage(out_file='test.nii', imsize=256) >>> sim_res = qinewimage.run() >>> print(sim_res.outputs) diff --git a/Python/QUIT/interfaces/fsl.py b/Python/qipype/interfaces/fsl.py similarity index 100% rename from Python/QUIT/interfaces/fsl.py rename to Python/qipype/interfaces/fsl.py diff --git a/Python/QUIT/interfaces/mt.py b/Python/qipype/interfaces/mt.py similarity index 56% rename from Python/QUIT/interfaces/mt.py rename to Python/qipype/interfaces/mt.py index c6ecdb5e..e9db011b 100644 --- a/Python/QUIT/interfaces/mt.py +++ b/Python/qipype/interfaces/mt.py @@ -8,12 +8,12 @@ from os import path from nipype.interfaces.base import TraitedSpec, DynamicTraitedSpec, File, traits, isdefined -from .. import base as QI +from .. import base ############################ qi_lineshape ############################ -class LineshapeInputSpec(QI.InputBaseSpec): +class LineshapeInputSpec(base.InputBaseSpec): # Options out_file = traits.File(argstr='%s', mandatory=True, exists=False, position=-1, desc='Output file') @@ -33,7 +33,7 @@ class LineshapeOutputSpec(TraitedSpec): out_file = File(desc='JSON Lineshape file') -class Lineshape(QI.BaseCommand): +class Lineshape(base.BaseCommand): """ Pre-calculate lineshapes and write them out for use with qMT """ @@ -50,91 +50,82 @@ def _list_outputs(self): ############################ qi_lorentzian ############################ -class LorentzianInputSpec(QI.FitInputSpec): +class LorentzianInputSpec(base.FitInputSpec): # Options - pools = traits.Int(argstr='--pools=%d', - desc='Number of pools, set from pool arg') additive = traits.Bool(argstr='--add', desc='Use an additive instead of subtractive model') Zref = traits.Float(argstr='--zref=%f', desc='Set reference Z-spectrum value (usually 1 or 0)') -class LorentzianOutputSpec(DynamicTraitedSpec): - rmse_map = File('LTZ_rmse.nii.gz', - desc='Path to residual map', usedefault=True) - - -class Lorentzian(QI.FitCommand): +def Lorentzian(pools): """ - Fit a Lorentzian function to a Z-spectrum + Returns a type for fitting a Lorentzian function to a Z-spectrum with the specified pools """ + pars = [] + for pool in pools: + poolname = pool['name'] + for par in ['f0', 'fwhm', 'A']: + pname = '{}_{}'.format(poolname, par) + pars.append(pname) + ospec = base.FitOutputSpec('LTZ', pars) + + def T_init(self, **kwargs): + base.FitCommand.__init__(self, **kwargs) + self._json['pools'] = pools - _cmd = 'qi lorentzian' - input_spec = LorentzianInputSpec - output_spec = LorentzianOutputSpec + attrs = {'_cmd': 'qi lorentzian --pools={}'.format(len(pools)), + 'input_spec': LorentzianInputSpec, + 'output_spec': ospec, + '__init__': T_init} - def __init__(self, pools={}, **kwargs): - super(Lorentzian, self).__init__(pools=len(pools), **kwargs) - self._json['pools'] = pools - for pool in pools: - pn = pool['name'] - setattr(self.output_spec, pn + '_f0', File('LTZ_%s_f0.nii.gz' % - pn, desc='Path to %s Δf0 map' % pn, usedefault=True)) - setattr(self.output_spec, pn + '_fwhm', File('LTZ_%s_fwhm.nii.gz' % - pn, desc='Path to %s FWHM map' % pn, usedefault=True)) - setattr(self.output_spec, pn + '_A', File('LTZ_%s_A.nii.gz' % - pn, desc='Path to %s A map' % pn, usedefault=True)) + name = 'Lorentzian{}Sim'.format(len(pools)) - def _list_outputs(self): - outputs = self.output_spec().get() - for pool in self._json['pools']: - pn = pool['name'] - outputs[pn + '_f0'] = 'LTZ_%s_f0.nii.gz' % pn - outputs[pn + '_fwhm'] = 'LTZ_%s_fwhm.nii.gz' % pn - outputs[pn + '_A'] = 'LTZ_%s_A.nii.gz' % pn - return self._add_prefixes(outputs) + T = type(name, (base.FitCommand,), attrs) + return T -class LorentzianSimInputSpec(QI.SimInputSpec): - # Options - # Options - pools = traits.Int(argstr='--pools=%d', - desc='Number of pools, set from pool arg') - additive = traits.Bool(argstr='--add', - desc='Use an additive instead of subtractive model') - Zref = traits.Float(argstr='--zref=%f', - desc='Set reference Z-spectrum value (usually 1 or 0)') +def LorentzianSim(pools): + """ + Returns a type for simulating a Lorentzian Z-spectrum with the specified pools + """ + pars = [] + for pool in pools: + poolname = pool['name'] + for par in ['f0', 'fwhm', 'A']: + pname = '{}_{}'.format(poolname, par) + pars.append(pname) + ispec = base.SimInputSpec('LTZ', pars, extras={'additive': traits.Bool(argstr='--add', + desc='Use an additive instead of subtractive model'), + 'Zref': traits.Float(argstr='--zref=%f', + desc='Set reference Z-spectrum value (usually 1 or 0)')}) + ospec = base.SimOutputSpec('LTZ') + + def T_init(self, **kwargs): + base.SimCommand.__init__(self, **kwargs) + self._json['pools'] = pools + attrs = {'_cmd': 'qi lorentzian --pools={}'.format(len(pools)), + 'input_spec': ispec, + 'output_spec': ospec, + '__init__': T_init} -class LorentzianSim(QI.SimCommand): - _cmd = 'qi lorentzian' - input_spec = LorentzianSimInputSpec - output_spec = QI.SimOutputSpec + name = 'Lorentzian{}Sim'.format(len(pools)) - def __init__(self, pools={}, **kwargs): - self._param_files = [] - for pool in pools: - pn = pool['name'] - self._param_files.append(pn + '_f0') - self._param_files.append(pn + '_fwhm') - self._param_files.append(pn + '_A') - super().__init__(pools=len(pools), **kwargs) - self._json['pools'] = pools + T = type(name, (base.SimCommand,), attrs) + return T -############################ qi_qmt ############################ +############################ qi_qmt ############################ -class qMTInputSpec(QI.FitInputSpec): - # Inputs - t1_map = File(exists=True, argstr='%s', mandatory=True, - position=-2, desc='Path to T1 map') - # Options +class qMTInputSpec(base.FitInputSpec): + T1_map = File(exists=True, argstr='--T1=%s', + mandatory=True, desc='Path to T1 map') + f0_map = File(desc='f0 map (Hertz)', argstr='--f0=%s', exists=True) + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') lineshape = traits.String(argstr='--lineshape=%s', mandatory=True, desc='Gauss/Lorentzian/SuperLorentzian/path to JSON file') - f0_map = File(desc='f0 map (Hertz)', argstr='--f0=%s', exists=True) - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') class qMTOutputSpec(TraitedSpec): @@ -153,38 +144,36 @@ class qMTOutputSpec(TraitedSpec): desc="Path to residual map", usedefault=True) -class qMT(QI.FitCommand): +class qMT(base.FitCommand): """ - Fit the Ramani model to a Z-spectrum + Fit the Ramani model to a Z-spectrum. The output parameters here are the interesting (derived) + parameters, not the intrinsic model parameters. """ _cmd = 'qi qmt' input_spec = qMTInputSpec - output_spec = qMTOutputSpec + output_spec = base.FitOutputSpec( + 'QMT', ['PD', 'T1_f', 'T2_f', 'k_bf', 'f_b']) -class qMTSimInputSpec(QI.SimInputSpec): - # Inputs - t1_map = File(argstr='%s', mandatory=True, - position=-2, desc='Path to T1 map') - - # Options - lineshape = traits.String(argstr='--lineshape=%s', mandatory=True, - desc='Gauss/Lorentzian/SuperLorentzian/path to JSON file') - f0_map = File(desc='f0 map (Hertz)', argstr='--f0=%s', exists=True) - b1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') - - -class qMTSim(QI.SimCommand): +class qMTSim(base.SimCommand): + """ + Simulate a Z-spectrum using the Ramani model. Note that the parameters here are the intrinsic + model parameters, which are different to the derived parameters output from the fitting + """ _cmd = 'qi qmt' - _param_files = ['M0_f', 'F_over_R1_f', 'T2_b', 'T1_f_over_T2_f', 'k'] - input_spec = qMTSimInputSpec - output_spec = QI.SimOutputSpec + input_spec = base.SimInputSpec('qMT', + varying=['M0_f', 'F_over_R1_f', + 'T2_b', 'T1_f_over_T2_f', 'k'], + fixed=['f0', 'B1', 'T1'], + extras={'lineshape': traits.String(argstr='--lineshape=%s', mandatory=True, + desc='Gauss/Lorentzian/SuperLorentzian/path to JSON file')}) + output_spec = base.SimOutputSpec('qMT') ############################ qi_zspec ############################ -class ZSpecInputSpec(QI.InputSpec): +class ZSpecInputSpec(base.InputSpec): # Input nifti in_file = File(exists=True, argstr='%s', mandatory=True, position=-1, desc='Input file') @@ -204,7 +193,7 @@ class ZSpecOutputSpec(TraitedSpec): out_file = File(desc="Path to interpolated Z-spectrum/MTA-spectrum") -class ZSpec(QI.BaseCommand): +class ZSpec(base.BaseCommand): """ Interpolate a Z-spectrum (with correction for off-resonance) @@ -230,7 +219,7 @@ def _list_outputs(self): ############################ qi_ssfp_emt ############################ -class eMTInputSpec(QI.InputSpec): +class eMTInputSpec(base.InputSpec): # Inputs G_map = File(exists=True, argstr='%s', mandatory=True, position=-3, desc='Path to G parameter map') @@ -240,26 +229,12 @@ class eMTInputSpec(QI.InputSpec): position=-1, desc='Path to b parameter map') # Options - b1map_file = File(desc='B1 map (ratio) file', argstr='--B1=%s') - f0map_file = File(desc='f0 map (Hertz) file', argstr='--f0=%s') + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + f0_map = File(desc='f0 map (Hertz) file', argstr='--f0=%s') T2b = traits.Float(desc='T2 of bound pool', argstr='--T2b=%f') -class eMTOutputSpec(TraitedSpec): - PD_map = File('EMT_PD.nii.gz', desc="Path to PD map", usedefault=True) - f_b_map = File('EMT_f_b.nii.gz', - desc="Path to bound-pool fraction map", usedefault=True) - k_bf_map = File('EMT_k_bf.nii.gz', - desc="Path to exchange rate map", usedefault=True) - T1_f_map = File('EMT_T1_f.nii.gz', - desc="Path to free-pool T1 map", usedefault=True) - T2_f_map = File('EMT_T2_f.nii.gz', - desc="Path to free-pool T2 map", usedefault=True) - rmse_map = File('EMT_rmse.nii.gz', - desc="Path to residual map", usedefault=True) - - -class eMT(QI.FitCommand): +class eMT(base.FitCommand): """ Fit MT parameters to ellipse parameters @@ -267,50 +242,29 @@ class eMT(QI.FitCommand): _cmd = 'qi ssfp_emt' input_spec = eMTInputSpec - output_spec = eMTOutputSpec - - -class eMTSimInputSpec(QI.SimInputBaseSpec): - # Options - b1map_file = File(desc='B1 map (ratio) file', argstr='--B1=%s') - f0map_file = File(desc='f0 map (Hertz) file', argstr='--f0=%s') - T2b = traits.Float(desc='T2 of bound pool', argstr='--T2b=%f') - G_file = File(argstr='%s', mandatory=True, - position=-3, desc='Output G file') - a_file = File(argstr='%s', mandatory=True, - position=-2, desc='Output a file') - b_file = File(argstr='%s', mandatory=True, - position=-1, desc='Output b file') - + output_spec = base.FitOutputSpec( + 'EMT', ['PD', 'T1_f', 'T2_f', 'k_bf', 'f_b']) -class eMTSimOutputSpec(TraitedSpec): - G_file = File(desc='Output G file') - a_file = File(desc='Output a file') - b_file = File(desc='Output b file') - -class eMTSim(QI.SimCommand): +class eMTSim(base.SimCommand): """ Simulate ellipse parameters from MT parameters """ _cmd = 'qi ssfp_emt' - _param_files = ['PD', 'f_b', 'k_bf', 'T1_f', 'T2_f'] - input_spec = eMTSimInputSpec - output_spec = eMTSimOutputSpec - - def _list_outputs(self): - outputs = self.output_spec().get() - outputs['G_file'] = path.abspath(self.inputs.G_file) - outputs['a_file'] = path.abspath(self.inputs.a_file) - outputs['b_file'] = path.abspath(self.inputs.b_file) - return outputs + input_spec = base.SimInputSpec('EMT', + varying=['PD', 'f_b', + 'k_bf', 'T1_f', 'T2_f'], + fixed=['f0', 'B1'], + out_files=['G', 'a', 'b'], + extras={'T2b': traits.Float(desc='T2 of bound pool', argstr='--T2b=%f')}) + output_spec = base.SimOutputSpec('EMT', ['G', 'a', 'b']) ############################ qi_mtsat ############################ -class MTSatInputSpec(QI.InputSpec): +class MTSatInputSpec(base.InputSpec): # Inputs pdw_file = File(exists=True, argstr='%s', mandatory=True, position=-3, desc='Path to PD-weighted data') @@ -333,7 +287,7 @@ class MTSatOutputSpec(TraitedSpec): desc='MTSat delta map', usedefault=True) -class MTSat(QI.FitCommand): +class MTSat(base.FitCommand): """ Runs qi_mtsat @@ -346,7 +300,7 @@ class MTSat(QI.FitCommand): ############################ qi mtr ############################ -class MTRInputSpec(QI.InputSpec): +class MTRInputSpec(base.InputSpec): # Options in_file = File(exists=True, argstr='%s', mandatory=True, position=-1, desc='Input file') @@ -356,7 +310,7 @@ class MTROutputSpec(DynamicTraitedSpec): pass -class MTR(QI.BaseCommand): +class MTR(base.BaseCommand): """ Calculate Magnetization Transfer Ratios """ diff --git a/Python/QUIT/interfaces/perfusion.py b/Python/qipype/interfaces/perfusion.py similarity index 73% rename from Python/QUIT/interfaces/perfusion.py rename to Python/qipype/interfaces/perfusion.py index f7ae61e8..5b01309e 100644 --- a/Python/QUIT/interfaces/perfusion.py +++ b/Python/qipype/interfaces/perfusion.py @@ -15,12 +15,12 @@ from os import path from nipype.interfaces.base import CommandLine, TraitedSpec, File, traits, isdefined -from .. import base as QI +from .. import base ############################ qi_asl ############################ -class ASLInputSpec(QI.FitInputSpec): +class ASLInputSpec(base.FitInputSpec): # Additional Options average = traits.Bool( desc='Average across time-series', argstr='--average') @@ -39,11 +39,7 @@ class ASLInputSpec(QI.FitInputSpec): desc='Apply slice-timing correction (number of PLDs and slices must match)', argstr='--slicetime') -class ASLOutputSpec(TraitedSpec): - CBF_map = File('CASL_CBF.nii.gz', desc="Path to CBF map", usedefault=True) - - -class ASL(QI.FitCommand): +class ASL(base.FitCommand): """ Calculate CBF from CASL data @@ -51,12 +47,12 @@ class ASL(QI.FitCommand): _cmd = 'qi asl' input_spec = ASLInputSpec - output_spec = ASLOutputSpec + output_spec = base.FitOutputSpec('CASL', 'CBF') ############################ qi_ase_oef ############################ -class ASEInputSpec(QI.FitInputSpec): +class ASEInputSpec(base.FitInputSpec): # Additional Options B0 = traits.Float( desc='Field-strength (Tesla), default 3', argstr='--B0=%f') @@ -84,7 +80,7 @@ class ASEOutputSpec(TraitedSpec): 'ASE_dHb.nii.gz', desc='Path to Deoxyhaemoglobin concentration map', usedefault=True) -class ASE(QI.FitCommand): +class ASE(base.FitCommand): """ Calculate the Oxygen Extraction Fraction from Asymmetric Spin Echo data @@ -92,48 +88,37 @@ class ASE(QI.FitCommand): _cmd = 'qi ase_oef' input_spec = ASEInputSpec - output_spec = ASEOutputSpec + output_spec = base.FitOutputSpec( + 'ASE', ['S0', 'dT', 'R2p', 'DBV', 'Tc', 'OEF', 'dHb']) - def _list_outputs(self): - outputs = self.output_spec().get() - if not isdefined(self.inputs.fix_DBV): - outputs.pop('DBV_map', None) - return self._add_prefixes(outputs) +class ASESim(base.SimCommand): + """ + Simulate the Oxygen Extraction Fraction from Asymmetric Spin Echo data with fixed DBV + """ -class ASESimInputSpec(QI.SimInputSpec): - # Additional Options - B0 = traits.Float( - desc='Field-strength (Tesla), default 3', argstr='--B0=%f') - fix_DBV = traits.Float( - desc='Fix Deoxygenated Blood Volume to value (fraction)', argstr='--DBV=%f') - gradz = traits.String( - desc='Field gradient in z/slice-direction for MFG correction', argstr='--gradz=%s') - slice_thickness = traits.Float( - desc='Actual slice-thickness for MFG correction if slice-gap was used', argstr='--slice=%f') + _cmd = 'qi ase_oef' + input_spec = base.SimInputSpec('ASE', + varying=['S0', 'dT', 'R2p'], + extras={'B0': traits.Float(desc='Field-strength (Tesla), default 3', argstr='--B0=%f'), + 'fix_DBV': traits.Float(desc='Fix Deoxygenated Blood Volume to value (fraction)', argstr='--DBV=%f', mandatory=True)}) + output_spec = base.SimOutputSpec('ASE') -class ASESim(QI.SimCommand): +class ASEDBVSim(base.SimCommand): """ - Simulate the Oxygen Extraction Fraction from Asymmetric Spin Echo data - + Simulate OEF with a varying DBV """ - _cmd = 'qi ase_oef' - input_spec = ASESimInputSpec - output_spec = QI.SimOutputSpec - - def __init__(self, **kwargs): - if not 'fix_DBV' in kwargs: - self._param_files = ['S0', 'dT', 'R2p', 'DBV'] - else: - self._param_files = ['S0', 'dT', 'R2p'] - super().__init__(**kwargs) + input_spec = base.SimInputSpec('ASE', + varying=['S0', 'dT', 'R2p', 'DBV'], + extras={'B0': traits.Float(desc='Field-strength (Tesla), default 3', argstr='--B0=%f')}) + output_spec = base.SimOutputSpec('ASE') ############################ qi_zshim ############################ -class ZShimInputSpec(QI.InputBaseSpec): +class ZShimInputSpec(base.InputBaseSpec): in_file = File(argstr='%s', mandatory=True, exists=True, position=-1, desc='Input file to fit polynomial to') zshims = traits.Int(argstr='--zshims=%d', desc='Number of Z-shims') @@ -148,7 +133,7 @@ class ZShimOutputSpec(TraitedSpec): out_file = File(desc="Shimmed Image") -class ZShim(QI.BaseCommand): +class ZShim(base.BaseCommand): """ Combine an EPI image with Z/Y-shimming """ diff --git a/Python/qipype/interfaces/relax.py b/Python/qipype/interfaces/relax.py new file mode 100644 index 00000000..54f5c1e3 --- /dev/null +++ b/Python/qipype/interfaces/relax.py @@ -0,0 +1,449 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- + +""" +Implementation of nipype interfaces for QUIT relaxometry tools +Requires that the QUIT tools are in your your system path +""" + +from os import path, getcwd +from nipype.interfaces.base import TraitedSpec, File, traits +from .. import base + + +############################ qidespot1 ############################ + +class DESPOT1InputSpec(base.FitInputSpec): + # Additional Options + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + algo = traits.String(desc="Choose algorithm (l/w/n)", argstr="--algo=%s") + iterations = traits.Int( + desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') + + +class DESPOT1(base.FitCommand): + """ + Run DESPOT1 analysis with qidespot1 + + Example + ------- + >>> from qipype.interfaces.relax import DESPOT1 + >>> seq = {'SPGR': {'TR':5E-3, 'FA':[5,10]} } + >>> d1 = DESPOT1(sequence=seq, in_file='SPGR.nii.gz') + >>> d1_res = d1.run() + >>> print(d1_res.outputs) + + """ + + _cmd = 'qi despot1' + input_spec = DESPOT1InputSpec + output_spec = base.FitOutputSpec('D1', ['PD', 'T1']) + +############################ qidespot1sim ############################ + + +class DESPOT1Sim(base.SimCommand): + """ + Run DESPOT1 simulation with qidespot1 + + Example with parameter dictionary + ------- + >>> from qipype.interfaces.relax import DESPOT1Sim + >>> seq = {'SPGR': {'TR':5E-3, 'FA':[5,10]}} + >>> d1sim = DESPOT1Sim(sequence=seq, PD_map='PD.nii.gz', T1_map='T1.nii.gz', out_file='SPGR.nii.gz') + >>> d1sim_res = d1.run() + >>> print(d1sim_res.outputs) + """ + + _cmd = 'qi despot1' + input_spec = base.SimInputSpec('D1', ['PD', 'T1'], ['B1']) + output_spec = base.SimOutputSpec('D1') + +############################ qidespot1hifi ############################ + + +class HIFIInputSpec(base.InputSpec): + # Inputs + spgr_file = File(exists=True, argstr='%s', mandatory=True, + position=-2, desc='Path to SPGR data') + + mprage_file = File(exists=True, argstr='%s', mandatory=True, + position=-1, desc='Path to MPRAGE data') + + # Options + residuals = traits.Bool( + desc='Write out residuals for each data-point', argstr='--resids') + clamp_T1 = traits.Float( + desc='Clamp T1 between 0 and value', argstr='--clamp=%f') + + +class HIFI(base.FitCommand): + """ + Calculate T1 & B1 map with the DESPOT1-HIFI method + """ + + _cmd = 'qi despot1hifi' + input_spec = HIFIInputSpec + output_spec = base.FitOutputSpec('HIFI', ['PD', 'T1', 'B1']) + +############################ qidespot1hifisim ############################ + + +class HIFISim(base.SimCommand): + """ + Simulate SPGR/FLASH and MPRAGE images using DESPOT1-HIFI model + """ + + _cmd = 'qi despot1hifi' + sim_files = ['spgr', 'mprage'] + input_spec = base.SimInputSpec( + 'HIFI', ['PD', 'T1', 'B1'], out_files=sim_files) + output_spec = base.SimOutputSpec('HIFI', sim_files) + +############################ qidespot2 ############################ + + +class DESPOT2InputSpec(base.FitInputSpec): + # Inputs + T1_map = File(exists=True, argstr='--T1=%s', mandatory=True, + position=-2, desc='Path to T1 map') + + # Options + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + algo = traits.Enum("LLS", "WLS", "NLS", + desc="Choose algorithm", argstr="--algo=%d") + ellipse = traits.Bool( + desc="Data is ellipse geometric solution", argstr='--gs') + iterations = traits.Int( + desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') + clamp_PD = traits.Float( + desc='Clamp PD between 0 and value', argstr='-f %f') + clamp_T2 = traits.Float( + desc='Clamp T2 between 0 and value', argstr='--clampT1=%f') + + +class DESPOT2(base.FitCommand): + """ + Run DESPOT2 analysis with qidespot2 + """ + + _cmd = 'qi despot2' + input_spec = DESPOT2InputSpec + output_spec = base.FitOutputSpec('D2', ['PD', 'T2']) + +############################ qidespot2sim ############################ + + +class DESPOT2Sim(base.SimCommand): + """ + Run DESPOT2 simulation + """ + + _cmd = 'qi despot2' + input_spec = base.SimInputSpec('D2', + varying=['PD', 'T2'], + fixed=['B1', 'T1'], + extras={'ellipse': traits.Bool( + desc="Data is ellipse geometric solution", argstr='--gs')}) + output_spec = base.SimOutputSpec('D2') + +############################ qidespot2fm ############################ + + +class FMInputSpec(base.FitInputSpec): + # Inputs + T1_map = File(exists=True, argstr='--T1=%s', + mandatory=True, desc='Path to T1 map') + + # Options + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + asym = traits.Bool( + desc="Fit asymmetric (+/-) off-resonance frequency", argstr='--asym') + algo = traits.Enum("LLS", "WLS", "NLS", + desc="Choose algorithm", argstr="--algo=%d") + + +class FM(base.FitCommand): + """ + Run DESPOT2-FM analysis + """ + + _cmd = 'qi despot2fm' + input_spec = FMInputSpec + output_spec = base.FitOutputSpec('FM', ['PD', 'T2', 'f0']) + + +class FMSim(base.SimCommand): + """ + Run DESPOT2-FM simulation + """ + + _cmd = 'qi despot2fm' + _param_files = ['PD', 'T2', 'f0'] + input_spec = base.SimInputSpec('FM', + varying=['PD', 'T2', 'f0'], + fixed=['B1', 'T1'], + extras={'asym': traits.Bool( + desc="Fit asymmetric (+/-) off-resonance frequency", argstr='--asym')}) + output_spec = base.SimOutputSpec('FM') + +############################ qi_jsr ############################ + + +class JSRInputSpec(base.InputSpec): + # Inputs + spgr_file = File(exists=True, argstr='%s', mandatory=True, + position=-2, desc='Path to SPGR data') + + ssfp_file = File(exists=True, argstr='%s', mandatory=True, + position=-1, desc='Path to SSFP data') + + # Options + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + residuals = traits.Bool( + desc='Write out residuals for each data-point', argstr='--resids') + npsi = traits.Int( + desc='Number of psi/off-resonance starts', argstr='--npsi=%d') + + +class JSR(base.FitCommand): + """ + Calculate T1 &T2 map with Joint System Relaxometry + + """ + + _cmd = 'qi jsr' + input_spec = JSRInputSpec + output_spec = base.FitOutputSpec('JSR', ['PD', 'T1', 'T2', 'df0']) + +############################ qimcdespot ############################ +# Status: Everything is there but not tested. Very much unsupported. + + +class QIMCDespotInputSpec(base.InputSpec): + # Inputs + spgr_file = File(exists=True, argstr='%s', mandatory=True, + position=0, desc='SPGR file') + ssfp_file = File(exists=True, argstr='%s', mandatory=True, + position=1, desc='SSFP file') + param_dict = traits.Dict(desc='Parameters as dictionary', position=2, + argstr='', mandatory=True, xor=['param_file']) + param_file = File(desc='Parameters as .json file', position=2, argstr='--json=%s', + xor=['param_dict'], mandatory=True, exists=True) + + # Options + f0_map = File(desc='f0 map (Hertz)', argstr='--f0=%s', exists=True) + B1_map = File(desc='B1 map (ratio)', argstr='--B1=%s', exists=True) + mask = File(desc='Only process voxels within the mask', + argstr='--mask=%s', exists=True) + residuals = traits.Bool( + desc='Write out residuals for each data-point', argstr='--resid') + model = traits.Enum("1", "2", "2nex", "3", "3_f0", "3nex", desc="Select model to fit - 1/2/2nex/3/3_f0/3nex, default 3", + argstr="--model=%d") + scale = traits.Bool( + desc='Normalize signals to mean (a good idea)', argstr='--scale') + algo = traits.Enum( + "S", "G", desc="Select (S)tochastic or (G)aussian Region Contraction", argstr="--algo=%d") + iterations = traits.Int( + desc='Max iterations, default 4', argstr='---its=%d') + field_strength = traits.Float( + desc='Specify field-strength for fitting regions - 3/7/u for user input', argstr='--tesla=%f') + threads = traits.Int( + desc='Use N threads (default=4, 0=hardware limit)', argstr='--threads=%d') + prefix = traits.String( + desc='Add a prefix to output filenames', argstr='--out=%s') + + +class QIMCDespot(base.FitCommand): + """ + Interace for qimcdespot + + Example 1 + ------- + >>> from qipype.interfaces.relax import QiMcDespot + >>> interface = QiMcDespot(prefix='nipype_', param_file='mcdespot_params.json') + """ + + _cmd = 'qi mcdespot' + input_spec = QIMCDespotInputSpec + output_spec = base.FitOutputSpec( + '3C', ['T1_m', 'T2_m', 'T1_ie', 'T2_ie', 'T1_csf', 'T2_csf', 'tau_m', 'f_m', 'f_csf', 'f0', 'B1']) + +############################ qimp2rage ############################ +# Implemented but not tested # + + +class MP2RAGEInputSpec(base.InputSpec): + # Inputs + in_file = File(exists=True, argstr='%s', mandatory=True, + position=0, desc='Path to complex MP-RAGE data') + + # Commonly used options + threads = traits.Int( + desc='Use N threads (default=hardware)', argstr='--threads=%d') + prefix = traits.String( + desc='Add a prefix to output filenames', argstr='--out=%s') + beta = traits.Float(desc='Regularisation paramter', argstr='--beta=%f') + + +class MP2RAGEOutputSpec(TraitedSpec): + # Specify which outputs there are + uni_file = File('MP2_UNI.nii.gz', + desc='The Uniform MP2 contrast image', usedefault=True) + T1_map = File('MP2_T1.nii.gz', desc='T1 Map', usedefault=True) + + +class MP2RAGE(base.FitCommand): + """ + Interface for qimp2rage + + Example 1 + ------- + >>> from qipype.interfaces.relax import QIMP2RAGE + >>> interface = QIMP2RAGE(prefix='nipype_', param_file='spgr_params.json') + + """ + + _cmd = 'qi mp2rage' + input_spec = MP2RAGEInputSpec + output_spec = MP2RAGEOutputSpec + +############################ qimultiecho ############################ + + +class MultiechoInputSpec(base.FitInputSpec): + # Options + algo = traits.String(desc="Choose algorithm (l/a/n)", argstr="--algo=%s") + iterations = traits.Int( + desc='Max iterations for WLLS/NLLS (default 15)', argstr='--its=%d') + thresh_PD = traits.Float( + desc='Only output maps when PD exceeds threshold value', argstr='-t=%f') + clamp_T2 = traits.Float( + desc='Clamp T2 between 0 and value', argstr='-p=%f') + + +class Multiecho(base.FitCommand): + """ + Run a multi-echo fit + """ + + _cmd = 'qi multiecho' + input_spec = MultiechoInputSpec + output_spec = base.FitOutputSpec('ME', ['PD', 'T2']) + + +class MultiechoSim(base.SimCommand): + """ + Run a multi-echo simulation + """ + + _cmd = 'qi multiecho' + input_spec = base.SimInputSpec('ME', ['PD', 'T2']) + output_spec = base.SimOutputSpec('ME') + +############################ qi_mpm_r2s ############################ + + +class MPMR2sInputSpec(base.InputSpec): + # Inputs + pdw_file = File(exists=True, argstr='%s', mandatory=True, + position=-3, desc='Path to PD-weighted data') + + t1w_file = File(exists=True, argstr='%s', mandatory=True, + position=-2, desc='Path to T1-weighted data') + + mtw_file = File(exists=True, argstr='%s', mandatory=True, + position=-1, desc='Path to MT-weighted data') + + # Options + residuals = traits.Bool( + desc='Write out residuals for each data-point', argstr='--resids') + + +class MPMR2s(base.FitCommand): + """ + Runs qi_mpm_r2s + """ + + _cmd = 'qi mpm_r2s' + input_spec = MPMR2sInputSpec + output_spec = base.FitOutputSpec( + 'MPM', ['R2s', 'S0_PDw', 'S0_T1w', 'S0_MTw']) + + +############################ qi_ssfp_elipse ############################ + + +class EllipseInputSpec(base.FitInputSpec): + # Additional Options + algo = traits.String(desc='Choose algorithm (h/d)', argstr='--algo=%s') + + +class Ellipse(base.FitCommand): + """ + Fit an ellipse to SSFP data + """ + + _cmd = 'qi ssfp_ellipse' + input_spec = EllipseInputSpec + output_spec = base.FitOutputSpec( + 'ES', ['G', 'a', 'b', 'theta_0', 'phi_rf']) + + +class EllipseSim(base.SimCommand): + """ + Simulate SSFP data from ellipse parameters + """ + + _cmd = 'qi ssfp_ellipse' + input_spec = base.SimInputSpec('ES', ['G', 'a', 'b', 'theta_0', 'phi_rf']) + output_spec = base.SimOutputSpec('ES') + +############################ qi_planet ############################ + + +class PLANETInputSpec(base.InputSpec): + # Inputs + G_map = File(exists=True, argstr='%s', mandatory=True, + position=-3, desc='Path to G parameter map') + a_map = File(exists=True, argstr='%s', mandatory=True, + position=-2, desc='Path to a parameter map') + b_map = File(exists=True, argstr='%s', mandatory=True, + position=-1, desc='Path to b parameter map') + + # Options + B1_map = File(desc='B1 map (ratio) file', argstr='--B1=%s') + residuals = traits.Bool( + desc='Write out residuals for each data-point', argstr='--resids') + + +class PLANET(base.FitCommand): + """ + Calculate T1/T2 from ellipse parameters + + """ + + _cmd = 'qi planet' + input_spec = PLANETInputSpec + output_spec = base.FitOutputSpec('PLANET', ['PD', 'T1', 'T2']) + + +class PLANETSim(base.SimCommand): + """ + Simulate ellipse parameters from T1/T2 + + """ + + _cmd = 'qi planet' + input_spec = base.SimInputSpec('PLANET', + varying=['PD', 'T1', 'T2'], + fixed=['B1'], + out_files=['G', 'a', 'b']) + output_spec = base.SimOutputSpec('PLANET', ['G', 'a', 'b']) + + def _list_outputs(self): + outputs = self.output_spec().get() + outputs['G_file'] = path.abspath(self.inputs.G_file) + outputs['a_file'] = path.abspath(self.inputs.a_file) + outputs['b_file'] = path.abspath(self.inputs.b_file) + return outputs diff --git a/Python/QUIT/interfaces/rufis.py b/Python/qipype/interfaces/rufis.py similarity index 100% rename from Python/QUIT/interfaces/rufis.py rename to Python/qipype/interfaces/rufis.py diff --git a/Python/QUIT/interfaces/stats.py b/Python/qipype/interfaces/stats.py similarity index 100% rename from Python/QUIT/interfaces/stats.py rename to Python/qipype/interfaces/stats.py diff --git a/Python/QUIT/interfaces/susceptibility.py b/Python/qipype/interfaces/susceptibility.py similarity index 100% rename from Python/QUIT/interfaces/susceptibility.py rename to Python/qipype/interfaces/susceptibility.py diff --git a/Python/QUIT/interfaces/utils.py b/Python/qipype/interfaces/utils.py similarity index 98% rename from Python/QUIT/interfaces/utils.py rename to Python/qipype/interfaces/utils.py index 189a3d35..b04eefea 100644 --- a/Python/QUIT/interfaces/utils.py +++ b/Python/qipype/interfaces/utils.py @@ -94,7 +94,7 @@ class RFProfile(QI.BaseCommand): Example 1 ------- - >>> from QUIT.nipype.utils import RFProfile + >>> from quit.nipype.utils import RFProfile >>> interface = RFProfile() """ @@ -156,7 +156,7 @@ class Affine(CommandLine): Example 1 ------- - >>> from QUIT.nipype.utils import Affine + >>> from quit.nipype.utils import Affine >>> interface = Affine() """ @@ -210,7 +210,7 @@ class Mask(QI.BaseCommand): Example 1 ------- - >>> from QUIT.nipype.utils import RFProfile + >>> from quit.nipype.utils import RFProfile >>> interface = RFProfile() """ @@ -313,7 +313,7 @@ class Complex(QI.BaseCommand): Example 1 ------- - >>> from QUIT.nipype.utils import Complex + >>> from quit.nipype.utils import Complex >>> interface = Complex() """ diff --git a/Python/QUIT/sims.py b/Python/qipype/sims.py similarity index 95% rename from Python/QUIT/sims.py rename to Python/qipype/sims.py index 6574a013..93f03769 100644 --- a/Python/QUIT/sims.py +++ b/Python/qipype/sims.py @@ -34,6 +34,7 @@ def init_brainweb(): url='https://brainweb.bic.mni.mcgill.ca/cgi/brainweb1', params=params) b1_file = open('rf20_C.mnc', 'wb') b1_file.write(response.content) + print('Finished') def make_phantom(parameters, subsamp=1, B1=True, f0=False): @@ -55,6 +56,7 @@ def make_phantom(parameters, subsamp=1, B1=True, f0=False): class_data[::subsamp, ::subsamp, ::subsamp], np.array(vals)).astype('float32') img = nib.nifti1.Nifti1Image(data, affine=classes.affine) nib.save(img, key + '.nii.gz') + print('Wrote', key + '.nii.gz') if B1: b1_img = nib.load('rf20_C.mnc') @@ -62,9 +64,11 @@ def make_phantom(parameters, subsamp=1, B1=True, f0=False): 'float32')[::subsamp, ::subsamp, ::subsamp] img = nib.nifti1.Nifti1Image(b1_data, affine=b1_img.affine) nib.save(img, 'B1.nii.gz') + print('Wrote B1.nii.gz') if f0: shape = [int(np.ceil(x / subsamp)) for x in class_data.shape] f0data = np.tile( np.linspace(-60, 60, shape[2]), [shape[0], shape[1], 1]) img = nib.nifti1.Nifti1Image(f0data, affine=classes.affine) nib.save(img, 'f0.nii.gz') + print('Wrote f0.nii.gz') diff --git a/Python/QUIT/workflows/cest.py b/Python/qipype/workflows/cest.py similarity index 98% rename from Python/QUIT/workflows/cest.py rename to Python/qipype/workflows/cest.py index 9b2cbe2e..04686f44 100644 --- a/Python/QUIT/workflows/cest.py +++ b/Python/qipype/workflows/cest.py @@ -1,8 +1,8 @@ import numpy as np from nipype import Workflow, Node, MapNode, IdentityInterface from nipype.interfaces.fsl import maths, BET, MCFLIRT, FLIRT, ApplyXFM, ImageMaths, Smooth -from QUIT.interfaces.mt import Lorentzian, LorentzianSim, ZSpec -from QUIT.interfaces.utils import PCA, Select +from quit.interfaces.mt import Lorentzian, LorentzianSim, ZSpec +from quit.interfaces.utils import PCA, Select def prep(zfrqs, dummies=0, pca_retain=0, name='CEST_prep'): @@ -105,7 +105,7 @@ def cert(zfrqs): name='cert_subtract') amide_index = (np.abs(zfrqs - 3.5)).argmin() amide = Node(Select(volumes=[amide_index], out_file='amide.nii.gz'), - name='select_amide') + name='select_amide') cert = Workflow(name='CERT') cert.connect([(inputnode, cert_sub, [('cert_360', 'in_file'), ('cert_180', 'in_file2')]), diff --git a/Python/QUIT/workflows/mpm.py b/Python/qipype/workflows/mpm.py similarity index 95% rename from Python/QUIT/workflows/mpm.py rename to Python/qipype/workflows/mpm.py index 0b9f8c58..9f32a32e 100644 --- a/Python/QUIT/workflows/mpm.py +++ b/Python/qipype/workflows/mpm.py @@ -2,8 +2,8 @@ from nipype import Workflow, Node, MapNode, IdentityInterface from nipype.interfaces.fsl import maths, BET, MCFLIRT, ExtractROI, FLIRT, Merge, ApplyXFM, ImageMaths, BinaryMaths, ConvertXFM import nipype.interfaces.utility as util -from QUIT.interfaces.relax import MPMR2s -from QUIT.interfaces.mt import MTSat +from quit.interfaces.relax import MPMR2s +from quit.interfaces.mt import MTSat from .interfaces import ApplyXfm4D diff --git a/Python/QUIT/workflows/prep.py b/Python/qipype/workflows/prep.py similarity index 98% rename from Python/QUIT/workflows/prep.py rename to Python/qipype/workflows/prep.py index 26c87965..e9a35e97 100644 --- a/Python/QUIT/workflows/prep.py +++ b/Python/qipype/workflows/prep.py @@ -3,8 +3,8 @@ from nipype.interfaces.fsl import maths, BET, MCFLIRT, ExtractROI, FLIRT, Merge, ApplyXFM, ImageMaths, BinaryMaths, ConvertXFM from nipype.interfaces.ants import Registration, ApplyTransforms import nipype.interfaces.utility as util -from QUIT.interfaces.utils import Mask, RFProfile, Complex, Filter, CoilCombine -from QUIT.interfaces.fsl import ApplyXfm4D +from quit.interfaces.utils import Mask, RFProfile, Complex, Filter, CoilCombine +from quit.interfaces.fsl import ApplyXfm4D def COMPOSER(verbose=False, is_bruker=False): diff --git a/Python/setup.py b/Python/setup.py index 8d8282eb..2af4164d 100644 --- a/Python/setup.py +++ b/Python/setup.py @@ -12,20 +12,22 @@ long_description = f.read() setup( - name='QUIT', - version='0.1.0', + name='qipype', + version='1.0', description='nipype interfaces to QUantitative Imaging Tools', + long_description=long_description, + long_description_content_type="text/markdown", url='https://github.com/spinicist/quit', author='Tobias Wood', author_email='tobias@spinicist.org.uk', - py_modules=['QUIT'], - install_requires=['nipype', - 'nibabel'], + py_modules=['qipype'], + install_requires=['nipype>=1.2.3', + 'nibabel>=2.5.1'], + python_requires='>=3', license='MPL', - classifiers=['Development Status :: 4 - Beta', - 'Topic :: Scientific/Engineering :: Visualization', + classifiers=['Topic :: Scientific/Engineering :: Physics', 'Programming Language :: Python :: 3', ], - keywords='neuroimaging nifti', + keywords='neuroimaging mri', packages=find_packages() ) diff --git a/README.md b/README.md index c87ff980..99002346 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ![Logo](Docs/logo.png) -![Build Status](https://github.com/actions/hello-world/workflows/Build/badge.svg) +![Build](https://github.com/spinicist/QUIT/workflows/Build/badge.svg) [![DOI](http://joss.theoj.org/papers/10.21105/joss.00656/status.svg)](https://doi.org/10.21105/joss.00656) [![DOI](https://zenodo.org/badge/37066948.svg)](https://zenodo.org/badge/latestdoi/37066948) @@ -12,9 +12,6 @@ file, You can obtain one at http://mozilla.org/MPL/2.0/. If you find the tools useful the author would love to hear from you. -This is the updated version of QUIT based on ITK http://www.itk.org -The previous version is here http://github.com/spinicist/old_quit - ## Brief Description A collection of programs for processing quantitative MRI data, originally DESPOT @@ -22,9 +19,8 @@ but now a wider variety of techniques. ## Thanks -Many thanks to Samuel Hurley for many, many suggestions. -Thanks to Anna Combes, Amy McDowell and Sjoerd Vos for finding bugs in early -versions. +Many thanks to Samuel Hurley, Erika Raven, Anna Combes, Amy McDowell and +Sjoerd Vos. ## Documentation @@ -33,36 +29,46 @@ in .rst format). ## Installation -Pre-compiled binaries are provided for Linux and Mac OS X in a .tar.gz archive -from http://github.com/spinicist/QUIT/releases. +### Executable + +Pre-compiled executables are provided for Linux and Mac OS X in a .tar.gz +archive from http://github.com/spinicist/QUIT/releases. Download the archive and +extract it with `tar -xzf qi-platform.tar.gz`. Then, move the resulting `qi` +file to somewhere on your `$PATH`, for instance `/usr/local/bin`. That's it. -Download the correct archive for your platform, untar it and then ensure that -the binaries can be found via your `PATH` environment variable. The Linux -binaries are built with Ubuntu 14.04 with GCC 6. If you need to run on an older -version of Linux with a previous version of `glibc` then you will need to -compile from source. +- MacOS Catalina or higher users should use `curl` to download the binary, i.e. + type `curl -L https://github.com/spinicist/QUIT/releases/download/v3.0/qi-macos.tar.gz` + This is because Safari now sets the quarantine attribute of all downloads, + which prevents them being run as the binary is unsigned. It is possible to + remove the quarantine flag with `xattr`, but downloading with `curl` is more + straightforward. +- The Linux executable is compiled on Ubuntu 16.04 with GLIBC version 2.3 and a + statically linked libc++. This means it will hopefully run on most modern + Linux distributions. Let me know if it doesn't. For instructions on how to compile from source, please see the developer page in the documentation. +### Python + +QUIT now comes with `nipype` wrappers. Install with `pip install QUIT`. + ## Usage -Example scripts for mcDESPOT processing are provided in the installation -archive. These print usage instructions if you call them with no arguments. -They take a set of filenames as input, and you will need to modify the scripts -with your particular flip-angles and TRs. +QUIT comes as a single executable file with multiple commands, similar to `git` +or `bart`. Type `qi` to see a list of all the available commands. The majority +of commands require an input 4D Nifti file containing the image data to fit, +and a `.json` file containing the sequence parameters (e.g. `TR`). See +examples for each command at http://quit.readthedocs.io. -Each product has some basic usage instructions that will be printed with either -the -h or --help options, e.g. `qimcdespot -h`. The majority of programs take -the input filenames as arguments on the command-line, and in addition will read -an input text file from `stdin`. For further details, see the -[documentation](https://spinicist.github.io/QUIT), which is also available in -the `Docs` folder. +For the `nipype` wrappers, there are example iPython notebooks available in this +directory https://github.com/spinicist/QUIT/blob/master/Python/. You can +download the one for T1 mapping here https://github.com/spinicist/QUIT/raw/master/Python/T1_mapping.ipynb. +Alternatively, look at the example workflows for CEST and MPM to see how to +build a complete pipeline. ## Getting Help If you can't find an answer to a problem in the documentation or help strings, -you can open an [issue](https://github.com/spinicist/QUIT/issues), post a -question on [Neurostars](https://neurostars.org) or find the main developer on -e-mail (tobias.wood@kcl.ac.uk) and Twitter -([@spinicist](https://twitter.com/spinicist)). +you can open an [issue](https://github.com/spinicist/QUIT/issues), or find the +developer on Twitter ([@spinicist](https://twitter.com/spinicist)). diff --git a/Source/Core/Args.h b/Source/Core/Args.h index 32c1312d..1bd66b71 100644 --- a/Source/Core/Args.h +++ b/Source/Core/Args.h @@ -37,7 +37,7 @@ template T CheckValue(args::ValueFlag &v) { if (v) { return v.Get(); } else { - QI::Fail("{} was not specified. Use --help to see usage.", v.Name()); + QI::Fail("{} was not specified but is required. Use --help to see usage.", v.Name()); } } diff --git a/Source/Core/SimulateModel.h b/Source/Core/SimulateModel.h index f33fe389..fdea9437 100644 --- a/Source/Core/SimulateModel.h +++ b/Source/Core/SimulateModel.h @@ -30,7 +30,7 @@ void SimulateModel(json & json, auto simulator = QI::ModelSimFilter::New(model, verbose, subRegion); simulator->SetNoise(noise); for (auto i = 0; i < Model::NV; i++) { - const std::string vname = model.varying_names[i]; + const std::string vname = fmt::format("{}_map", model.varying_names[i]); const std::string vfile = json.at(vname).get(); QI::Log(verbose, "Reading {} from file: {}", vname, vfile); simulator->SetVarying(i, QI::ReadImage(vfile, false)); diff --git a/Source/MT/qi_qmt.cpp b/Source/MT/qi_qmt.cpp index 7fa5b9c8..d2c4bab1 100644 --- a/Source/MT/qi_qmt.cpp +++ b/Source/MT/qi_qmt.cpp @@ -123,9 +123,9 @@ using RamaniFitFunction = QI::ScaledAutoDiffFit; // Main //****************************************************************************** int qmt_main(args::Subparser &parser) { - args::Positional T1(parser, "T1", "T1 map (seconds) file"); args::Positional mtsat_path(parser, "MTSAT FILE", "Path to MT-Sat data"); QI_COMMON_ARGS; + args::ValueFlag T1(parser, "T1", "T1 map (seconds) file ** REQUIRED **", {"T1"}); args::ValueFlag f0(parser, "f0", "f0 map (Hz) file", {'f', "f0"}); args::ValueFlag B1(parser, "B1", "B1 map (ratio) file", {'b', "B1"}); args::ValueFlag lineshape_arg( @@ -162,7 +162,7 @@ int qmt_main(args::Subparser &parser) { if (simulate) { QI::SimulateModel(input, model, - {f0.Get(), B1.Get(), T1.Get()}, + {f0.Get(), B1.Get(), QI::CheckValue(T1)}, {mtsat_path.Get()}, mask.Get(), verbose, @@ -173,7 +173,8 @@ int qmt_main(args::Subparser &parser) { auto fit_filter = QI::ModelFitFilter::New( &fit, verbose, covar, resids, subregion.Get()); - fit_filter->ReadInputs({mtsat_path.Get()}, {f0.Get(), B1.Get(), T1.Get()}, mask.Get()); + fit_filter->ReadInputs( + {mtsat_path.Get()}, {f0.Get(), B1.Get(), QI::CheckValue(T1)}, mask.Get()); fit_filter->Update(); fit_filter->WriteOutputs(prefix.Get() + "QMT_"); QI::Log(verbose, "Finished."); diff --git a/Source/Perfusion/qi_ase_oef.cpp b/Source/Perfusion/qi_ase_oef.cpp index f45cee22..35d77379 100644 --- a/Source/Perfusion/qi_ase_oef.cpp +++ b/Source/Perfusion/qi_ase_oef.cpp @@ -150,11 +150,6 @@ int ase_oef_main(args::Subparser &parser) { QI_COMMON_ARGS; args::ValueFlag B0(parser, "B0", "Field-strength (Tesla), default 3", {'B', "B0"}, 3.0); args::ValueFlag DBV(parser, "DBV", "Fix DBV and only fit R2'", {'d', "DBV"}, 0.0); - args::ValueFlag slice_arg( - parser, - "SLICE THICKNESS", - "Slice-thickness for MFG calculation (useful if there was a slice gap)", - {'s', "slice"}); parser.Parse(); json input = json_file ? QI::ReadJSON(json_file.Get()) : QI::ReadJSON(std::cin); diff --git a/Source/Relaxometry/qidespot2.cpp b/Source/Relaxometry/qidespot2.cpp index 8e608b4e..fc89a2ad 100644 --- a/Source/Relaxometry/qidespot2.cpp +++ b/Source/Relaxometry/qidespot2.cpp @@ -229,9 +229,9 @@ struct DESPOT2NLLS : DESPOT2Fit { // Main //****************************************************************************** int despot2_main(args::Subparser &parser) { - args::Positional t1_path(parser, "T1 MAP", "Path to T1 map"); args::Positional ssfp_path(parser, "SSFP FILE", "Path to SSFP data"); QI_COMMON_ARGS; + args::ValueFlag t1_path(parser, "T1 MAP", "Path to T1 map **REQUIRED**", {"T1"}); args::ValueFlag B1(parser, "B1", "B1 map (ratio) file", {'b', "B1"}); args::ValueFlag algorithm(parser, "ALGO", "Choose algorithm (l/w/n)", {'a', "algo"}, 'l'); args::Flag gs_arg( @@ -249,7 +249,7 @@ int despot2_main(args::Subparser &parser) { model.elliptical = true; QI::SimulateModel(input, model, - {QI::CheckPos(t1_path), B1.Get()}, + {QI::CheckValue(t1_path), B1.Get()}, {QI::CheckPos(ssfp_path)}, mask.Get(), verbose, @@ -276,7 +276,7 @@ int despot2_main(args::Subparser &parser) { d2->model.elliptical = true; } auto fit = QI::ModelFitFilter::New(d2, verbose, covar, resids, subregion.Get()); - fit->ReadInputs({QI::CheckPos(ssfp_path)}, {QI::CheckPos(t1_path), B1.Get()}, mask.Get()); + fit->ReadInputs({QI::CheckPos(ssfp_path)}, {QI::CheckValue(t1_path), B1.Get()}, mask.Get()); fit->Update(); fit->WriteOutputs(prefix.Get() + "D2_"); QI::Log(verbose, "Finished."); diff --git a/Source/Relaxometry/qidespot2fm.cpp b/Source/Relaxometry/qidespot2fm.cpp index a6a5a9e4..fdf1ec37 100644 --- a/Source/Relaxometry/qidespot2fm.cpp +++ b/Source/Relaxometry/qidespot2fm.cpp @@ -160,9 +160,9 @@ struct FMNLLS : FMFit { // Main //****************************************************************************** int despot2fm_main(args::Subparser &parser) { - args::Positional t1_path(parser, "T1_MAP", "Input T1 map"); args::Positional ssfp_path(parser, "SSFP_FILE", "Input SSFP file"); QI_COMMON_ARGS; + args::ValueFlag t1_path(parser, "T1_MAP", "Input T1 map ** REQUIRED **", {"T1"}); args::ValueFlag B1(parser, "B1", "B1 map (ratio) file", {'b', "B1"}); args::ValueFlag its( parser, "ITERS", "Max iterations for NLLS (default 75)", {'i', "its"}, 75); @@ -176,7 +176,7 @@ int despot2fm_main(args::Subparser &parser) { if (simulate) { QI::SimulateModel(input, model, - {QI::CheckPos(t1_path), B1.Get()}, + {QI::CheckValue(t1_path), B1.Get()}, {QI::CheckPos(ssfp_path)}, mask.Get(), verbose, @@ -189,7 +189,7 @@ int despot2fm_main(args::Subparser &parser) { auto fit_filter = QI::ModelFitFilter::New(&fm, verbose, covar, resids, subregion.Get()); fit_filter->ReadInputs( - {QI::CheckPos(ssfp_path)}, {QI::CheckPos(t1_path), B1.Get()}, mask.Get()); + {QI::CheckPos(ssfp_path)}, {QI::CheckValue(t1_path), B1.Get()}, mask.Get()); fit_filter->Update(); fit_filter->WriteOutputs(prefix.Get() + "FM_"); QI::Log(verbose, "Finished.");