From f4a5eaec8f605cb253417ce2708c8c16028743ea Mon Sep 17 00:00:00 2001 From: Helge <47348963+HJZollner@users.noreply.github.com> Date: Fri, 8 Mar 2024 15:47:48 -0500 Subject: [PATCH] Add LCModel support for GABA-edited MEGA-PRESS Full implementation of LCModel wrapper for GABA-edited MEGA-PRESS with the optimal settings and MM models as in Osprey --- GUI/osp_updateFitWindow.m | 28 +- exampledata/sdat/MEGA/jobSDAT_MEGA.m | 8 +- exampledata/sdat/MEGA/jobSDAT_MEGA_LCModel.m | 348 +++++++++++++ fit/code/LCMwrapper/LCModelWrapper.m | 3 +- fit/code/LCMwrapper/readLCMFitParams.m | 8 +- fit/code/osp_fitInitialise.m | 98 +++- fit/code/osp_fitMEGA.m | 463 +++++++++--------- fit/code/osp_fit_Quality.m | 12 +- .../FID-A/inputOutput/io_loadspec_niimrs.m | 20 +- libraries/FID-A/processingTools/op_rmempty.m | 4 + load/OspreyLoad.m | 16 + overview/OspreyOverview.m | 17 +- plot/osp_plotLoad.m | 9 +- process/osp_onOffClassifyMEGA.m | 4 +- process/osp_saveLCM.m | 10 +- quantify/OspreyQuantify.m | 81 +-- 16 files changed, 823 insertions(+), 306 deletions(-) create mode 100644 exampledata/sdat/MEGA/jobSDAT_MEGA_LCModel.m diff --git a/GUI/osp_updateFitWindow.m b/GUI/osp_updateFitWindow.m index 19519f72..91ecfa10 100644 --- a/GUI/osp_updateFitWindow.m +++ b/GUI/osp_updateFitWindow.m @@ -70,7 +70,7 @@ function osp_updateFitWindow(gui) % No info panel string for the water fit range waterFitRangeString = ''; % Where are the metabolite names stored? - basisSetNames = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected}.name; + basisSetNames = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected,subspectrum}.name; % Smaller fonts for the results resultsFontSize = 9; case 'Osprey' @@ -99,10 +99,10 @@ function osp_updateFitWindow(gui) switch MRSCont.opts.fit.method case 'LCModel' if strcmp(gui.fit.Names{gui.fit.Selected}, 'ref') || strcmp(gui.fit.Names{gui.fit.Selected}, 'w') - RawAmpl = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected}.h2oarea .* MRSCont.fit.scale{1,gui.controls.Selected}; + RawAmpl = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected,subspectrum}.h2oarea .* MRSCont.fit.scale{1,gui.controls.Selected}; else - RawAmpl = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected}.ampl .* MRSCont.fit.scale{1,gui.controls.Selected}; - CRLB = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected}.CRLB; + RawAmpl = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected,subspectrum}.ampl .* MRSCont.fit.scale{1,gui.controls.Selected}; + CRLB = MRSCont.fit.results.(gui.fit.Style).fitParams{1,gui.controls.Selected,subspectrum}.CRLB; end case 'Osprey' RawAmpl = MRSCont.fit.results.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.ampl .* MRSCont.fit.scale{1,gui.controls.Selected}; @@ -125,19 +125,19 @@ function osp_updateFitWindow(gui) switch MRSCont.opts.fit.method case 'LCModel' if strcmp(gui.fit.Names{gui.fit.Selected}, 'ref') || strcmp(gui.fit.Names{gui.fit.Selected}, 'w') - RawAmpl = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected}.h2oarea .* MRSCont.fit.scale{gui.controls.Selected}; + RawAmpl = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.h2oarea .* MRSCont.fit.scale{gui.controls.Selected}; else - RawAmpl = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; - CRLB = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected}.CRLB; + RawAmpl = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; + CRLB = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.CRLB; end case 'Osprey' RawAmpl = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; end - ph0 = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.ph0; - ph1 = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.ph1; + ph0 = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.ph0; + ph1 = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.ph1; if ~strcmp(gui.fit.Names{gui.fit.Selected}, 'ref') && ~strcmp(gui.fit.Names{gui.fit.Selected}, 'w') - refShift = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.refShift; - refFWHM = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.refFWHM; + refShift = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.refShift; + refFWHM = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.refFWHM; switch MRSCont.opts.fit.method case 'Osprey' iniph0 = MRSCont.fit.results{1,gui.controls.act_x}.(gui.fit.Style).fitParams{basis,gui.controls.Selected}.prelimParams.ph0; @@ -151,10 +151,10 @@ function osp_updateFitWindow(gui) switch MRSCont.opts.fit.method case 'LCModel' if strcmp(gui.fit.Names{gui.fit.Selected}, 'ref') || strcmp(gui.fit.Names{gui.fit.Selected}, 'w') - RawAmpl = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected}.h2oarea .* MRSCont.fit.scale{gui.controls.Selected}; + RawAmpl = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.h2oarea .* MRSCont.fit.scale{gui.controls.Selected}; else - RawAmpl = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; - CRLB = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected}.CRLB; + RawAmpl = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; + CRLB = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{gui.controls.Selected,subspectrum}.CRLB; end case 'Osprey' RawAmpl = MRSCont.fit.results{gui.controls.act_x,gui.controls.act_y}.(gui.fit.Style).fitParams{basis,gui.controls.Selected,subspectrum}.ampl .* MRSCont.fit.scale{gui.controls.Selected}; diff --git a/exampledata/sdat/MEGA/jobSDAT_MEGA.m b/exampledata/sdat/MEGA/jobSDAT_MEGA.m index adf28f21..2d41a247 100644 --- a/exampledata/sdat/MEGA/jobSDAT_MEGA.m +++ b/exampledata/sdat/MEGA/jobSDAT_MEGA.m @@ -111,7 +111,7 @@ opts.SubSpecAlignment.mets = 'L2Norm'; % OPTIONS: - 'L2Norm' (default) % - 'L1Norm' % - 'none' - +opts.UnstableWater = 0; %Perform eddy-current correction on the metabolite data (raw) or metabolite %-nulled data (mm). This can either be done similar for all data sets by %supplying a single value or specified for each dataset individually by supplying @@ -126,10 +126,10 @@ % - [] array % Save LCModel-exportable files for each spectrum? -opts.saveLCM = 0; % OPTIONS: - 0 (no, default) +opts.saveLCM = 1; % OPTIONS: - 0 (no, default) % - 1 (yes) % Save jMRUI-exportable files for each spectrum? -opts.savejMRUI = 0; % OPTIONS: - 0 (no, default) +opts.savejMRUI = 1; % OPTIONS: - 0 (no, default) % - 1 (yes) % Save processed spectra in vendor-specific format (SDAT/SPAR, RDA, P)? @@ -184,7 +184,7 @@ % - 1 (yes, default) % How do you want to model the co-edited macromolecules at 3 ppm for GABA-edited MRS? -opts.fit.coMM3 = 'freeGauss'; % OPTIONS: - {'3to2MM'} (default) +opts.fit.coMM3 = '3to2MM'; % OPTIONS: - {'3to2MM'} (default) % - {'3to2MMsoft'} % - {'1to1GABA'} % - {'1to1GABAsoft'} diff --git a/exampledata/sdat/MEGA/jobSDAT_MEGA_LCModel.m b/exampledata/sdat/MEGA/jobSDAT_MEGA_LCModel.m new file mode 100644 index 00000000..ed1b6dd9 --- /dev/null +++ b/exampledata/sdat/MEGA/jobSDAT_MEGA_LCModel.m @@ -0,0 +1,348 @@ +%% jobSDAT.m +% This function describes an Osprey job defined in a MATLAB script. +% +% A valid Osprey job contains four distinct classes of items: +% 1. basic information on the MRS sequence used +% 2. several settings for data handling and modeling +% 3. a list of MRS (and, optionally, structural imaging) data files +% to be loaded +% 4. an output folder to store the results and exported files +% +% The list of MRS and structural imaging files is provided in the form of +% cell arrays. They can simply be provided explicitly, or from a more +% complex script that automatically determines file names from a given +% folder structure. +% +% Osprey distinguishes between four sets of data: +% - metabolite (water-suppressed) data +% (MANDATORY) +% Defined in cell array "files" +% - water reference data acquired with the SAME sequence as the +% metabolite data, just without water suppression RF pulses. This +% data is used to determine complex coil combination +% coefficients, and perform eddy current correction. +% (OPTIONAL) +% Defined in cell array "files_ref" +% - additional water data used for water-scaled quantification, +% usually from short-TE acquisitions due to reduced T2-weighting +% (OPTIONAL) +% Defined in cell array "files_w" +% - Structural image data used for co-registration and tissue class +% segmentation (usually a T1 MPRAGE). These files need to be +% provided in the NIfTI format (*.nii) or, for GE data, as a +% folder containing DICOM Files (*.dcm). +% (OPTIONAL) +% Defined in cell array "files_nii" +% - External segmentation results. These files need to be +% provided in the NIfTI format (*.nii or *.nii.gz). +% (OPTIONAL) +% Defined in cell array "files_seg" with 1 x 3 cell for each +% subject or 1 x 1 cell if a single 4D NIfTI is supplied. +% +% Files in the formats +% - .7 (GE) +% - .SDAT, .DATA/.LIST, .RAW/.SIN/.LAB (Philips) +% - .DAT (Siemens) +% - .nii, .nii.gz (NIfTI-MRS) +% usually contain all of the acquired data in a single file per scan. GE +% systems store water reference data in the same .7 file, so there is no +% need to specify it separately under files_ref. +% +% Files in the formats +% - .DCM (any) +% - .IMA, .RDA (Siemens) +% may contain separate files for each average. Instead of providing +% individual file names, please specify folders. Metabolite data, water +% reference data, and water data need to be located in separate folders. +% +% In the example script at hand the MATLAB functions strrep and which are +% used to generate a relative path, which allows you to run the examples +% on your machine directly. To set up your own Osprey job supply the +% specific locations as described above. +% +% AUTHOR: +% Dr. Georg Oeltzschner (Johns Hopkins University, 2019-07-15) +% goeltzs1@jhmi.edu +% +% HISTORY: +% 2019-07-15: First version of the code. + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% 1. SPECIFY SEQUENCE INFORMATION %%% + +% Specify sequence type +seqType = 'MEGA'; % OPTIONS: - 'unedited' (default) + % - 'MEGA' + % - 'HERMES' + % - 'HERCULES' + +% Specify editing targets +editTarget = {'GABA'}; % OPTIONS: - {'none'} (default if 'unedited') + % - {'GABA'}, {'GSH'}, {'Lac'}, {'PE322'}, {'PE398'} (for 'MEGA') + % - {'GABA', 'GSH'}, {'GABA', 'Lac'}, {'NAA', 'NAAG'} (for 'HERMES'and 'HERCULES') + + % Specify data scenario +dataScenario = 'invivo'; % OPTIONS: - 'invivo' (default) + % - 'phantom' + % - 'PRIAM' + % - 'MRSI' +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% 2. SPECIFY DATA HANDLING AND MODELING OPTIONS %%% +% Which spectral registration method should be used? Robust spectral +% registration is default, a frequency restricted spectral registration +% method is also availaible and is linked to the fit range. +opts.SpecReg = 'RobSpecReg'; % OPTIONS: - 'RobSpecReg' (default) Spectral aligment with Water/Lipid removal, using simialrity meric, and weighted averaging + % - 'ProbSpecReg' Probabilistic spectral aligment to median target and weighted averaging + % - 'RestrSpecReg' Frequency restricted (fit range) spectral aligment, using simialrity meric, and weighted averaging + % - 'none' + +% Which algorithm do you want to align the sub spectra? L2 norm +% optimization is the default. This is only used for edited MRS! +% Which algorithm do you want to align the sub spectra? L2 norm +% optimization is the default. This is only used for edited MRS! +%Perform correction on the metabolite data (raw) or metabolite +%-nulled data (mm). +opts.SubSpecAlignment.mets = 'L2Norm'; % OPTIONS: - 'L2Norm' (default) + % - 'L1Norm' + % - 'none' + +%Perform eddy-current correction on the metabolite data (raw) or metabolite +%-nulled data (mm). This can either be done similar for all data sets by +%supplying a single value or specified for each dataset individually by supplying +% multiple entries (number has to match the number of datasets) e.g. to perform ECC +% for the second dataset only: +% opts.ECC.raw = [0 1]; +% opts.ECC.mm = [0 1]; + + +opts.ECC.raw = 1; % OPTIONS: - '1' (default) +opts.ECC.mm = 1; % - '0' (no) + % - [] array + +% Save LCModel-exportable files for each spectrum? +opts.saveLCM = 1; % OPTIONS: - 0 (no, default) + % - 1 (yes) +% Save jMRUI-exportable files for each spectrum? +opts.savejMRUI = 0; % OPTIONS: - 0 (no, default) + % - 1 (yes) + +% Save processed spectra in vendor-specific format (SDAT/SPAR, RDA, P)? +opts.saveVendor = 0; % OPTIONS: - 0 (no, default) + % - 1 (yes) + +% Save processed spectra in NIfTI-MRS format? +opts.saveNII = 1; % OPTIONS: - 0 (no) + % - 1 (yes, default) + +% Save PDF output for all Osprey modules and subjects? +opts.savePDF = 0; % OPTIONS: - 0 (no, default) + % - 1 (yes) + +% Save mat file of the compiled fitting parameters? +opts.exportParams.flag = 0; % Options: - 0 (no, default) + % - 1 (yes) +opts.exportParams.path = ''; % Replace with string for the path + % to the save directory + +% Choose the fitting algorithm +opts.fit.method = 'LCModel'; % OPTIONS: - 'Osprey' (default) + +% Select the metabolites to be included in the basis set as a cell array, +% with entries separates by commas. +% With default Osprey basis sets, you can select the following metabolites: +% Ala, Asc, Asp, bHB, bHG, Cit, Cr, Cystat, CrCH2, EtOH, GABA, GPC, GSH, Glc, Gln, +% Glu, Gly, H2O, mI, Lac, NAA, NAAG, PCh, PCr, PE, Phenyl, sI, Ser, +% Tau, Tyros, MM09, MM12, MM14, MM17, MM20, Lip09, Lip13, Lip20. +% If you enter 'default', the basis set will include all of the above +% except for Ala, bHB, bHG, Cit, Cystat, EtOH, Glc, Gly, Phenyl, Ser, and Tyros. +opts.fit.includeMetabs = {'default'}; % OPTIONS: - {'default'} + % - {'full'} + % - {custom} + +% Choose the fitting style for difference-edited datasets (MEGA, HERMES, HERCULES) +% (only available for the Osprey fitting method) +opts.fit.style = 'Separate'; % OPTIONS: - 'Concatenated' (default) - will fit DIFF and SUM simultaneously) + % - 'Separate' - will fit DIFF and OFF separately + +% Determine fitting range (in ppm) for the metabolite and water spectra +opts.fit.range = [0.5 4]; % [ppm] Default: [0.2 4.2] +opts.fit.rangeWater = [2.0 7.4]; % [ppm] Default: [2.0 7.4] +opts.fit.GAP.A = []; +opts.fit.GAP.diff1 = [1.2 1.95]; % [ppm] Default: [1.2 1.95] + +% Determine the baseline knot spacing (in ppm) for the metabolite spectra +opts.fit.bLineKnotSpace = Inf; % [ppm] Default: Inf. + % Inf sets option to nobaseline in LCModel,i.e., nobase = T + +% Add macromolecule and lipid basis functions to the fit? +opts.fit.fitMM = 1; % OPTIONS: - 0 (no) + % - 1 (yes, default) + +% How do you want to model the co-edited macromolecules at 3 ppm for GABA-edited MRS? +opts.fit.coMM3 = '3to2MM'; % OPTIONS: - {'3to2MM'} (default) + % - {'3to2MMsoft'} + % - {'1to1GABA'} + % - {'1to1GABAsoft'} + % - {'fixedGauss'} + % - {'none'} + +opts.fit.FWHMcoMM3 = 14; + +% Optional: In case the automatic basisset picker is not working you can manually +% select the path to the basis set in the osprey/fit/basis, i.e.: +% opts.fit.basisSetFile = 'osprey/fit/basis/3T/philips/mega/press/gaba68/basis_philips_megapress_gaba68.mat'; + +% Optional: Deface the strucutral images in the Coreg/Seg figures for HIPAA +% compliance +opts.img.deface = 0; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% 3. SPECIFY MRS DATA AND STRUCTURAL IMAGING FILES %% +% When using single-average Siemens RDA or DICOM files, specify their +% folders instead of single files! + +% Clear existing files +clear files files_ref files_w files_nii files_mm + +% Data folder in BIDS format +% The filparts(which()) comment is needed to find the data on your machine. If you set +% up the jobFile for your own data you can set a direct path to your data +% folder e.g., data_folder = /Volumes/MyProject/data/' + +data_folder = fileparts(which(fullfile('exampledata','sdat','MEGA','jobSDAT_MEGA.m'))); + +% The following lines perform an automated set-up of the jobFile which +% takes advatage of the BIDS foramt. If you are not using BIDS (highly +% recommended) you can look at the definitions below the loop to see how to +% set up direct path links to your data. + +subs = dir(data_folder); +subs(1:2) = []; +subs = subs([subs.isdir]); +subs = subs(contains({subs.name},'sub')); +counter = 1; + +for kk = 1:length(subs) + % Loop over sessions + sess = dir([subs(kk).folder filesep subs(kk).name]); + sess(1:2) = []; + sess = sess([sess.isdir]); + sess = sess(contains({sess.name},'ses')); + for ll = 1:length(sess) + + % Specify metabolite data + % (MANDATORY) + dir_metabolite = dir([sess(ll).folder filesep sess(ll).name filesep 'mrs' filesep subs(kk).name '_' sess(ll).name '_megapress' filesep '*.SDAT']); + files(counter) = {[dir_metabolite(end).folder filesep dir_metabolite(end).name]}; + + % Specify water reference data for eddy-current correction (same sequence as metabolite data!) + % (OPTIONAL) + % Leave empty for GE P-files (.7) - these include water reference data by + % default. + dir_ref = dir([sess(ll).folder filesep sess(ll).name filesep 'mrs' filesep subs(kk).name '_' sess(ll).name '_megapress-ref' filesep '*.SDAT']); + files_ref(counter) = {[dir_ref(end).folder filesep dir_ref(end).name]}; + + % Specify water data for quantification (e.g. short-TE water scan) + % (OPTIONAL) + dir_w = dir([sess(ll).folder filesep sess(ll).name filesep 'mrs' filesep subs(kk).name '_' sess(ll).name '_press-ref' filesep '*.SDAT']); + files_w(counter) = {[dir_w(end).folder filesep dir_w(end).name]}; + + % Specify metabolite-nulled data for quantification + % (OPTIONAL) + files_mm = {}; + + % Specify T1-weighted structural imaging data + % (OPTIONAL) + % Link to single NIfTI (*.nii) files for Siemens and Philips data + % Link to DICOM (*.dcm) folders for GE data + files_nii(counter) = {[sess(ll).folder filesep sess(ll).name filesep 'anat' filesep subs(kk).name filesep subs(kk).name '_' sess(ll).name '_T1w.nii.gz']}; + + % External segmentation results + % (OPTIONAL) + % Link to NIfTI (*.nii or *.nii.gz) files with segmentation results + % Add supply gray matter, white matter, and CSF as 1 x 3 cell within a + % cell array or a single 4D file in the same order supplied as 1 x 1 cell; +% files_seg(counter) = {{[sess(ll).folder filesep sess(ll).name filesep 'anat' filesep subs(kk).name filesep 'c1' sess(ll).name '_T1w.nii.gz'],... +% [sess(ll).folder filesep sess(ll).name filesep 'anat' filesep subs(kk).name filesep 'c2' sess(ll).name '_T1w.nii.gz'],... +% [sess(ll).folder filesep sess(ll).name filesep 'anat' filesep subs(kk).name filesep 'c3' sess(ll).name '_T1w.nii.gz']}}; + +% files_seg(counter) = {{[sess(ll).folder filesep sess(ll).name filesep 'anat' filesep subs(kk).name filesep '4D' sess(ll).name '_T1w.nii.gz']}}; + + counter = counter + 1; + end +end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Definitions without using BIDS + +% You can always supply direct path to each of the files within +% the cell array. For example: + +% Specify metabolite data +% (MANDATORY) +% files(counter) = {'/Volumes/MyProject/data/sub-01/mrs/MEGAPRESS_act.SDAT',... +% '/Volumes/MyProject/data/sub-02/mrs/MEGAPRESS_act.SDAT'}; + +% Specify water reference data for eddy-current correction (same sequence as metabolite data!) +% (OPTIONAL) +% Leave empty for GE P-files (.7) - these include water reference data by +% default. +% files_ref(counter) = {'/Volumes/MyProject/data/sub-01/mrs/MEGAPRESS_ref.SDAT',... +% '/Volumes/MyProject/data/sub-02/mrs/MEGAPRESS_ref.SDAT'}; + +% Specify water data for quantification (e.g. short-TE water scan) +% (OPTIONAL) +% files_w = = {'/Volumes/MyProject/data/sub-01/mrs/PRESS_ref.SDAT',... +% '/Volumes/MyProject/data/sub-02/mrs/PRESS_ref.SDAT'}; + +% Specify metabolite-nulled data for quantification +% (OPTIONAL) +% files_mm = {}; + +% Specify T1-weighted structural imaging data +% (OPTIONAL) +% Link to single NIfTI (*.nii.gz or #.nii) files for GE, Siemens and Philips data +% files_nii = {'/Volumes/MyProject/data/sub-01/anat/T1w.nii.gz',... +% '/Volumes/MyProject/data/sub-02/anat/T1w.nii.gz'}; + +% External segmentation results +% (OPTIONAL) +% Link to NIfTI (*.nii or *.nii.gz) files with segmentation results +% Add supply gray matter, white matter, and CSF as 1 x 3 cell within a +% cell array or a single 4D file in the same order supplied as 1 x 1 cell; +% files_seg(counter) = {{'/Volumes/MyProject/data/sub-01/anat/c1T1w.nii.gz',... +% '/Volumes/MyProject/data/sub-01/anat/c2T1w.nii.gz',... +% '/Volumes/MyProject/data/sub-01/anat/c3T1w.nii.gz'},... +% {'/Volumes/MyProject/data/sub-02/anat/c1T1w.nii.gz',... +% '/Volumes/MyProject/data/sub-02/anat/c2T1w.nii.gz',... +% '/Volumes/MyProject/data/sub-02/anat/c3T1w.nii.gz'}}; +% files_seg(counter) = {{'/Volumes/MyProject/data/sub-01/anat/4DT1w.nii.gz'},... +% {'/Volumes/MyProject/data/sub-02/anat/4DT1w.nii.gz'}}; + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% 4. SPECIFY STAT FILE %%% +% Supply location of a csv file, which contains possible correlation +% measures and group variables. Each column must start with the name of the +% measure. For the grouping variable use 'group' and numbers between 1 and +% the number of included groups. If no group is supplied the data will be +% treated as one group. (You can always use the direct path) + +file_stat = fullfile(data_folder, 'stat.csv'); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% 5. SPECIFY OUTPUT FOLDER %% +% The Osprey data container will be saved as a *.mat file in the output +% folder that you specify below. In addition, any exported files (for use +% with jMRUI, TARQUIN, or LCModel) will be saved in sub-folders. + +% Specify output folder (you can always use the direct path) +% (MANDATORY) +outputFolder = fullfile(data_folder, 'derivativesLCM'); diff --git a/fit/code/LCMwrapper/LCModelWrapper.m b/fit/code/LCMwrapper/LCModelWrapper.m index 6b2ae753..110b077d 100644 --- a/fit/code/LCMwrapper/LCModelWrapper.m +++ b/fit/code/LCMwrapper/LCModelWrapper.m @@ -65,9 +65,10 @@ % Save the parameters and information about the basis set - MRSCont.fit.results.metab.fitParams{kk} = readLCMFitParams(MRSCont, 'A', kk); + MRSCont.fit.results.metab.fitParams{1,kk,1} = readLCMFitParams(MRSCont, 'A', kk); if MRSCont.flags.isMEGA callLCModel(MRSCont, MRSCont.opts.fit.lcmodel.controlfileDiff1{kk},[pathLCModelBinary filesep bin]); + MRSCont.fit.results.metab.fitParams{1,kk,2} = readLCMFitParams(MRSCont, 'diff1', kk); end end \ No newline at end of file diff --git a/fit/code/LCMwrapper/readLCMFitParams.m b/fit/code/LCMwrapper/readLCMFitParams.m index 33f44129..5ce90539 100644 --- a/fit/code/LCMwrapper/readLCMFitParams.m +++ b/fit/code/LCMwrapper/readLCMFitParams.m @@ -24,7 +24,7 @@ end % Read the fit results from the .table files -tab = mrs_readLcmodelTABLE([MRSCont.opts.fit.lcmodel.outputfileA{kk} '.table']); +tab = mrs_readLcmodelTABLE([MRSCont.opts.fit.lcmodel.(lcmOutputFile){kk} '.table']); fitParams.name = tab.name; fitParams.CRLB = tab.SDpct; fitParams.relConc = tab.relative_conc; @@ -60,8 +60,8 @@ fitParams.residual = fitParams.data - fitParams.completeFit; % Add Nan values for nicer plots -if ~isempty(MRSCont.opts.fit.GAP.A) - idx = find(x_ppm < MRSCont.opts.fit.GAP.A(1)); +if ~isempty(MRSCont.opts.fit.GAP.(which)) + idx = find(x_ppm < MRSCont.opts.fit.GAP.(which)(1)); idx_1 = idx(1); idx_2 = idx_1-1; fitParams.data(idx_1) = NaN; @@ -88,7 +88,7 @@ if idxMatch fitParams.indivMets(:,rr) = spectra_metabolites(:,idxMatch) - fitParams.baseline; else - if ~isempty(MRSCont.opts.fit.GAP.A) + if ~isempty(MRSCont.opts.fit.GAP.(which)) fitParams.indivMets(idx_1,rr) = NaN; fitParams.indivMets(idx_2,rr) = NaN; end diff --git a/fit/code/osp_fitInitialise.m b/fit/code/osp_fitInitialise.m index 13c2f095..0c93f12b 100644 --- a/fit/code/osp_fitInitialise.m +++ b/fit/code/osp_fitInitialise.m @@ -317,7 +317,20 @@ MRSCont.opts.fit = rmfield(MRSCont.opts.fit,'basisSetFile'); MRSCont.opts.fit.basisSetFile{1} = [path filesep Bo '_' seq '_' MRSCont.vendor '_' te 'ms_noMM_A.BASIS']; if MRSCont.flags.isMEGA + if strcmp(MRSCont.opts.fit.coMM3,'1to1GABA') + MRSCont.opts.fit.CrFactor = 1; + basisSetDiff1 = basisSet; + basisSetDiff1.fids = basisSetDiff1.fids(:,:,3); + basisSetDiff1.specs = basisSetDiff1.specs(:,:,3); + basisSetOff = basisSet; + basisSetOff.fids = basisSetOff.fids(:,:,1); + basisSetOff.specs = basisSetOff.specs(:,:,1); + [basisSetDiff1] = osp_addDiffMMPeaks(basisSetDiff1,basisSetOff,MRSCont.opts.fit); + basisSet.fids(:,:,3)=basisSetDiff1.fids; + basisSet.specs(:,:,3)=basisSetDiff1.specs; + end [~] = io_writelcmBASIS(basisSet,[path filesep Bo '_' seq '_' MRSCont.vendor '_' te 'ms_noMM_diff1.BASIS'], MRSCont.vendor, seq, 3); + MRSCont.opts.fit.basisSetFile{2} = [path filesep Bo '_' seq '_' MRSCont.vendor '_' te 'ms_noMM_diff1.BASIS']; end elseif strcmpi(exten, '.basis') @@ -335,8 +348,13 @@ metabsToInclude{1} = fit_createMetabList(MRSCont.opts.fit.includeMetabs); metabsInBasis{1} = fit_readLCMBasisSetMetabs(basisSetFile{1}); if MRSCont.flags.isMEGA - % metabsToInclude{2} = metabsToInclude{1}; - metabsToInclude{2} = fit_createMetabList({'GABA','GSH','Gln','Glu','NAAG','NAA'}); + basisSetFile{2} = MRSCont.opts.fit.basisSetFile{2}; + switch MRSCont.opts.fit.coMM3 + case {'3to2MM','none','1to1GABA'} + metabsToInclude{2} = fit_createMetabList({'GABA','GSH','Gln','Glu','NAAG','NAA','MM09'}); + case {'fixedGauss','1to1GABAsoft','3to2MMsoft'} + metabsToInclude{2} = fit_createMetabList({'GABA','GSH','Gln','Glu','NAAG','NAA','MM09','MM30'}); + end metabsInBasis{2} = fit_readLCMBasisSetMetabs(basisSetFile{2}); end % Loop over metabolites in the basis set @@ -358,6 +376,30 @@ end end + % MM and Lipids may not be part of the basis set file but + % still be ignored + MMList = {'MM09','MM12','MM14','MM17','MM20','MM22','MM27','MM30','MM32','Lip09','Lip20',... + 'MM37','MM38','MM40','MM42'}; + if ~metabsToInclude{nn}.Lip13 + MMList{end+1} = 'Lip13a'; + MMList{end+1} = 'Lip13b'; + end + for qq = 1:length(MMList) + % Take current name + currentName = MMList{qq}; + % Locate it in the list + idx = find(ismember(fieldnames(metabsToInclude{nn}),currentName)); + if ~isempty(idx) + % If it's a match, check whether it should be included + if metabsToInclude{nn}.(currentName) == 1 + else + chOmitList{1,nn}{end+1} = ['''' currentName '''']; + end + else + chOmitList{1,nn}{end+1} = ['''' currentName '''']; + end + + end end end @@ -376,7 +418,7 @@ LCMparam = osp_editControlParameters(LCMparam, 'ltable', '7'); LCMparam = osp_editControlParameters(LCMparam, 'lcsv', '11'); LCMparam = osp_editControlParameters(LCMparam, 'neach', '99'); - LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau'''}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau'''}); % Augment list of metabolites to omit if isfield(LCMparam, 'chomit') chOmitList{1} = unique(horzcat(chOmitList{1}, LCMparam.chomit)); @@ -477,20 +519,52 @@ LCMparam = osp_editControlParameters(LCMparam, 'key', '210387309'); LCMparam = osp_editControlParameters(LCMparam, 'owner', '''Osprey processed spectra'''); LCMparam = osp_editControlParameters(LCMparam, 'filbas', ['''' basisSetFile{2} '''']); - % LCMparam = osp_editControlParameters(LCMparam, 'dkntmn', '0.15'); + LCMparam = osp_editControlParameters(LCMparam, 'dkntmn', '0.15'); LCMparam = osp_editControlParameters(LCMparam, 'neach', '99'); LCMparam = osp_editControlParameters(LCMparam, 'sptype', '''mega-press-3'''); - LCMparam = osp_editControlParameters(LCMparam, 'nobase', 'T'); - %LCMparam = osp_editControlParameters(LCMparam, 'wdline', '0'); - LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '1'); - LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {'''MM09 @ 0.915 +- .02 FWHM= .14 < .17 +- .015 AMP= 3. @ 3.0 FWHM= .19 AMP= 2.'''}); - LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau'''}); + if ~isinf(MRSCont.opts.fit.bLineKnotSpace) + LCMparam = osp_editControlParameters(LCMparam, 'nobase', 'F'); + else + LCMparam = osp_editControlParameters(LCMparam, 'dkntmn', '0.6'); + end + switch MRSCont.opts.fit.coMM3 + case {'3to2MM'} + LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '1'); + LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {['''' sprintf('MM09 @ 0.915 +- .02 FWHM= .085 < .1 +- .35 AMP= 3. @ 3.0 FWHM= %4.2f AMP= 2.',MRSCont.opts.fit.FWHMcoMM3/(MRSCont.processed.metab{kk}.txfrq/1000000)) '''']}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau''','''GABA+MM09'''}); + case {'none','1to1GABA'} + LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '1'); + LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {'''MM09 @ 0.915 +- .02 FWHM= .085 < .1 +- .35 AMP= 3.'''}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau'''}); + case {'fixedGauss'} + LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '2'); + LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {'''MM09 @ 0.915 +- .02 FWHM= .085 < .1 +- .35 AMP= 3.''',... + ['''' sprintf('MM30 @ 3.0 +- .02 FWHM= .085 < %4.2f +- .35 AMP= 2.',MRSCont.opts.fit.FWHMcoMM3/(MRSCont.processed.metab{kk}.txfrq/1000000)) '''']}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau''','''GABA+MM30'''}); + case {'1to1GABAsoft'} + LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '2'); + LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {'''MM09 @ 0.915 +- .02 FWHM= .085 < .1 +- .35 AMP= 3.''',... + ['''' sprintf('MM30 @ 3.0 +- .02 FWHM= .085 < %4.2f +- .35 AMP= 2.',MRSCont.opts.fit.FWHMcoMM3/(MRSCont.processed.metab{kk}.txfrq/1000000)) '''']}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau''','''GABA+MM30'''}); + LCMparam = osp_editControlParameters(LCMparam, 'chrato', {'''GABA/MM30 = 1.0 +- .1'''}); + case {'3to2MMsoft'} + LCMparam = osp_editControlParameters(LCMparam, 'nsimul', '2'); + LCMparam = osp_editControlParameters(LCMparam, 'chsimu', {'''MM09 @ 0.915 +- .02 FWHM= .085 < .1 +- .35 AMP= 3.''',... + ['''' sprintf('MM30 @ 3.0 +- .02 FWHM= .085 < %4.2f +- .35 AMP= 2.',MRSCont.opts.fit.FWHMcoMM3/(MRSCont.processed.metab{kk}.txfrq/1000000)) '''']}); + LCMparam = osp_editControlParameters(LCMparam, 'chcomb', {'''PCh+GPC''','''Cr+PCr''','''NAA+NAAG''','''Glu+Gln''','''Glc+Tau''','''GABA+MM30'''}); + LCMparam = osp_editControlParameters(LCMparam, 'chrato', {'''MM30/MM09 = 0.66 +- .2'''}); + end + LCMparam = osp_editControlParameters(LCMparam, 'chomit', chOmitList{2}); - LCMparam = osp_editControlParameters(LCMparam, 'namrel', '''Cr+PCr'''); + LCMparam = osp_editControlParameters(LCMparam, 'namrel', '''NAA+NAAG'''); LCMparam = osp_editControlParameters(LCMparam, 'ppmst', ['' sprintf('%4.2f', MRSCont.opts.fit.range(2)) '']); LCMparam = osp_editControlParameters(LCMparam, 'ppmend', ['' sprintf('%4.2f', MRSCont.opts.fit.range(1)) '']); - LCMparam = osp_editControlParameters(LCMparam, 'ppmgap11', ['' sprintf('%4.2f', 1.95) '']); - LCMparam = osp_editControlParameters(LCMparam, 'ppmgap21', ['' sprintf('%4.2f', 1.2) '']); + if isfield(MRSCont.opts.fit, 'GAP') && isfield(MRSCont.opts.fit.GAP, 'diff1') + if ~isempty(MRSCont.opts.fit.GAP.diff1) + LCMparam = osp_editControlParameters(LCMparam, 'ppmgap11', ['' sprintf('%4.2f', MRSCont.opts.fit.GAP.diff1(2)) '']); + LCMparam = osp_editControlParameters(LCMparam, 'ppmgap21', ['' sprintf('%4.2f', MRSCont.opts.fit.GAP.diff1(1)) '']); + end + end % Add water-scaling-related flags only if water reference % data has been provided if MRSCont.flags.hasRef || MRSCont.flags.hasWater diff --git a/fit/code/osp_fitMEGA.m b/fit/code/osp_fitMEGA.m index d01cf12b..94718493 100644 --- a/fit/code/osp_fitMEGA.m +++ b/fit/code/osp_fitMEGA.m @@ -54,242 +54,253 @@ end for kk = 1:MRSCont.nDatasets(1) - [~] = printLog('OspreyFit',kk,1,MRSCont.nDatasets,progressText,MRSCont.flags.isGUI ,MRSCont.flags.isMRSI); - - %%% 1. DETERMINE THE FITTING STYLE %%% - % Extract fit options - fitOpts = MRSCont.opts.fit; - fitModel = fitOpts.method; - fitStyle = fitOpts.style; - - %%% 2. SEPARATE FIT %%% - % For the separate (classic) MEGA fit, model the EDIT-OFF and DIFF - % spectra separately. - if strcmp(fitStyle, 'Separate') - if ~(MRSCont.flags.didFit == 1 && MRSCont.flags.speedUp && isfield(MRSCont, 'fit') && (kk > size(MRSCont.fit.results.metab.fitParams,2))) || ~strcmp(MRSCont.ver.Osp,MRSCont.ver.CheckOsp) - %%% 2a. FIT OFF-SPECTRUM - % Apply scaling factor to the data - dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'A'); - basisSetOff = MRSCont.fit.basisSet; - basisSetOff.names{1} = 'A'; - basisSetOff.fids = basisSetOff.fids(:,:,1); - basisSetOff.specs = basisSetOff.specs(:,:,1); - dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); - fitOpts.GAP = MRSCont.opts.fit.GAP.A; - - - % Call the fit function - [fitParamsOff, resBasisSetOff] = fit_runFit(dataToFit, basisSetOff, fitModel, fitOpts); - - % Save back the basis set and fit parameters to MRSCont - MRSCont.fit.resBasisSet.metab{1,kk,1} = resBasisSetOff; - MRSCont.fit.results.metab.fitParams{1,kk,1} = fitParamsOff; - - - %%% 2b. FIT DIFF1-SPECTRUM - % Apply scaling factor to the data - fitOpts = MRSCont.opts.fit; - dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'diff1'); - basisSetDiff1 = MRSCont.fit.basisSet; - basisSetDiff1.fids = basisSetDiff1.fids(:,:,3); - basisSetDiff1.specs = basisSetDiff1.specs(:,:,3); - dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); - dataToFit.refShift = fitParamsOff.refShift; - dataToFit.refFWHM = fitParamsOff.refFWHM; - - if ~strcmp(fitOpts.coMM3, 'none') - fitOpts.CrFactor = 1; - [basisSetDiff1] = osp_addDiffMMPeaks(basisSetDiff1,basisSetOff,fitOpts); - end - - fitOpts.GAP = MRSCont.opts.fit.GAP.diff1; - basisSetDiff1.names{1} = 'diff1'; - % Call the fit function - [fitParamsDiff1, resBasisSetDiff1] = fit_runFit(dataToFit, basisSetDiff1, fitModel, fitOpts); - - % Save back the basis set and fit parameters to MRSCont - MRSCont.fit.resBasisSet.metab{1,kk,2} = resBasisSetDiff1; - MRSCont.fit.results.metab.fitParams{1,kk,2} = fitParamsDiff1; - - %Modeling MM spectra after the main spectrum - if MRSCont.flags.hasMM == 1 - dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'A'); - dataToFit_mm = op_ampScale(dataToFit_mm, 1/MRSCont.fit.scale{kk}); - [refShift_mm, ~] = fit_OspreyReferencingMM(dataToFit_mm); - [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); - %add some info from the metabolite fit - dataToFit_mm.lineShape = fitParamsOff.lineShape; - dataToFit_mm.refFWHM = fitParamsOff.refFWHM; - dataToFit_mm.refShift = 0; - % Extract fit options - fitOpts_mm = MRSCont.opts.fit; - fitModel_mm = fitOpts.method; - fitOpts_mm.sequence = 'unedited'; - - % Call the fit function - [fitParams_mm, resBasisSet_mm] = fit_runFitMM(dataToFit_mm, basisSet_mm_A, fitModel_mm, fitOpts_mm); - % Save back the basis set and fit parameters to MRSCont - MRSCont.fit.basisSet_mm = basisSet_mm_A; - MRSCont.fit.resBasisSet.mm{1,kk,1} = resBasisSet_mm; - MRSCont.fit.results.mm.fitParams{1,kk,1} = fitParams_mm; - - % Create full spectralwidth clean MM spectrum - disp('Create clean metabolite-nulled off spectrum...'); - [mm_clean] = fit_OspreyCleanMM(dataToFit_mm, resBasisSet_mm, fitOpts_mm,MRSCont.fit.scale{kk},fitParams_mm); - mm_clean.names = {'A_clean'}; - mm_clean.flags.isUnEdited = 1; - mm_clean.flags.isMEGA = 0; -% [mm_clean,~] = op_autophase(mm_clean,0.8,1); - [refShift_mm, ~] = fit_OspreyReferencingMM(mm_clean); - [mm_clean] = op_freqshift(mm_clean,-refShift_mm); - mm_clean = op_ampScale(mm_clean, MRSCont.fit.scale{kk}); - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean); - dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'A'); - [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},dataToFit_mm); - - - % Create full spectralwidth clean MM spectrum spline model - disp('Create clean metabolite-nulled off spline model...'); - %Fit them with a spline? - dataToFit_mm_clean = op_takesubspec(MRSCont.processed.mm{kk},'A_clean'); - fitOpts_mm_clean = MRSCont.opts.fit; - fitOpts_mm_clean.range=[dataToFit_mm_clean.ppm(1) dataToFit_mm_clean.ppm(end)]; - [~,mm_clean_spline] = fit_Osprey_SplineOnly(dataToFit_mm_clean, 0.1, fitOpts_mm_clean.range); - mm_clean_spline.names = {'A_spline'}; - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean_spline); - -% % Gaussian fit with(out) baseline -% disp('Fit clean metabolite-nulled spectrum...'); -% fitOpts_mm.Params2 = fitParams_mm; -% dataToFit_mm_clean = op_ampScale(dataToFit_mm_clean, 1/MRSCont.fit.scale{kk}); -% [fitParams_mm_Gauss, resBasisSet_mm_Gauss] = fit_runFitMM(dataToFit_mm_clean, basisSetMM, 'OspreyMMGaussian', fitOpts_mm); -% MRSCont.fit.basisSet_mm_Gauss = basisSetMM; -% MRSCont.fit.resBasisSet.diff1_mm_Gauss{kk} = resBasisSet_mm_Gauss; -% MRSCont.fit.results.diff1_mm.fitParams{kk}.fitParams_mm_Gauss = fitParams_mm_Gauss; - -% % Fit off with individual MM model - disp('Fit off with metabolite-nulled spline model...'); - fitOpts.mm_clean_spline = mm_clean_spline; - fitOpts.prelimParams = fitParamsOff.prelimParams; + % ----- Osprey fit pipeline ----- + if strcmpi(MRSCont.opts.fit.method, 'Osprey') + [~] = printLog('OspreyFit',kk,1,MRSCont.nDatasets,progressText,MRSCont.flags.isGUI ,MRSCont.flags.isMRSI); + + %%% 1. DETERMINE THE FITTING STYLE %%% + % Extract fit options + fitOpts = MRSCont.opts.fit; + fitModel = fitOpts.method; + fitStyle = fitOpts.style; + + %%% 2. SEPARATE FIT %%% + % For the separate (classic) MEGA fit, model the EDIT-OFF and DIFF + % spectra separately. + if strcmp(fitStyle, 'Separate') + if ~(MRSCont.flags.didFit == 1 && MRSCont.flags.speedUp && isfield(MRSCont, 'fit') && (kk > size(MRSCont.fit.results.metab.fitParams,2))) || ~strcmp(MRSCont.ver.Osp,MRSCont.ver.CheckOsp) + %%% 2a. FIT OFF-SPECTRUM + % Apply scaling factor to the data dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'A'); + basisSetOff = MRSCont.fit.basisSet; + basisSetOff.names{1} = 'A'; + basisSetOff.fids = basisSetOff.fids(:,:,1); + basisSetOff.specs = basisSetOff.specs(:,:,1); dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); + fitOpts.GAP = MRSCont.opts.fit.GAP.A; + + % Call the fit function - [fitParamsOffExpMM, resBasisSetOffExpMM] = fit_runFit(dataToFit, basisSetOff, fitModel, fitOpts); - MRSCont.fit.resBasisSet.metab{2,kk,1} = resBasisSetOffExpMM; - MRSCont.fit.results.metab.fitParams{2,kk,1} = fitParamsOffExpMM; - - %Clean diff1 mm spectrum - dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'diff1'); - dataToFit_mm = op_ampScale(dataToFit_mm, 1/MRSCont.fit.scale{kk}); - [refShift_mm, ~] = fit_OspreyReferencingMM(dataToFit_mm); - [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); - %add some info from the metabolite fit - dataToFit_mm.lineShape = fitParamsDiff1.lineShape; - dataToFit_mm.refFWHM = fitParamsDiff1.refFWHM; - dataToFit_mm.refShift = 0; - % Extract fit options - fitOpts_mm = MRSCont.opts.fit; - fitModel_mm = fitOpts.method; - fitOpts_mm.sequence = 'MEGA'; - - % Call the fit function - [fitParams_mm, resBasisSet_mm] = fit_runFitMM(dataToFit_mm, basisSet_mm_diff1, fitModel_mm, fitOpts_mm); - % Save back the basis set and fit parameters to MRSCont - MRSCont.fit.basisSet_mm = basisSet_mm_diff1; - MRSCont.fit.resBasisSet.mm{1,kk,2} = resBasisSet_mm; - MRSCont.fit.results.mm.fitParams{1,kk,2} = fitParams_mm; - - % Create full spectralwidth clean MM spectrum - disp('Create clean metabolite-nulled diff1 spectrum...'); - [mm_clean] = fit_OspreyCleanMM(dataToFit_mm, resBasisSet_mm, fitOpts_mm,MRSCont.fit.scale{kk},fitParams_mm); - mm_clean.names = {'diff1_clean'}; -% [mm_clean,~] = op_autophase(mm_clean,0.8,1); - [refShift_mm, ~] = fit_OspreyReferencingMM(mm_clean); - [mm_clean] = op_freqshift(mm_clean,-refShift_mm); - mm_clean = op_ampScale(mm_clean, MRSCont.fit.scale{kk}); - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean); - dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'diff1'); - [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},dataToFit_mm); - - - % Create full spectralwidth clean MM spectrum spline model - disp('Create clean metabolite-nulled diff1 spline model...'); - %Fit them with a spline? - dataToFit_mm_clean = op_takesubspec(MRSCont.processed.mm{kk},'diff1_clean'); - fitOpts_mm_clean = MRSCont.opts.fit; - fitOpts_mm_clean.range=[dataToFit_mm_clean.ppm(1) dataToFit_mm_clean.ppm(end)]; - [~,mm_clean_spline] = fit_Osprey_SplineOnly(dataToFit_mm_clean, 0.1, fitOpts_mm_clean.range); - mm_clean_spline.names = {'diff1_spline'}; - MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean_spline); - -% % Gaussian fit with(out) baseline -% disp('Fit clean metabolite-nulled spectrum...'); -% fitOpts_mm.Params2 = fitParams_mm; -% dataToFit_mm_clean = op_ampScale(dataToFit_mm_clean, 1/MRSCont.fit.scale{kk}); -% [fitParams_mm_Gauss, resBasisSet_mm_Gauss] = fit_runFitMM(dataToFit_mm_clean, basisSetMM, 'OspreyMMGaussian', fitOpts_mm); -% MRSCont.fit.basisSet_mm_Gauss = basisSetMM; -% MRSCont.fit.resBasisSet.diff1_mm_Gauss{kk} = resBasisSet_mm_Gauss; -% MRSCont.fit.results.diff1_mm.fitParams{kk}.fitParams_mm_Gauss = fitParams_mm_Gauss; - -% % Fit diff1 with individual MM model - disp('Fit diff1 with metabolite-nulled spline model...'); - fitOpts.mm_clean_spline = mm_clean_spline; - fitOpts.prelimParams = fitParamsDiff1.prelimParams; - fitOpts.coMM3 = 'none'; + [fitParamsOff, resBasisSetOff] = fit_runFit(dataToFit, basisSetOff, fitModel, fitOpts); + + % Save back the basis set and fit parameters to MRSCont + MRSCont.fit.resBasisSet.metab{1,kk,1} = resBasisSetOff; + MRSCont.fit.results.metab.fitParams{1,kk,1} = fitParamsOff; + + + %%% 2b. FIT DIFF1-SPECTRUM + % Apply scaling factor to the data + fitOpts = MRSCont.opts.fit; dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'diff1'); - dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); - % Call the fit function - [fitParamsDiff1ExpMM, resBasisSetDiff1ExpMM] = fit_runFit(dataToFit, basisSetDiff1, fitModel, fitOpts); - MRSCont.fit.resBasisSet.metab{2,kk,2} = resBasisSetDiff1ExpMM; - MRSCont.fit.results.metab.fitParams{2,kk,2} = fitParamsDiff1ExpMM; - - end - - end - end - - - %%% 3. CONCATENATED FIT %%% - % For the concatenated MEGA fit, model the DIFF1 and SUM spectra - % together. - if strcmp(fitStyle, 'Concatenated') - if ~(MRSCont.flags.didFit == 1 && MRSCont.flags.speedUp && isfield(MRSCont, 'fit') && (kk > length(MRSCont.fit.results.conc.fitParams))) || ~strcmp(MRSCont.ver.Osp,MRSCont.ver.CheckOsp) - %%% 3a. FIT CONCATENATED SPECTRUM - % Select the difference and sum spectrum to put into the - % concatenated fit - dataToFit = {op_takesubspec(MRSCont.processed.metab{kk},3) op_takesubspec(MRSCont.processed.metab{kk},4)}; - basisSetConc = MRSCont.fit.basisSet; - % Create the basis set with difference and sum basis functions - basisSetConc.fids(:,:,1) = basisSetConc.fids(:,:,3); - basisSetConc.specs(:,:,1) = basisSetConc.specs(:,:,3); - basisSetConc.fids(:,:,2) = basisSetConc.fids(:,:,4); - basisSetConc.specs(:,:,2) = basisSetConc.specs(:,:,4); - basisSetConc.fids(:,:,3:4) = []; - basisSetConc.specs(:,:,3:4) = []; - - if ~strcmp(fitOpts.coMM3, 'none') basisSetDiff1 = MRSCont.fit.basisSet; basisSetDiff1.fids = basisSetDiff1.fids(:,:,3); basisSetDiff1.specs = basisSetDiff1.specs(:,:,3); - fitOpts.CrFactor = 1; - [basisSetDiff1] = osp_addDiffMMPeaks(basisSetDiff1,fitOpts); - basisSetConc.fids(:,:,1) = basisSetDiff1.fids(:,:); - basisSetConc.specs(:,:,1) = basisSetDiff1.specs(:,:); - end + dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); + dataToFit.refShift = fitParamsOff.refShift; + dataToFit.refFWHM = fitParamsOff.refFWHM; + + if ~strcmp(fitOpts.coMM3, 'none') + fitOpts.CrFactor = 1; + [basisSetDiff1] = osp_addDiffMMPeaks(basisSetDiff1,basisSetOff,fitOpts); + end + + fitOpts.GAP = MRSCont.opts.fit.GAP.diff1; + basisSetDiff1.names{1} = 'diff1'; - % Apply scaling factor to the data - for rr = 1:length(dataToFit) - dataToFit{rr} = op_ampScale(dataToFit{rr}, 1/MRSCont.fit.scale{kk}); - end - % Call the multi-spectrum fit function - [fitParamsConc, resBasisSetConc] = fit_runFitMultiplex(dataToFit, basisSetConc, fitModel, fitOpts); + % metabList = fit_createMetabList({'GABA','GSH','Gln','Glu','NAAG','NAA','MM09'}); + % Create the modified basis set + % basisSetDiff1 = fit_selectMetabs(basisSetDiff1, metabList, 1); - % Save back the basis set and fit parameters to MRSCont - MRSCont.fit.resBasisSet.conc{kk} = resBasisSetConc; - MRSCont.fit.results.conc.fitParams{kk} = fitParamsConc; + % Call the fit function + [fitParamsDiff1, resBasisSetDiff1] = fit_runFit(dataToFit, basisSetDiff1, fitModel, fitOpts); + + % Save back the basis set and fit parameters to MRSCont + MRSCont.fit.resBasisSet.metab{1,kk,2} = resBasisSetDiff1; + MRSCont.fit.results.metab.fitParams{1,kk,2} = fitParamsDiff1; + + %Modeling MM spectra after the main spectrum + if MRSCont.flags.hasMM == 1 + dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'A'); + dataToFit_mm = op_ampScale(dataToFit_mm, 1/MRSCont.fit.scale{kk}); + [refShift_mm, ~] = fit_OspreyReferencingMM(dataToFit_mm); + [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); + %add some info from the metabolite fit + dataToFit_mm.lineShape = fitParamsOff.lineShape; + dataToFit_mm.refFWHM = fitParamsOff.refFWHM; + dataToFit_mm.refShift = 0; + % Extract fit options + fitOpts_mm = MRSCont.opts.fit; + fitModel_mm = fitOpts.method; + fitOpts_mm.sequence = 'unedited'; + + % Call the fit function + [fitParams_mm, resBasisSet_mm] = fit_runFitMM(dataToFit_mm, basisSet_mm_A, fitModel_mm, fitOpts_mm); + % Save back the basis set and fit parameters to MRSCont + MRSCont.fit.basisSet_mm = basisSet_mm_A; + MRSCont.fit.resBasisSet.mm{1,kk,1} = resBasisSet_mm; + MRSCont.fit.results.mm.fitParams{1,kk,1} = fitParams_mm; + + % Create full spectralwidth clean MM spectrum + disp('Create clean metabolite-nulled off spectrum...'); + [mm_clean] = fit_OspreyCleanMM(dataToFit_mm, resBasisSet_mm, fitOpts_mm,MRSCont.fit.scale{kk},fitParams_mm); + mm_clean.names = {'A_clean'}; + mm_clean.flags.isUnEdited = 1; + mm_clean.flags.isMEGA = 0; + % [mm_clean,~] = op_autophase(mm_clean,0.8,1); + [refShift_mm, ~] = fit_OspreyReferencingMM(mm_clean); + [mm_clean] = op_freqshift(mm_clean,-refShift_mm); + mm_clean = op_ampScale(mm_clean, MRSCont.fit.scale{kk}); + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean); + dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'A'); + [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},dataToFit_mm); + + + % Create full spectralwidth clean MM spectrum spline model + disp('Create clean metabolite-nulled off spline model...'); + %Fit them with a spline? + dataToFit_mm_clean = op_takesubspec(MRSCont.processed.mm{kk},'A_clean'); + fitOpts_mm_clean = MRSCont.opts.fit; + fitOpts_mm_clean.range=[dataToFit_mm_clean.ppm(1) dataToFit_mm_clean.ppm(end)]; + [~,mm_clean_spline] = fit_Osprey_SplineOnly(dataToFit_mm_clean, 0.1, fitOpts_mm_clean.range); + mm_clean_spline.names = {'A_spline'}; + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean_spline); + + % % Gaussian fit with(out) baseline + % disp('Fit clean metabolite-nulled spectrum...'); + % fitOpts_mm.Params2 = fitParams_mm; + % dataToFit_mm_clean = op_ampScale(dataToFit_mm_clean, 1/MRSCont.fit.scale{kk}); + % [fitParams_mm_Gauss, resBasisSet_mm_Gauss] = fit_runFitMM(dataToFit_mm_clean, basisSetMM, 'OspreyMMGaussian', fitOpts_mm); + % MRSCont.fit.basisSet_mm_Gauss = basisSetMM; + % MRSCont.fit.resBasisSet.diff1_mm_Gauss{kk} = resBasisSet_mm_Gauss; + % MRSCont.fit.results.diff1_mm.fitParams{kk}.fitParams_mm_Gauss = fitParams_mm_Gauss; + + % % Fit off with individual MM model + disp('Fit off with metabolite-nulled spline model...'); + fitOpts.mm_clean_spline = mm_clean_spline; + fitOpts.prelimParams = fitParamsOff.prelimParams; + dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'A'); + dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); + % Call the fit function + [fitParamsOffExpMM, resBasisSetOffExpMM] = fit_runFit(dataToFit, basisSetOff, fitModel, fitOpts); + MRSCont.fit.resBasisSet.metab{2,kk,1} = resBasisSetOffExpMM; + MRSCont.fit.results.metab.fitParams{2,kk,1} = fitParamsOffExpMM; + + %Clean diff1 mm spectrum + dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'diff1'); + dataToFit_mm = op_ampScale(dataToFit_mm, 1/MRSCont.fit.scale{kk}); + [refShift_mm, ~] = fit_OspreyReferencingMM(dataToFit_mm); + [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); + %add some info from the metabolite fit + dataToFit_mm.lineShape = fitParamsDiff1.lineShape; + dataToFit_mm.refFWHM = fitParamsDiff1.refFWHM; + dataToFit_mm.refShift = 0; + % Extract fit options + fitOpts_mm = MRSCont.opts.fit; + fitModel_mm = fitOpts.method; + fitOpts_mm.sequence = 'MEGA'; + + % Call the fit function + [fitParams_mm, resBasisSet_mm] = fit_runFitMM(dataToFit_mm, basisSet_mm_diff1, fitModel_mm, fitOpts_mm); + % Save back the basis set and fit parameters to MRSCont + MRSCont.fit.basisSet_mm = basisSet_mm_diff1; + MRSCont.fit.resBasisSet.mm{1,kk,2} = resBasisSet_mm; + MRSCont.fit.results.mm.fitParams{1,kk,2} = fitParams_mm; + + % Create full spectralwidth clean MM spectrum + disp('Create clean metabolite-nulled diff1 spectrum...'); + [mm_clean] = fit_OspreyCleanMM(dataToFit_mm, resBasisSet_mm, fitOpts_mm,MRSCont.fit.scale{kk},fitParams_mm); + mm_clean.names = {'diff1_clean'}; + % [mm_clean,~] = op_autophase(mm_clean,0.8,1); + [refShift_mm, ~] = fit_OspreyReferencingMM(mm_clean); + [mm_clean] = op_freqshift(mm_clean,-refShift_mm); + mm_clean = op_ampScale(mm_clean, MRSCont.fit.scale{kk}); + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean); + dataToFit_mm = op_takesubspec(MRSCont.processed.mm{kk},'diff1'); + [dataToFit_mm] = op_freqshift(dataToFit_mm,-refShift_mm); + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},dataToFit_mm); + + + % Create full spectralwidth clean MM spectrum spline model + disp('Create clean metabolite-nulled diff1 spline model...'); + %Fit them with a spline? + dataToFit_mm_clean = op_takesubspec(MRSCont.processed.mm{kk},'diff1_clean'); + fitOpts_mm_clean = MRSCont.opts.fit; + fitOpts_mm_clean.range=[dataToFit_mm_clean.ppm(1) dataToFit_mm_clean.ppm(end)]; + [~,mm_clean_spline] = fit_Osprey_SplineOnly(dataToFit_mm_clean, 0.1, fitOpts_mm_clean.range); + mm_clean_spline.names = {'diff1_spline'}; + MRSCont.processed.mm{kk} = op_mergesubspec(MRSCont.processed.mm{kk},mm_clean_spline); + + % % Gaussian fit with(out) baseline + % disp('Fit clean metabolite-nulled spectrum...'); + % fitOpts_mm.Params2 = fitParams_mm; + % dataToFit_mm_clean = op_ampScale(dataToFit_mm_clean, 1/MRSCont.fit.scale{kk}); + % [fitParams_mm_Gauss, resBasisSet_mm_Gauss] = fit_runFitMM(dataToFit_mm_clean, basisSetMM, 'OspreyMMGaussian', fitOpts_mm); + % MRSCont.fit.basisSet_mm_Gauss = basisSetMM; + % MRSCont.fit.resBasisSet.diff1_mm_Gauss{kk} = resBasisSet_mm_Gauss; + % MRSCont.fit.results.diff1_mm.fitParams{kk}.fitParams_mm_Gauss = fitParams_mm_Gauss; + + % % Fit diff1 with individual MM model + disp('Fit diff1 with metabolite-nulled spline model...'); + fitOpts.mm_clean_spline = mm_clean_spline; + fitOpts.prelimParams = fitParamsDiff1.prelimParams; + fitOpts.coMM3 = 'none'; + dataToFit = op_takesubspec(MRSCont.processed.metab{kk},'diff1'); + dataToFit = op_ampScale(dataToFit, 1/MRSCont.fit.scale{kk}); + % Call the fit function + [fitParamsDiff1ExpMM, resBasisSetDiff1ExpMM] = fit_runFit(dataToFit, basisSetDiff1, fitModel, fitOpts); + MRSCont.fit.resBasisSet.metab{2,kk,2} = resBasisSetDiff1ExpMM; + MRSCont.fit.results.metab.fitParams{2,kk,2} = fitParamsDiff1ExpMM; + + end + + end + end + + + %%% 3. CONCATENATED FIT %%% + % For the concatenated MEGA fit, model the DIFF1 and SUM spectra + % together. + if strcmp(fitStyle, 'Concatenated') + if ~(MRSCont.flags.didFit == 1 && MRSCont.flags.speedUp && isfield(MRSCont, 'fit') && (kk > length(MRSCont.fit.results.conc.fitParams))) || ~strcmp(MRSCont.ver.Osp,MRSCont.ver.CheckOsp) + %%% 3a. FIT CONCATENATED SPECTRUM + % Select the difference and sum spectrum to put into the + % concatenated fit + dataToFit = {op_takesubspec(MRSCont.processed.metab{kk},3) op_takesubspec(MRSCont.processed.metab{kk},4)}; + basisSetConc = MRSCont.fit.basisSet; + % Create the basis set with difference and sum basis functions + basisSetConc.fids(:,:,1) = basisSetConc.fids(:,:,3); + basisSetConc.specs(:,:,1) = basisSetConc.specs(:,:,3); + basisSetConc.fids(:,:,2) = basisSetConc.fids(:,:,4); + basisSetConc.specs(:,:,2) = basisSetConc.specs(:,:,4); + basisSetConc.fids(:,:,3:4) = []; + basisSetConc.specs(:,:,3:4) = []; + + if ~strcmp(fitOpts.coMM3, 'none') + basisSetDiff1 = MRSCont.fit.basisSet; + basisSetDiff1.fids = basisSetDiff1.fids(:,:,3); + basisSetDiff1.specs = basisSetDiff1.specs(:,:,3); + fitOpts.CrFactor = 1; + [basisSetDiff1] = osp_addDiffMMPeaks(basisSetDiff1,fitOpts); + basisSetConc.fids(:,:,1) = basisSetDiff1.fids(:,:); + basisSetConc.specs(:,:,1) = basisSetDiff1.specs(:,:); + end + + % Apply scaling factor to the data + for rr = 1:length(dataToFit) + dataToFit{rr} = op_ampScale(dataToFit{rr}, 1/MRSCont.fit.scale{kk}); + end + % Call the multi-spectrum fit function + [fitParamsConc, resBasisSetConc] = fit_runFitMultiplex(dataToFit, basisSetConc, fitModel, fitOpts); + + % Save back the basis set and fit parameters to MRSCont + MRSCont.fit.resBasisSet.conc{kk} = resBasisSetConc; + MRSCont.fit.results.conc.fitParams{kk} = fitParamsConc; + end end + % ----- LCModel wrapper fit pipeline ----- + elseif strcmpi(MRSCont.opts.fit.method, 'LCModel') + [MRSCont] = LCModelWrapper(MRSCont,kk,progressText); end end %Fit again with mean MM basis function diff --git a/fit/code/osp_fit_Quality.m b/fit/code/osp_fit_Quality.m index 781952a5..a0e3d495 100644 --- a/fit/code/osp_fit_Quality.m +++ b/fit/code/osp_fit_Quality.m @@ -78,6 +78,9 @@ end else FitSpecNamesStruct.(FitSpecNames{1}){1} = 'A'; + if MRSCont.flags.isMEGA + FitSpecNamesStruct.(FitSpecNames{1}){2} = 'diff1'; + end end @@ -140,8 +143,13 @@ end case 'LCModel' dataToPlot = MRSCont.processed.(FitSpecNames{ss}){kk}; - fitParams = MRSCont.fit.results.(FitSpecNames{ss}).fitParams{bf,kk,sf}; - + if sf ==1 + dataToPlot = op_takesubspec(dataToPlot,'A'); + else + dataToPlot = op_takesubspec(dataToPlot,'diff1'); + end + fitParams = MRSCont.fit.results.(FitSpecNames{ss}).fitParams{1,kk,sf}; + % Get the LCModel plots we previously extracted from .coord % etc. [ModelOutput] = fit_LCModelParamsToModel(fitParams); diff --git a/libraries/FID-A/inputOutput/io_loadspec_niimrs.m b/libraries/FID-A/inputOutput/io_loadspec_niimrs.m index 891881df..1730b851 100644 --- a/libraries/FID-A/inputOutput/io_loadspec_niimrs.m +++ b/libraries/FID-A/inputOutput/io_loadspec_niimrs.m @@ -279,10 +279,13 @@ end end - if length(sqzDims) > length(size(fids)) % The case for 1 average with multiple subspectra - if dims.averages < dims.subSpecs + if length(sqzDims) > length(size(fids)) + if dims.averages < dims.subSpecs % The case for 1 average with multiple subspectra dims.subSpecs = 0; end + if dims.averages < dims.extras % The case for 1 average with multiple subspectra + dims.extras = 0; + end sqzDims = {}; dimsFieldNames = fieldnames(dims); for rr = 1:length(dimsFieldNames) @@ -292,9 +295,6 @@ end end subspecs = 0; - if rawSubspecs > 0 - subspecs = rawSubspecs; - end end if length(sqzDims)==5 @@ -344,9 +344,17 @@ sz=size(fids); %Remove phase cycle for Philips data - if isfield(hdr_ext, 'Manufacturer') && (strcmp(hdr_ext.Manufacturer,'Philips') || strcmp(hdr_ext.Manufacturer,'GE')) && undoPhaseCycle + if isfield(hdr_ext, 'ProtocolName') && contains(hdr_ext.ProtocolName,'HBCD') + undoPhaseCycle = 1; + end + if isfield(hdr_ext, 'Manufacturer') && strcmp(hdr_ext.Manufacturer,'Philips') && undoPhaseCycle fids = fids .* repmat(conj(fids(1,:,:,:))./abs(fids(1,:,:,:)),[size(fids,1) 1]); end + if isfield(hdr_ext, 'Manufacturer') && strcmp(hdr_ext.Manufacturer,'GE') && undoPhaseCycle + if hdr_ext.WaterSuppressed && hdr_ext.EchoTime == 0.080 + fids(:,:,2:2:end,:) = -fids(:,:,2:2:end,:); + end + end %Compared to NIfTI MRS, FID-A needs the conjugate fids = conj(fids); diff --git a/libraries/FID-A/processingTools/op_rmempty.m b/libraries/FID-A/processingTools/op_rmempty.m index 3c14212d..be3ba553 100644 --- a/libraries/FID-A/processingTools/op_rmempty.m +++ b/libraries/FID-A/processingTools/op_rmempty.m @@ -24,6 +24,10 @@ idx_nonzeroFID = find(std(in.fids,0,1)>1e-40); nonzero_fids = in.fids(:,idx_nonzeroFID); + %reshape if still not coil combined + if length(in.sz) > 2 + nonzero_fids = reshape(nonzero_fids,in.sz); + end %re-calculate Specs using fft nonzero_specs=fftshift(fft(nonzero_fids,[],in.dims.t),in.dims.t); diff --git a/load/OspreyLoad.m b/load/OspreyLoad.m index 6d71cd63..c0d5915c 100644 --- a/load/OspreyLoad.m +++ b/load/OspreyLoad.m @@ -226,6 +226,22 @@ switch MRSCont.datatype case 'P' [MRSCont] = osp_LoadP(MRSCont); + if MRSCont.flags.hasWater + MRSCont.opts.MultipleSpectra.w = fliplr(size(MRSCont.files_w)); + if length(MRSCont.opts.MultipleSpectra.w) > 1 + tempDatasets(end+1) = MRSCont.opts.MultipleSpectra.w(2); + else + MRSCont.opts.MultipleSpectra.w(2) = 1; + end + if MRSCont.opts.MultipleSpectra.w(2) < maxDatasets + MRSCont.opts.MultipleSpectra.w = ones(1,maxDatasets); + else if maxDatasets > 1 + MRSCont.opts.MultipleSpectra.w = [1:MRSCont.opts.MultipleSpectra.w(2)]; + else + MRSCont.opts.MultipleSpectra.w = [1:MRSCont.nDatasets(1)]; + end + end + end otherwise msg = 'Data type not supported. Please contact the Osprey team (gabamrs@gmail.com).'; fprintf(msg); diff --git a/overview/OspreyOverview.m b/overview/OspreyOverview.m index 73f309f1..7fc84e1c 100755 --- a/overview/OspreyOverview.m +++ b/overview/OspreyOverview.m @@ -58,7 +58,9 @@ if MRSCont.flags.hasWater OrderNames = horzcat(OrderNames, 'w'); - OrderNamesFit = horzcat(OrderNamesFit, 'w'); + if ~strcmp(MRSCont.opts.fit.method,'LCModel') + OrderNamesFit = horzcat(OrderNamesFit, 'w'); + end end S = orderfields(MRSCont.processed,OrderNames); dataPlotNames = fieldnames(S)'; @@ -90,6 +92,9 @@ end else FitSpecNamesStruct.(FitSpecNames{1}){1} = 'A'; + if MRSCont.flags.isMEGA + FitSpecNamesStruct.(FitSpecNames{1}){2} = 'diff1'; + end end MRSCont.overview.FitSpecNamesStruct = FitSpecNamesStruct; end @@ -490,6 +495,16 @@ end MRSCont.overview.Osprey.(['all_models_voxel_' num2str(rr)]).(ModelCombs{sc}){bf,kk,sf}.ppm = ppmRangeData'; MRSCont.overview.Osprey.(['all_models_voxel_' num2str(rr)]).(ModelCombs{sc}){bf,kk,sf}.res = MRSCont.overview.Osprey.(['all_models_voxel_' num2str(rr)]).(ModelCombs{sc}){bf,kk,sf}.data-MRSCont.overview.Osprey.(['all_models_voxel_' num2str(rr)]).(ModelCombs{sc}){bf,kk,sf}.fit; + if ~isempty(MRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})) + MRSCont.overview.Osprey.(['all_models_voxel_' num2str(rr)]).(ModelCombs{sc}){bf,kk,sf}.fit(ppmRangeData>MRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})(1) & ppmRangeDataMRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})(1) & ppmRangeDataMRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})(1) & ppmRangeDataMRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})(1) & ppmRangeDataMRSCont.opts.fit.GAP.(FitSpecNamesStruct.metab{sf})(1) & ppmRangeData max_diffBA diff --git a/process/osp_saveLCM.m b/process/osp_saveLCM.m index 92c2f69c..f2c43fd5 100644 --- a/process/osp_saveLCM.m +++ b/process/osp_saveLCM.m @@ -143,8 +143,14 @@ % Check if short-TE water scans exist, if so, write LCM .RAW file if MRSCont.flags.hasWater % Get TE and the input file name - te_w = MRSCont.processed.w{kk}.te; - [path_w,filename_w,~] = fileparts(MRSCont.files_w{kk}); + if strcmpi(MRSCont.vendor, 'GE') || strcmp(MRSCont.datatype,'DATA') + te_w = MRSCont.processed.w{kk}.te; + [path_w,filename_w,~] = fileparts(MRSCont.files{kk}); + else + te_w = MRSCont.processed.w{kk}.te; + [path_w,filename_w,~] = fileparts(MRSCont.files_w{kk}); + end + % For batch analysis, get the last two sub-folders (e.g. site and % subject) to augment the filename, avoiding duplicate output filenames diff --git a/quantify/OspreyQuantify.m b/quantify/OspreyQuantify.m index e5b04486..5c86add6 100644 --- a/quantify/OspreyQuantify.m +++ b/quantify/OspreyQuantify.m @@ -128,11 +128,24 @@ end end else - MRSCont.quantify.names.metab{1,1} = MRSCont.fit.results.metab.fitParams{1, 1}.name; - MRSCont.quantify.names.SubSpectra{1,1} = MRSCont.processed.metab{1, 1}.names{1}; + MRSCont.quantify.names.metab{1,ss} = MRSCont.fit.results.metab.fitParams{1, 1, ss}.name; + common_LCModelNames = {'PCh_GPC','Cr_PCr','NAA_NAAG','Glu_Gln','GABA_MM30','GABA_MM09'}; + Osprey_Names = {'tCho','tCr','tNAA','Glx','GABAplus','GABAplus'}; + for nn = 1 : length(common_LCModelNames) + idx = find(strcmp(MRSCont.quantify.names.metab{1,ss},common_LCModelNames{nn})); + if ~isempty(idx) + MRSCont.quantify.names.metab{1,ss}{idx} = Osprey_Names{nn}; + end + end + if ss == 1 + MRSCont.quantify.names.SubSpectra{1,ss} = MRSCont.processed.metab{1, 1}.names{1}; + else + MRSCont.quantify.names.SubSpectra{1,ss} = MRSCont.processed.metab{1, 1}.names{3}; + end end end + for ss = 1 : SubSpectraFitted for kk = 1:MRSCont.nDatasets(1) for mm = 1 : BasisSetsFitted @@ -149,20 +162,22 @@ end if strcmp(MRSCont.opts.fit.method, 'LCModel') - for kk = 1:MRSCont.nDatasets(1) - if ~iscell(MRSCont.fit.results) %Is SVS %Is SVS - MRSCont.quantify.CRLB{kk}.metab = MRSCont.fit.results.metab.fitParams{kk}.CRLB'; - - if qtfyH2O - MRSCont.quantify.h2oarea{kk}.metab = MRSCont.fit.results.metab.fitParams{kk}.h2oarea; - end - else %Is DualVoxel - MRSCont.quantify.CRLB{kk}.metab(:,1) = MRSCont.fit.results{1}.metab.fitParams{kk}.CRLB'; - MRSCont.quantify.CRLB{kk}.metab(:,2) = MRSCont.fit.results{2}.metab.fitParams{kk}.CRLB'; - - if qtfyH2O - MRSCont.quantify.h2oarea{kk}.metab(:,1) = MRSCont.fit.results{1}.metab.fitParams{kk}.h2oarea; - MRSCont.quantify.h2oarea{kk}.metab(:,2) = MRSCont.fit.results{2}.metab.fitParams{kk}.h2oarea; + for ss = 1 : SubSpectraFitted + for kk = 1:MRSCont.nDatasets(1) + if ~iscell(MRSCont.fit.results) %Is SVS %Is SVS + MRSCont.quantify.CRLB{1,kk,ss}.metab = MRSCont.fit.results.metab.fitParams{1,kk,ss}.CRLB'; + + if qtfyH2O + MRSCont.quantify.h2oarea{1,kk,ss}.metab = MRSCont.fit.results.metab.fitParams{1,kk,ss}.h2oarea; + end + else %Is DualVoxel + MRSCont.quantify.CRLB{kk}.metab(:,1) = MRSCont.fit.results{1}.metab.fitParams{1,kk,ss}.CRLB'; + MRSCont.quantify.CRLB{kk}.metab(:,2) = MRSCont.fit.results{2}.metab.fitParams{1,kk,ss}.CRLB'; + + if qtfyH2O + MRSCont.quantify.h2oarea{kk}.metab(:,1) = MRSCont.fit.results{1}.metab.fitParams{1,kk,ss}.h2oarea; + MRSCont.quantify.h2oarea{kk}.metab(:,2) = MRSCont.fit.results{2}.metab.fitParams{1,kk,ss}.h2oarea; + end end end end @@ -221,20 +236,32 @@ else % Get WCONC, ATTMET, and ATTH2O from control file LCMparam = osp_readlcm_control(MRSCont.opts.fit.lcmodel.controlfileA{kk}); - if isfield(LCMparam, 'WCONC') % User-supplied WCONC - amplWater = str2double(LCMparam.WCONC); + if isfield(LCMparam, 'WCONC') || isfield(LCMparam, 'wconc') % User-supplied WCONC + try + amplWater = str2double(LCMparam.WCONC); + catch + amplWater = str2double(LCMparam.wconc); + end else amplWater = 35880; %LCModel default WCONC assumes pure WM end - if isfield(LCMparam, 'ATTMET') % User-supplied ATTMET - amplWater = amplWater * str2double(LCMparam.ATTMET); + if isfield(LCMparam, 'ATTMET') || isfield(LCMparam, 'attmet')% User-supplied ATTMET + try + amplWater = amplWater * str2double(LCMparam.ATTMET); + catch + amplWater = amplWater * str2double(LCMparam.attmet); + end else %Do nothing as LCModel default ATTMET is 1 end - if isfield(LCMparam, 'ATTH2O') % User-supplied ATTH2O - amplWater = amplWater / str2double(LCMparam.ATTH2O); + if isfield(LCMparam, 'ATTH2O') || isfield(LCMparam, 'atth2o') % User-supplied ATTH2O + try + amplWater = amplWater / str2double(LCMparam.ATTH2O); + catch + amplWater = amplWater / str2double(LCMparam.atth2o); + end else amplWater = amplWater / 0.7; %LCModel default ATTH2O is 0.7 end @@ -804,9 +831,7 @@ for ss = 1 : size(amplMets,3) for AlphaMets = 1 : length(metabNames) idx = find(strcmp(metsName.metab{mm,ss},metabNames{AlphaMets})); - if isempty(idx) - idx = find(strcmp(metsName.metab{mm,ss},'Glu_Gln')); - end + [T1_Metab_GM, T1_Metab_WM, T2_Metab_GM, T2_Metab_WM] = lookUpRelaxTimes(metsName.metab{mm,ss}{idx},Bo); % average across GM and WM T1_Metab = mean([T1_Metab_GM T1_Metab_WM]); @@ -979,11 +1004,7 @@ conc = zeros(MRSCont.nDatasets(1),length(names)); for kk = 1:MRSCont.nDatasets(1) - if (strcmp(qtfyType, 'h2oarea') || strcmp(qtfyType, 'CRLB')) - conc(kk,:) = MRSCont.quantify.(qtfyType){kk}.metab(:,rr); - else - conc(kk,:) = MRSCont.quantify.(qtfyType){mm,kk,ss}.metab(:,rr)'; - end + conc(kk,:) = MRSCont.quantify.(qtfyType){mm,kk,ss}.metab(:,rr)'; end % Save back to Osprey data container if isfield(MRSCont, 'exclude')