diff --git a/GUI/Osprey.m b/GUI/Osprey.m index fa19fa20..88f9d034 100644 --- a/GUI/Osprey.m +++ b/GUI/Osprey.m @@ -15,6 +15,9 @@ % % HISTORY: % 2019-07-11: First version of the code. +%% Get version +VersionStruct = getCurrentVersion; + %% Check for available add-ons [~,~] = osp_Toolbox_Check ('OspreyGUI',0); %% Set up Layout @@ -30,7 +33,7 @@ logoFcn = @()imread('osprey.png', 'BackgroundColor', gui.colormap.Background); logoBanner = uiw.utility.loadIcon(logoFcn); % Here the intro banner is created -gui.d = uiw.dialog.About('Name', 'Osprey','Version','2.6.0','Date', 'July 15, 2024',... +gui.d = uiw.dialog.About('Name', 'Osprey','Version',VersionStruct.Version,'Date', VersionStruct.Date,... 'Timeout', 3,'CustomText', 'Osprey is provided by Johns Hopkins University.',... 'ContactInfo', 'gabamrs@gmail.com','LogoCData', logoBanner); diff --git a/GUI/osp_Toolbox_Check.m b/GUI/osp_Toolbox_Check.m index 0e0d6144..17bc60c3 100644 --- a/GUI/osp_Toolbox_Check.m +++ b/GUI/osp_Toolbox_Check.m @@ -31,8 +31,8 @@ % HISTORY: % 2020-05-15: First version of the code. %% % 1. SAVE OSPREY VERSION%%% -%%% 1. SAVE OSPREY VERSION%%% -OspreyVersion = 'Osprey 2.6.0'; +VersionStruct = getCurrentVersion; +OspreyVersion = ['Osprey ' VersionStruct.Version]; fprintf(['Timestamp %s ' OspreyVersion ' ' Module '\n'], datestr(now,'mmmm dd, yyyy HH:MM:SS')); hasSPM = 1; % For the compiled GUI %% % 2. GET SPMPATH AND TOOLBOXES%%% diff --git a/docs_source/05-tutorial_cmd.Rmd b/docs_source/05-tutorial_cmd.Rmd index c7fe999f..e50b139a 100644 --- a/docs_source/05-tutorial_cmd.Rmd +++ b/docs_source/05-tutorial_cmd.Rmd @@ -14,7 +14,7 @@ MRSCont = OspreyJob('Path/To/Job/File'); MRSCont will now contain the options and file locations for your data, as specified in the job file. Next, the data can be loaded: ```octave -MRSCont = OspreyLoad(MRScont); +MRSCont = OspreyLoad(MRSCont); ``` The substructure, MRSCont.raw, now contains the raw MRS data, as extracted from their native format. We may now preprocess these data: @@ -56,4 +56,4 @@ At any point during a command-line analysis, it is possible to visualize the dat ```octave OspreyGUI(MRSCont); -``` \ No newline at end of file +``` diff --git a/job/OspreyJob.m b/job/OspreyJob.m index b78b0190..8aedca03 100755 --- a/job/OspreyJob.m +++ b/job/OspreyJob.m @@ -44,6 +44,7 @@ %%% 1. INITIALISE DATA CONTAINER WITH DEFAULT SETTINGS [MRSCont] = OspreySettings; +VersionStruct = getCurrentVersion; %%% 2. CHECK JOB INPUT FILE FORMAT %%% [~,~,ext] = fileparts(jobFile); @@ -210,6 +211,16 @@ fprintf('Adding macromolecule and lipid basis functions to the fit (default). Please indicate otherwise in the csv-file or the GUI \n'); MRSCont.opts.fit.fitMM = 1; end + if isfield(jobStruct,'RelaxationAtlas') + MRSCont.opts.quantify.RelaxationAtlas = jobStruct(1).RelaxationAtlas; + else + MRSCont.opts.quantify.RelaxationAtlas = 0; + end + if isfield(jobStruct,'RelaxationAtlasAge') + MRSCont.opts.quantify.RelaxationAtlasAge = jobStruct(1).RelaxationAtlasAge; + else + MRSCont.opts.quantify.RelaxationAtlasAge = 0; + end if isfield(jobStruct,'deface') MRSCont.opts.img.deface = jobStruct.deface; else @@ -397,6 +408,16 @@ else MRSCont.opts.fit.fitMM = 1; end + if isfield(jobStruct,'RelaxationAtlas') + MRSCont.opts.quantify.RelaxationAtlas = jobStruct.RelaxationAtlas; + else + MRSCont.opts.quantify.RelaxationAtlas = 0; + end + if isfield(jobStruct,'RelaxationAtlasAge') + MRSCont.opts.quantify.RelaxationAtlasAge = jobStruct.RelaxationAtlasAge; + else + MRSCont.opts.quantify.RelaxationAtlasAge = 0; + end if isfield(jobStruct,'basisSet') if isfolder(jobStruct.basisSet) opts.fit.basissetFolder = jobStruct.basisSet; @@ -781,7 +802,7 @@ %%% 7. SET FLAGS AND VERSION %%% MRSCont.flags.didJob = 1; MRSCont.loadedJob = jobFile; -MRSCont.ver.Osp = 'Osprey 2.6.0'; +MRSCont.ver.Osp = ['Osprey ' VersionStruct.Version]; %%% 8. CHECK IF OUTPUT STRUCTURE ALREADY EXISTS IN OUTPUT FOLDER %%% diff --git a/libraries/FID-A/inputOutput/io_loadspec_dicom.m b/libraries/FID-A/inputOutput/io_loadspec_dicom.m index 5268af5f..e0dbf862 100644 --- a/libraries/FID-A/inputOutput/io_loadspec_dicom.m +++ b/libraries/FID-A/inputOutput/io_loadspec_dicom.m @@ -3,35 +3,97 @@ %Georg Oeltzschner, Johns Hopkins University 2019. % % USAGE: -% out=io_loadspec_dicom(folder); -% +% out=io_loadspec_dicom(source); +% % DESCRIPTION: % Loads a DICOM (.dcm, .ima) file into matlab structure format. -% +% % INPUTS: -% filename = Name of a folder with DICOM (.dcm, .ima) data to load. +% source = Name of a folder (or file) with DICOM (.dcm, .ima) data to load. % % OUTPUTS: % out = Input dataset in FID-A structure format. -function out=io_loadspec_dicom(folder); +function out=io_loadspec_dicom(source); -% If path to a .dcm is provided, then extract the folder location: -if ~isfolder(folder) && isfile(folder) - folder =fileparts(folder); +% ARC170724 : Source may be a folder, or an individual file. We later determine +% how this should be interpreted. +if ~isfolder(source) && isfile(source) + [folder,filename,ext] =fileparts(source); + requested_file=[filename ext]; +else + folder=source; + requested_file=[]; end % Create list of complete filenames (incl. path) in the folder dirFolder = dir(folder); filesInFolder = dirFolder(~[dirFolder.isdir]); hidden = logical(ones(1,length(filesInFolder))); -for jj = 1:length(filesInFolder) +for jj = 1:length(filesInFolder) if strcmp(filesInFolder(jj).name(1),'.') hidden(jj) = 0; end end -filesInFolder = filesInFolder(hidden);%delete hidden files -filesInFolder = fullfile(folder, {filesInFolder.name}); +filesInFolder = filesInFolder(hidden);%delete hidden files + +if length(filesInFolder)==0 + % No files survived. This will not end well. + error([ 'io_loadspec_dicom did not find any files in ' folder ]) +end + +requested_file_ix=find(strcmp(requested_file,{filesInFolder.name})); % NB, may be empty + +if isfolder(folder) + filesInFolder = fullfile(folder, {filesInFolder.name}); +else + % this can happen if the function is called with a wildcard, eg /path/to/abc*dcm + filesInFolder = fullfile(dirname(folder), {filesInFolder.name}); +end + +% ARC170724 : If a single file is specified (rather than a directory), check +% whether this file should be considered in isolation (pre-averaged data) or +% if it should be combined with other items in the folder. We do this by +% inspecting the Series Number header. + +if ~isempty(requested_file_ix) + % check header of the first, last and requested items; if Series Number + % does not match, don't try to combine + check_header_ix = unique([1, requested_file_ix, length(filesInFolder) ]); + encountered_series_numbers = []; + encountered_problems = []; + + for ix = check_header_ix + try + dh = dicominfo(filesInFolder{ix}); + encountered_series_numbers(end+1) = dh.SeriesNumber; + catch + encountered_problems(end+1) = ix; + end + end + + encountered_series_numbers=unique(encountered_series_numbers); + + if any(encountered_problems == requested_file_ix) + + warning([ 'io_loadspec_dicom could not determine SeriesNumber of the requested file: ' source ]); + % in this case, fall back on default behaviour + + elseif length(encountered_series_numbers)>1 || length(encountered_problems)>0 + + % if we found more than one series, OR we found some unreadable files + % (probably non-DICOM), then just retain that one specific file the + % caller originally asked for + + if length(encountered_problems)>0 + warning([ 'io_loadspec_dicom encountered some non-DICOM data in ' source ]); + end + + filesInFolder={source}; + end + + % in any other scenario, full back on the default behaviour +end % Get the header of the first file to make some decisions. DicomHeader = read_dcm_header(filesInFolder{1}); @@ -41,7 +103,7 @@ seqorig = DicomHeader.seqorig; else % deal with missing seqorig field for dicom data load on Siemens Minnesota -% sequences (Auerbach and Deelchand versions, VB17A & VE11C) +% sequences (Auerbach and Deelchand versions, VB17A & VE11C) if contains(DicomHeader.sequenceFileName, 'slaser_dkd') %Deelchand/Oz seqorig = 'CMRR'; elseif contains(DicomHeader.sequenceFileName, 'eja_svs') %Auerbach/Marjanska @@ -80,7 +142,7 @@ end % Collect all FIDs and sort them into fids array for kk = 1:length(filesInFolder) - + % First, attempt to retrieve the FID from the DICOM header: infoDicom = dicominfo(filesInFolder{kk}); if isfield(infoDicom, 'SpectroscopyData') @@ -95,7 +157,7 @@ fids(:,kk) = dicom_get_spectrum_siemens(fd); fclose(fd); end - + end @@ -115,16 +177,16 @@ else out.flags.averaged = 0; end - + % Currently, the DICOM recon of the universal sequence is flawed. % Kick out empty lines here and see if data can be reconstructed. if strcmp(seqorig, 'Universal') && out.flags.averaged == 0 fids(:,size(fids,2)/2+1:end) = []; end - + % Rearrange into subspecs fids = reshape(fids,[size(fids,1) size(fids,2)/2 2]); - + elseif strcmp(seqtype,'PRESS') || strcmp(seqtype,'STEAM') || strcmp(seqtype,'sLASER') % If the number of stored FIDs does not match the number of averages % stored in the DICOM header, the data are averaged. @@ -137,9 +199,9 @@ elseif strcmp(seqtype,'HERMES') - out.flags.averaged = 0; + out.flags.averaged = 0; % Currently, the DICOM recon of the universal sequence is flawed. - % Kick out empty lines here and see if data can be reconstructed. + % Kick out empty lines here and see if data can be reconstructed. if strcmp(seqorig, 'Universal') && out.flags.averaged == 0 fids(:,size(fids,2)/2+1:end) = []; end @@ -222,7 +284,7 @@ %Find the number of averages. 'averages' will specify the current number %of averages in the dataset as it is processed, which may be subject to -%change. 'rawAverages' will specify the original number of acquired +%change. 'rawAverages' will specify the original number of acquired %averages in the dataset, which is unchangeable. if dims.subSpecs ~=0 if dims.averages~=0 @@ -244,7 +306,7 @@ %Find the number of subspecs. 'subspecs' will specify the current number %of subspectra in the dataset as it is processed, which may be subject to -%change. 'rawSubspecs' will specify the original number of acquired +%change. 'rawSubspecs' will specify the original number of acquired %subspectra in the dataset, which is unchangeable. if dims.subSpecs ~=0 subspecs=sz(dims.subSpecs); @@ -288,8 +350,8 @@ out.fids=fids; out.specs=specs; out.sz=sz; -out.ppm=ppm; -out.t=t; +out.ppm=ppm; +out.t=t; out.spectralwidth=spectralwidth; out.dwelltime=dwelltime; out.txfrq=txfrq; diff --git a/libraries/FID-A/inputOutput/io_loadspec_twix.m b/libraries/FID-A/inputOutput/io_loadspec_twix.m index e57d6e64..375b5e4f 100644 --- a/libraries/FID-A/inputOutput/io_loadspec_twix.m +++ b/libraries/FID-A/inputOutput/io_loadspec_twix.m @@ -77,7 +77,8 @@ isjnseq=~isempty(strfind(sequence,'jn_')); %Is this another one of Jamie Near's sequences? isWIP529=~isempty(strfind(sequence,'edit_529'));%Is this WIP 529 (MEGA-PRESS)? ismodWIP=((~isempty(strfind(sequence,'\svs_edit'))||... - ~isempty(strfind(sequence,'\wip_svs_edit'))) && isempty(strfind(sequence,'edit_859'))); %Modified WIP + ~isempty(strfind(sequence,'\wip_svs_edit'))) && ... + (isempty(strfind(sequence,'edit_859')) && isempty(strfind(sequence,'univ')))); %Modified WIP isWIP859=~isempty(strfind(sequence,'edit_859'))||...;%Is this WIP 859 (MEGA-PRESS)? ~isempty(strfind(sequence,'WIP_859'));%Is this WIP 859 (MEGA-PRESS)? Other naming isTLFrei=~isempty(strfind(sequence,'md_svs_edit')) ||... %Is Thomas Lange's MEGA-PRESS sequence @@ -353,9 +354,8 @@ else % GO 6/2024: Encountered a VD13 Siemens product sequence that had % transients in Set, not in Ave... - if strcmp(sqzDims,'Ave') - dims.averages=find(strcmp(sqzDims,'Ave')); - else + dims.averages=find(strcmp(sqzDims,'Ave')); + if isempty(dims.averages) dims.averages=find(strcmp(sqzDims,'Set')); end end diff --git a/process/OspreyProcess.m b/process/OspreyProcess.m index 810364a9..992a5714 100644 --- a/process/OspreyProcess.m +++ b/process/OspreyProcess.m @@ -745,18 +745,18 @@ end end if MRSCont.flags.hasRef - MRSCont.QM.SNR.ref(ref_ll,kk) = op_getSNR(MRSCont.processed.ref{kk}); - MRSCont.QM.FWHM.ref(ref_ll,kk) = op_getLW(MRSCont.processed.ref{kk},4.2,5.2); + MRSCont.QM.SNR.ref(ref_ll,kk) = op_getSNR(MRSCont.processed.ref{kk},3.7,5.7); + MRSCont.QM.FWHM.ref(ref_ll,kk) = op_getLW(MRSCont.processed.ref{kk},3.7,5.7); MRSCont.processed.ref{kk}.QC_names = {'water'}; end if MRSCont.flags.hasWater - MRSCont.QM.SNR.w(w_ll,kk) = op_getSNR(MRSCont.processed.w{kk}); - MRSCont.QM.FWHM.w(w_ll,kk) = op_getLW(MRSCont.processed.w{kk},4.2,5.2); + MRSCont.QM.SNR.w(w_ll,kk) = op_getSNR(MRSCont.processed.w{kk},3.7,5.7); + MRSCont.QM.FWHM.w(w_ll,kk) = op_getLW(MRSCont.processed.w{kk},3.7,5.7); MRSCont.processed.w{kk}.QC_names = {'water'}; end if MRSCont.flags.hasMMRef - MRSCont.QM.SNR.mm_ref(mm_ref_ll,kk) = op_getSNR(MRSCont.processed.mm_ref{kk}); - MRSCont.QM.FWHM.mm_ref(mm_ref_ll,kk) = op_getLW(MRSCont.processed.mm_ref{kk},4.2,5.2); + MRSCont.QM.SNR.mm_ref(mm_ref_ll,kk) = op_getSNR(MRSCont.processed.mm_ref{kk},3.7,5.7); + MRSCont.QM.FWHM.mm_ref(mm_ref_ll,kk) = op_getLW(MRSCont.processed.mm_ref{kk},3.7,5.7); MRSCont.processed.mm_ref{kk}.QC_names = {'water'}; end end diff --git a/quantify/OspreyQuantify.m b/quantify/OspreyQuantify.m index 5c86add6..8c3f03bd 100644 --- a/quantify/OspreyQuantify.m +++ b/quantify/OspreyQuantify.m @@ -112,6 +112,12 @@ mkdir(saveDestination); end +% Handle path to Atlas-based resources for relaxation correction +if isfield(MRSCont.opts, 'quantify') && isfield(MRSCont.opts.quantify, 'RelaxationAtlas') && MRSCont.opts.quantify.RelaxationAtlas && ~(ismcc || isdeployed) + AtlasPath = fileparts(which('quantify/atlas/AtlasLookUpTable.mat')); +elseif isfield(MRSCont.opts, 'quantify') && isfield(MRSCont.opts.quantify, 'RelaxationAtlas') && MRSCont.opts.quantify.RelaxationAtlas && (ismcc || isdeployed) + AtlasPath = spm_select(1,'dir','Select the folder that contains the relaxometry atlas resources ("quantify/atlas")',{},pwd); +end % Add combinations of metabolites to the basisset for ss = 1 : SubSpectraFitted @@ -312,19 +318,51 @@ % Apply tissue correction fGM = MRSCont.seg.tissue.fGM(kk,:); fWM = MRSCont.seg.tissue.fWM(kk,:); - TissCorrWaterScaled = quantTiss(metsName, amplMets, amplWater, metsTR, waterTR, metsTE, waterTE, fGM, fWM, fCSF,Bo); - % Save back to Osprey data container - for ss = 1 :SubSpectraFitted - for mm = 1 :BasisSetsFitted - if ~isempty( TissCorrWaterScaled{mm,ss}) - MRSCont.quantify.metab.TissCorrWaterScaled{mm,kk,ss} = TissCorrWaterScaled{mm,ss}.metab; + if ~isfield(MRSCont.opts,'quantify') || ~isfield(MRSCont.opts.quantify, 'RelaxationAtlas') || ~MRSCont.opts.quantify.RelaxationAtlas + TissCorrWaterScaled = quantTiss(metsName, amplMets, amplWater, metsTR, waterTR, metsTE, waterTE, fGM, fWM, fCSF,Bo); + % Save back to Osprey data container + for ss = 1 :SubSpectraFitted + for mm = 1 :BasisSetsFitted + if ~isempty( TissCorrWaterScaled{mm,ss}) + MRSCont.quantify.metab.TissCorrWaterScaled{mm,kk,ss} = TissCorrWaterScaled{mm,ss}.metab; + end end end + else % Use relaxation atlas for region- (& age-)specifc water relaxation reference values + if isfield(MRSCont.opts.quantify, 'RelaxationAtlas') && MRSCont.opts.quantify.RelaxationAtlas + if isfield(MRSCont.opts.quantify, 'RelaxationAtlasAge') && MRSCont.opts.quantify.RelaxationAtlasAge % If age-dependent correction enabled, do it + if MRSCont.flags.hasStatfile % Load statfile: + statFile = readtable(MRSCont.file_stat, 'Delimiter', ',','ReadVariableNames',1); % Load CSV input + name = statFile.Properties.VariableNames; + if ~any(matches(name,{'age'})) % Check 'age' is present + error('No "age" column in the stat file. This is required for using the age-dependent atlas!'); + else + Age = statFile.age(kk); + end + else + error('Missing stat file!') + end + else % If age-dependent correction not enabled, then make dummy variables: + Age = []; + MRSCont.opts.quantify.RelaxationAtlasAge = 0; + end + + % Run correction: + [TissCorrWaterScaled,AtlasRelax] = quantTissAtlas(metsName, amplMets, amplWater, metsTR, waterTR, metsTE, waterTE, fGM, fWM, fCSF,Bo, MRSCont.coreg.vol_mask{kk}.fname,AtlasPath,Age,MRSCont.opts.quantify.RelaxationAtlasAge); + + for ss = 1 :SubSpectraFitted + for mm = 1 :BasisSetsFitted + if ~isempty( TissCorrWaterScaled{mm,ss}) + MRSCont.quantify.metab.TissCorrWaterScaled{mm,kk,ss} = TissCorrWaterScaled{mm,ss}.metab; + end + end + end + MRSCont.quantify.metab.AtlasT1T2(kk,1:4) = AtlasRelax; % Save relaxation times (per dataset) back in the container [T1_GM T1_WM T2_GM T2_WM] + end end end - %%% 6. GET ALPHA CORRECTION (THIS IS FOR GABA ONLY AT THIS POINT) %%% if qtfyAlpha @@ -766,6 +804,169 @@ %%% /Calculate CSF-corrected water-scaled estimates %%% +%%% Calculate tissue-corrected water-scaled estimates using the relaxometry atlas %%% +function [TissCorrWaterScaled, AtlasRelax] = quantTissAtlas(metsName, amplMets, amplWater, metsTR, waterTR, metsTE, waterTE, fGM, fWM, fCSF, Bo, CoregPath, AtlasPath, Age, Age_correct) +% This function calculates water-scaled, tissue-corrected metabolite +% estimates in molal units, according to Gasparovic et al, Magn Reson Med +% 55:1219-26 (2006). + +% Additionally it performs atlas-based regional (and age) dependent +% correction: Simegn et al 2024 + +% Define Constants +switch Bo + case '3T' + atlas_filename = fullfile(AtlasPath,'atlas_130.nii'); % EVE ATLAS, symmetrizwed by L/R parcels, in SPM space + ATLASvol = spm_vol(atlas_filename); + lookup_filename = fullfile(AtlasPath,'AtlasLookUpTable.mat'); % Load the lookup table for parcel indices + load(lookup_filename) + + %load voxel mask + [Path,FName,Ext] = fileparts(CoregPath); + voxel_mask_filename = CoregPath; + Rvoxel_mask_filename = fullfile(Path,['r',FName,Ext]); + + if(isfile(Rvoxel_mask_filename)~=1) + % Need to reslice the voxel mask to the atlas raster + if ~exist(voxel_mask_filename) && exist([voxel_mask_filename,'.gz']) + gunzip([voxel_mask_filename,'.gz']); + end + clear matlabbatch + matlabbatch{1}.spm.spatial.coreg.write.ref = {atlas_filename}; + matlabbatch{1}.spm.spatial.coreg.write.source = {voxel_mask_filename}; + matlabbatch{1}.spm.spatial.coreg.write.roptions.interp = 0; + matlabbatch{1}.spm.spatial.coreg.write.roptions.wrap = [0 0 0]; + matlabbatch{1}.spm.spatial.coreg.write.roptions.mask = 0; + matlabbatch{1}.spm.spatial.coreg.write.roptions.prefix = 'r'; + spm_jobman('run',matlabbatch); + end + + %load resliced voxel mask + vol_mask = spm_vol(Rvoxel_mask_filename); + Parcels_vol = ATLASvol(1).private.dat(:,:,:,1) .* vol_mask(1).private.dat(:,:,:); + + %count the unique parcel elements in the voxel (ignoring zeros) + unique_parcels = unique(Parcels_vol(:)); + unique_parcels = round(unique_parcels(unique_parcels~=0)); %contains parcel indices for parcels in MRS voxel + n_unique_parcels = histc(Parcels_vol(:),unique_parcels); + n_unique_parcels = n_unique_parcels(unique_parcels~=0); %contains voxel-counts for parcels in MRS voxel + + %Extract T1 and T2 for each parcel + T1_voxel_mask_parcels=zeros(size(unique_parcels)); + T2_voxel_mask_parcels=zeros(size(unique_parcels)); + GM =zeros(size(unique_parcels)); + WM =zeros(size(unique_parcels)); + for ii=1:length(unique_parcels) + if(Age_correct) + T1_voxel_mask_parcels(ii) = AtlasLookUpTable.T1(unique_parcels(ii)) + (Age-30)*AtlasLookUpTable.T1_slope(unique_parcels(ii)); + T2_voxel_mask_parcels(ii) = AtlasLookUpTable.T2(unique_parcels(ii)) + (Age-30)*AtlasLookUpTable.T2_slope(unique_parcels(ii)); + else + T1_voxel_mask_parcels(ii) = AtlasLookUpTable.T1(unique_parcels(ii)); + T2_voxel_mask_parcels(ii) = AtlasLookUpTable.T2(unique_parcels(ii)); + end + GM(ii) = AtlasLookUpTable.GM(unique_parcels(ii)); + WM(ii) = AtlasLookUpTable.WM(unique_parcels(ii)); + end + %Calculate weighted mean of T1 and T2 for GM and WM + T1w_GM= sum(T1_voxel_mask_parcels(GM==1).*n_unique_parcels(GM==1))./sum(n_unique_parcels(GM==1)); + T2w_GM= sum(T2_voxel_mask_parcels(GM==1).*n_unique_parcels(GM==1))./sum(n_unique_parcels(GM==1)); + T1w_WM= sum(T1_voxel_mask_parcels(WM==1).*n_unique_parcels(WM==1))./sum(n_unique_parcels(WM==1)); + T2w_WM= sum(T2_voxel_mask_parcels(WM==1).*n_unique_parcels(WM==1))./sum(n_unique_parcels(WM==1)); + AtlasRelax = [T1w_GM T1w_WM T2w_GM T2w_WM]; + T1w_CSF = 3.817; + T2w_CSF = 0.503; + case '7T' + error('No atlas for 7T yet!') +end + +% Determine concentration of water in GM, WM and CSF +% Gasparovic et al. 2006 (MRM) uses relative densities, ref to +% Ernst et al. 1993 (JMR) +% fGM = 0.78 +% fWM = 0.65 +% fCSF = 0.97 +% such that +% concw_GM = 0.78 * 55.51 mol/kg = 43.30 +% concw_WM = 0.65 * 55.51 mol/kg = 36.08 +% concw_CSF = 0.97 * 55.51 mol/kg = 53.84 +concW_GM = 43.30*1e3; +concW_WM = 36.08*1e3; +concW_CSF = 53.84*1e3; +molal_concW = 55.51*1e3; + +% Gasparovic et al. method +% Calculate molal fractions from volume fractions (equivalent to eqs. 5-7 in Gasparovic et al., 2006) +molal_fGM = (fGM*concW_GM) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); +molal_fWM = (fWM*concW_WM) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); +molal_fCSF = (fCSF*concW_CSF) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); + +for mm = 1 : size(amplMets,1) + for ss = 1 : size(amplMets,3) + % Metabolites + for kk = 1:length(metsName.metab{mm,ss}) + [T1_Metab_GM(kk), T1_Metab_WM(kk), T2_Metab_GM(kk), T2_Metab_WM(kk)] = lookUpRelaxTimes(metsName.metab{mm,ss}{kk},Bo); + % average across GM and WM + T1_Metab(kk) = mean([T1_Metab_GM(kk) T1_Metab_WM(kk)]); + T2_Metab(kk) = mean([T2_Metab_GM(kk) T2_Metab_WM(kk)]); + + % Calculate water-scaled, tissue-corrected molal concentration + % estimates + if ~isempty(amplMets{mm,1,ss}) + TissCorrWaterScaled{mm,ss}.metab(kk,:) = (amplMets{mm,ss}.metab(kk,:) ./ amplWater) .* molal_concW ... + .* (molal_fGM * (1 - exp(-waterTR/T1w_GM)) * exp(-waterTE/T2w_GM) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk))) + ... + molal_fWM * (1 - exp(-waterTR/T1w_WM)) * exp(-waterTE/T2w_WM) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk))) + ... + molal_fCSF * (1 - exp(-waterTR/T1w_CSF)) * exp(-waterTE/T2w_CSF) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk)))) ./ ... + (1 - molal_fCSF); + end + end + end +end + +% Determine concentration of water in GM, WM and CSF +% Gasparovic et al. 2006 (MRM) uses relative densities, ref to +% Ernst et al. 1993 (JMR) +% fGM = 0.78 +% fWM = 0.65 +% fCSF = 0.97 +% such that +% concw_GM = 0.78 * 55.51 mol/kg = 43.30 +% concw_WM = 0.65 * 55.51 mol/kg = 36.08 +% concw_CSF = 0.97 * 55.51 mol/kg = 53.84 +concW_GM = 43.30*1e3; +concW_WM = 36.08*1e3; +concW_CSF = 53.84*1e3; +molal_concW = 55.51*1e3; + +% Gasparovic et al. method +% Calculate molal fractions from volume fractions (equivalent to eqs. 5-7 in Gasparovic et al., 2006) +molal_fGM = (fGM*concW_GM) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); +molal_fWM = (fWM*concW_WM) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); +molal_fCSF = (fCSF*concW_CSF) ./ (fGM*concW_GM + fWM*concW_WM + fCSF*concW_CSF); + +for mm = 1 : size(amplMets,1) + for ss = 1 : size(amplMets,3) + % Metabolites + for kk = 1:length(metsName.metab{mm,ss}) + [T1_Metab_GM(kk), T1_Metab_WM(kk), T2_Metab_GM(kk), T2_Metab_WM(kk)] = lookUpRelaxTimes(metsName.metab{mm,ss}{kk},Bo); + % average across GM and WM + T1_Metab(kk) = mean([T1_Metab_GM(kk) T1_Metab_WM(kk)]); + T2_Metab(kk) = mean([T2_Metab_GM(kk) T2_Metab_WM(kk)]); + + % Calculate water-scaled, tissue-corrected molal concentration + % estimates + if ~isempty(amplMets{mm,1,ss}) + TissCorrWaterScaled{mm,ss}.metab(kk,:) = (amplMets{mm,ss}.metab(kk,:) ./ amplWater) .* molal_concW ... + .* (molal_fGM * (1 - exp(-waterTR/T1w_GM)) * exp(-waterTE/T2w_GM) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk))) + ... + molal_fWM * (1 - exp(-waterTR/T1w_WM)) * exp(-waterTE/T2w_WM) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk))) + ... + molal_fCSF * (1 - exp(-waterTR/T1w_CSF)) * exp(-waterTE/T2w_CSF) / ((1 - exp(-metsTR/T1_Metab(kk))) * exp(-metsTE/T2_Metab(kk)))) ./ ... + (1 - molal_fCSF); + end + end + end +end +end +%%%/ Calculate tissue-corrected water-scaled estimates using the relaxometry atlas %%% + %%% Calculate alpha-corrected water-scaled GABA estimates %%% function [AlphaCorrWaterScaled, AlphaCorrWaterScaledGroupNormed] = quantAlpha(metsName, amplMets, amplWater, metsTR, waterTR, metsTE, waterTE, fGM, fWM, fCSF, meanfGM, meanfWM,coMM3,Bo) @@ -962,7 +1163,6 @@ end end - %%% / Lookup function for metabolite relaxation times %%% %%% Function to create metabolite overview in MATLAB table format %%% diff --git a/quantify/atlas/AtlasLookUpTable.mat b/quantify/atlas/AtlasLookUpTable.mat new file mode 100644 index 00000000..9728c31b Binary files /dev/null and b/quantify/atlas/AtlasLookUpTable.mat differ diff --git a/quantify/atlas/atlas_130.nii b/quantify/atlas/atlas_130.nii new file mode 100644 index 00000000..2a60270a Binary files /dev/null and b/quantify/atlas/atlas_130.nii differ diff --git a/settings/version.json b/settings/version.json new file mode 100644 index 00000000..ef09467c --- /dev/null +++ b/settings/version.json @@ -0,0 +1,3 @@ +{"Version": "2.7.0", + "Date": "August 15, 2024", +"Description": "First digit for major releases, second digit for regular compiled releases, third digit for each commit to the develop branch. "} diff --git a/utilities/CompileOspreyStandalone.m b/utilities/CompileOspreyStandalone.m index 64c64b6a..1bae299a 100644 --- a/utilities/CompileOspreyStandalone.m +++ b/utilities/CompileOspreyStandalone.m @@ -70,6 +70,7 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat end end +VersionStruct = getCurrentVersion; if ~exist(OutputDir) mkdir(OutputDir); @@ -202,7 +203,7 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat 'OutputDir',OutputDirHBCD,... 'ExecutableIcon',fullfile(OspreyDir, 'graphics','osprey.gif'),... 'ExecutableSplashScreen',fullfile(OspreyDir, 'graphics','osprey.gif'),... - 'ExecutableVersion','2.6.0',... + 'ExecutableVersion',VersionStruct.Version,... 'ExecutableName','OspreyHBCD',... 'AdditionalFiles',{ fullfile(SPM12Dir),... fullfile(OspreyDir,'coreg'),... @@ -243,7 +244,7 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat 'OutputDir',OutputDirCmd,... 'ExecutableIcon',fullfile(OspreyDir, 'graphics','osprey.gif'),... 'ExecutableSplashScreen',fullfile(OspreyDir, 'graphics','osprey.gif'),... - 'ExecutableVersion','2.6.0',... + 'ExecutableVersion',VersionStruct.Version,... 'ExecutableName','OspreyCMD',... 'AdditionalFiles',{ fullfile(SPM12Dir),... fullfile(OspreyDir,'coreg'),... @@ -273,15 +274,17 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat buildResults = compiler.build.standaloneApplication (opts); -end -if PackageInstaller - compiler.package.installer(buildResults,'OutputDir',fullfile(OspreyDir,'OspreyCMD_installer'), 'InstallerName', 'OspreyCMD_install', 'RuntimeDelivery', 'installer',... - 'AuthorCompany','The Johns Hopkins University', 'Description', 'Installer for Osprey CMD version',... - 'InstallationNotes', 'Thank you for downloading Osprey. This installer will allow you to run the Osprey CMD as a standalone application. Please copy the "basissets" folder into the "application" folder after the installation.',... - 'Version', '2.6.0', 'ApplicationName','OspreyCMD') + if PackageInstaller + compiler.package.installer(buildResults,'OutputDir',fullfile(OspreyDir,'OspreyCMD_installer'), 'InstallerName', 'OspreyCMD_install', 'RuntimeDelivery', 'installer',... + 'AuthorCompany','The Johns Hopkins University', 'Description', 'Installer for Osprey CMD version',... + 'InstallationNotes', 'Thank you for downloading Osprey. This installer will allow you to run the Osprey CMD as a standalone application. Please copy the "basissets" folder into the "application" folder after the installation.',... + 'Version', VersionStruct.Version, 'ApplicationName','OspreyCMD') + end end + + %% 4. GUI export % Compile GUI Osprey which is only working if you supplied the path to the % widgets and GUI Layout Toolboxes @@ -292,7 +295,7 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat 'OutputDir',OutputDirGUI,... 'ExecutableIcon',fullfile(OspreyDir, 'graphics','osprey.gif'),... 'ExecutableSplashScreen',fullfile(OspreyDir, 'graphics','osprey.gif'),... - 'ExecutableVersion','2.6.0',... + 'ExecutableVersion',VersionStruct.Version,... 'ExecutableName','Osprey',... 'AdditionalFiles',{fullfile(WidgetsDir),... fullfile(GUILayoutDir),... @@ -329,13 +332,13 @@ function CompileOspreyStandalone(OutputDir,SPM12Dir,CreateHBCD, CreateCMD, Creat buildResults = compiler.build.standaloneApplication(opts); -end -if PackageInstaller - compiler.package.installer(buildResults,'OutputDir',fullfile(OutputDirGUI,'OspreyGUI_installer'), 'InstallerName', 'OspreyGUI_install', 'RuntimeDelivery', 'web',... - 'AuthorCompany','The Johns Hopkins University', 'Description', 'Installer for Osprey GUI version',... - 'InstallationNotes', 'Thank you for downloading Osprey. This installer will allow you to run the Osprey GUI as a standalone application. Please copy the "basissets" folder into the "application" folder after the installation',... - 'Version', '2.6.0', 'ApplicationName','OspreyGUI') + if PackageInstaller + compiler.package.installer(buildResults,'OutputDir',fullfile(OutputDirGUI,'OspreyGUI_installer'), 'InstallerName', 'OspreyGUI_install', 'RuntimeDelivery', 'web',... + 'AuthorCompany','The Johns Hopkins University', 'Description', 'Installer for Osprey GUI version',... + 'InstallationNotes', 'Thank you for downloading Osprey. This installer will allow you to run the Osprey GUI as a standalone application. Please copy the "basissets" folder into the "application" folder after the installation',... + 'Version', VersionStruct.Version, 'ApplicationName','OspreyGUI') + end end end \ No newline at end of file diff --git a/utilities/getCurrentVersion.m b/utilities/getCurrentVersion.m new file mode 100644 index 00000000..f1667f89 --- /dev/null +++ b/utilities/getCurrentVersion.m @@ -0,0 +1,38 @@ +function VersionStruct = getCurrentVersion() +%% VersionStruct = getCurrentVersion() +% This function reads the current version and release date from the +% version.json file in the settings folder +% +% USAGE: +% VersionStruct = getCurrentVersion() +% +% INPUTS: none +% +% +% OUTPUTS: +% VersionStruct = Struct with version number and release date. +% +% AUTHOR: +% Dr. Helge Zoellner (Johns Hopkins University, 2024-07-26) +% hzoelln2@jhmi.edu +% +% HISTORY: +% 2024-07-26: First version of the code. + %% Get version + fid = fopen(which(fullfile('settings','version.json'))); % Open and read + raw = fread(fid,inf); + str = char(raw'); + fclose(fid); + + if strcmp('win',osp_platform('filesys')) % Correct single backslashes in Windows paths + str = strrep(str,'\','\\'); + end + + % Following line replaces white space characters + pattern = '[ \t\n]*"'; % Match zero or more spaces, tabs, or newlines, followed by a double quote + replacement = '"'; % Replace the matched string with just a double quote + str = regexprep(str, pattern, replacement); + + % Return struct + VersionStruct = jsondecode(str); +end \ No newline at end of file diff --git a/utilities/osp_CheckRunPreviousModule.m b/utilities/osp_CheckRunPreviousModule.m index 6c5e9751..60188da8 100644 --- a/utilities/osp_CheckRunPreviousModule.m +++ b/utilities/osp_CheckRunPreviousModule.m @@ -21,8 +21,8 @@ % AUTHOR: % Georg Oeltzschner (Johns Hopkins University, 2021-07-07) % goeltzs1@jhmi.edu -% - +%% +VersionStruct = getCurrentVersion; %Get current version switch module case 'OspreyLoad' requiredModules = {'OspreyJob'}; @@ -74,7 +74,7 @@ end %Do the toolbox check here -OspreyVersion = 'Osprey 2.6.0'; +OspreyVersion = ['Osprey ' VersionStruct.Version]; hasSPM = 1; [hasSPM,OspreyVersion ] = osp_Toolbox_Check (module,MRSCont.flags.isGUI);