From f64a568069525ae8d21afb42eb89cf666088df2d Mon Sep 17 00:00:00 2001 From: Helge <47348963+HJZollner@users.noreply.github.com> Date: Mon, 23 May 2022 15:35:02 -0400 Subject: [PATCH] [FEATURE REQUEST] - Add pre/post alignment - Eric Porges - Added pre/post alignment spectra to the OspreyProcess plot --- GUI/osp_updateProWindow.m | 2 +- plot/osp_plotProcess.m | 48 +++++++++++++++++++++++++++++++-------- process/OspreyProcess.m | 19 +++++++++------- 3 files changed, 50 insertions(+), 19 deletions(-) diff --git a/GUI/osp_updateProWindow.m b/GUI/osp_updateProWindow.m index 83f25207..cac8c193 100644 --- a/GUI/osp_updateProWindow.m +++ b/GUI/osp_updateProWindow.m @@ -102,7 +102,7 @@ function osp_updateProWindow(gui) set( temp.Children(1).Children, 'Parent', gui.layout.proDrift.Children ); % Update drift plot set( gui.layout.proDrift.Children,'Children',flipud(gui.layout.proDrift.Children.Children)); set( gui.layout.proDrift.Children, 'YLim', temp.Children(1).YLim); - set( temp.Children(2).Children, 'Parent', gui.layout.proAlgn.Children ); % Update aligned and averaged plot + set( flipud(temp.Children(2).Children), 'Parent', gui.layout.proAlgn.Children ); % Update aligned and averaged plot set( gui.layout.proAlgn.Children, 'XLim', temp.Children(2).XLim); set( gui.layout.proAlgn.Children, 'YLim', temp.Children(2).YLim); set( temp.Children(3).Children, 'Parent', gui.layout.proPost.Children ); % Update post alignment plot diff --git a/plot/osp_plotProcess.m b/plot/osp_plotProcess.m index 9a1783d9..94d809ce 100755 --- a/plot/osp_plotProcess.m +++ b/plot/osp_plotProcess.m @@ -215,6 +215,12 @@ if MRSCont.flags.isHERMES || MRSCont.flags.isHERCULES || MRSCont.flags.isMEGA procDataToPlot = op_takesubspec(procDataToPlot,find(strcmp(which_sub_spec,procDataToPlot.names))); end + if isfield(MRSCont,'processed_no_align') + NoAlignProcDataToPlot = op_takeextra(MRSCont.processed_no_align.(which_spec){kk},ExtraIndex); + if MRSCont.flags.isHERMES || MRSCont.flags.isHERCULES || MRSCont.flags.isMEGA + NoAlignProcDataToPlot = op_takesubspec(NoAlignProcDataToPlot,find(strcmp(which_sub_spec,NoAlignProcDataToPlot.names))); + end + end % Get sub-spectra, depending on whether they are stored as such if MRSCont.flags.isHERMES || MRSCont.flags.isHERCULES if raw.subspecs == 4 @@ -252,6 +258,10 @@ rawDataToPlot = MRSCont.raw{ExtraIndex,kk}; procDataToPlot = op_takeextra(MRSCont.processed.(which_spec){kk},ExtraIndex); procDataToPlot = op_takesubspec(procDataToPlot,find(strcmp(which_sub_spec,procDataToPlot.names))); + if isfield(MRSCont,'processed_no_align') + NoAlignProcDataToPlot = op_takeextra(MRSCont.processed_no_align.(which_spec){kk},ExtraIndex); + NoAlignProcDataToPlot = op_takesubspec(NoAlignProcDataToPlot,find(strcmp(which_sub_spec,NoAlignProcDataToPlot.names))); + end if MRSCont.flags.isHERMES || MRSCont.flags.isHERCULES temp_spec = cat(3,rawDataToPlot.specs(:,:,procDataToPlot.commuteOrder(1)),rawDataToPlot.specs(:,:,procDataToPlot.commuteOrder(2)),... rawDataToPlot.specs(:,:,procDataToPlot.commuteOrder(3)),rawDataToPlot.specs(:,:,procDataToPlot.commuteOrder(4))); @@ -848,6 +858,24 @@ %%% 6. PLOT PROCESSED %%% % Add the data and plot hold(ax_proc, 'on'); +if exist('NoAlignProcDataToPlot','var') + if isfield(MRSCont,'plot') && isfield(MRSCont.plot, 'processed') && (MRSCont.plot.processed.match == 1) + if MRSCont.flags.hasRef || MRSCont.flags.hasWater + if MRSCont.flags.hasRef + plot(ax_proc, NoAlignProcDataToPlot.ppm, real(NoAlignProcDataToPlot.specs)/MRSCont.plot.processed.ref.max(kk), 'Color',MRSCont.colormap.LightAccent, 'LineWidth', 1); + else + plot(ax_proc, NoAlignProcDataToPlot.ppm, real(NoAlignProcDataToPlot.specs)/MRSCont.plot.processed.w.max(kk), 'Color',MRSCont.colormap.LightAccent, 'LineWidth', 1); + end + + else + plot(ax_proc, NoAlignProcDataToPlot.ppm, real(NoAlignProcDataToPlot.specs(:,1)), 'Color',MRSCont.colormap.LightAccent, 'LineWidth', 1); + end + else + plot(ax_proc, NoAlignProcDataToPlot.ppm, real(NoAlignProcDataToPlot.specs(:,1))/max(real(NoAlignProcDataToPlot.specs(NoAlignProcDataToPlot.ppm>ppmmin&NoAlignProcDataToPlot.ppm 1 - crDriftPre = MRSCont.QM{VoxelIndex}.drift.pre.(which_sub_spec){kk} + MRSCont.QM{VoxelIndex}.freqShift.(which_sub_spec)(kk)/applyDataToPlot.txfrq*1e6; - crDriftPost = MRSCont.QM{VoxelIndex}.drift.post.(which_sub_spec){kk} + MRSCont.QM{VoxelIndex}.freqShift.(which_sub_spec)(kk)/applyDataToPlot.txfrq*1e6; + crDriftPre = MRSCont.QM{VoxelIndex}.drift.pre.(which_sub_spec){kk}; + crDriftPost = MRSCont.QM{VoxelIndex}.drift.post.(which_sub_spec){kk}; hold(ax_drift, 'on'); colors = ones(length(crDriftPre),1).*colormap.Foreground; for dots = 1 : length(crDriftPre) @@ -983,8 +1011,8 @@ elseif (isfield(MRSCont.flags,'isMRSI') && (MRSCont.flags.isMRSI == 1)) if isfield(MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.pre, which_sub_spec) if length(MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.pre.(which_sub_spec){kk}) > 1 - crDriftPre = MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.pre.(which_sub_spec){kk} + MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.freqShift.(which_sub_spec)(kk)/applyDataToPlot.txfrq*1e6; - crDriftPost = MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.post.(which_sub_spec){kk} + MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.freqShift.(which_sub_spec)(kk)/applyDataToPlot.txfrq*1e6; + crDriftPre = MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.pre.(which_sub_spec){kk}; + crDriftPost = MRSCont.QM{VoxelIndex(1), VoxelIndex(2)}.drift.post.(which_sub_spec){kk}; hold(ax_drift, 'on'); colors = ones(length(crDriftPre),1).*colormap.Foreground; for dots = 1 : length(crDriftPre) @@ -1018,8 +1046,8 @@ else if isfield(MRSCont.QM.drift.pre, which_sub_spec) && ~strcmp(which_spec,'mm') if length(MRSCont.QM.drift.pre.(which_sub_spec){kk}) > 1 - crDriftPre = MRSCont.QM.drift.pre.(which_sub_spec){kk} + MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/applyDataToPlot.txfrq*1e6; - crDriftPost = MRSCont.QM.drift.post.(which_sub_spec){kk} + MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/applyDataToPlot.txfrq*1e6; + crDriftPre = MRSCont.QM.drift.pre.(which_sub_spec){kk}; + crDriftPost = MRSCont.QM.drift.post.(which_sub_spec){kk}; hold(ax_drift, 'on'); colors = ones(length(crDriftPre),1).*colormap.Foreground; for dots = 1 : length(crDriftPre) @@ -1027,11 +1055,11 @@ colors(dots,2) = colors(dots,2) + (1 - colors(dots,2)) * (1-weights(dots,1)); colors(dots,3) = colors(dots,3) + (1 - colors(dots,3)) * (1-weights(dots,1)); end - scatter(ax_drift, [1:length(crDriftPre)],crDriftPre'-(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6),36,ones(length(crDriftPre),1).*colormap.LightAccent); - scatter(ax_drift, [1:length(crDriftPost)],crDriftPost'-(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6),36,colors,'filled','MarkerEdgeColor',colormap.Foreground); + scatter(ax_drift, [1:length(crDriftPre)],crDriftPre+(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6),36,ones(length(crDriftPre),1).*colormap.LightAccent); + scatter(ax_drift, [1:length(crDriftPost)],crDriftPost+(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6),36,colors,'filled','MarkerEdgeColor',colormap.Foreground); - text(ax_drift, length(crDriftPre)*1.05, crDriftPre(end)-(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6), 'Pre', 'Color', colormap.LightAccent); - text(ax_drift, length(crDriftPost)*1.05, crDriftPost(end)-(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6), 'Post', 'Color', colormap.Foreground); + text(ax_drift, length(crDriftPre)*1.05, crDriftPre(end)+(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6), 'Pre', 'Color', colormap.LightAccent); + text(ax_drift, length(crDriftPost)*1.05, crDriftPost(end)+(MRSCont.QM.freqShift.(which_spec)(1,kk,SubSpectraIndex)/procDataToPlot.txfrq*1e6), 'Post', 'Color', colormap.Foreground); set(ax_drift, 'YLim', [3.028-0.1 3.028+0.1]); yticks([3.028-0.08 3.028-0.04 3.028 3.028+0.04 3.028+0.08]); yticklabels({'2.94' '2.98' '3.02' '3.06' '3.10'}); diff --git a/process/OspreyProcess.m b/process/OspreyProcess.m index 3556df94..d29d27ac 100644 --- a/process/OspreyProcess.m +++ b/process/OspreyProcess.m @@ -458,6 +458,7 @@ % raw struct raw=op_HadamardScans(raw,[-1 1],'diff1'); raw=op_HadamardScans(raw,[1 1],'sum'); + raw_no_subspec_aling = raw; if ~raw.flags.isMRSI % [raw,~] = op_phaseCrCho(raw, 1); [refShift_SubSpecAlign, ~] = osp_XReferencing(raw,[3.03 3.22],[1 1],[1.85 4.2]);% determine frequency shift @@ -480,9 +481,6 @@ case 'L2Norm' [raw] = osp_editSubSpecAlign(raw, seq, target,MRSCont.opts.UnstableWater); otherwise - if ~MRSCont.flags.isMRSI - raw_B = op_addphase(raw_B, -ph*180/pi, 0, raw_B.centerFreq, 1); - end end if MRSCont.flags.hasMM @@ -512,6 +510,7 @@ % Create the sum spectrum raw=op_HadamardScans(raw,[1 1],'sum'); raw.target = target; + raw = op_freqshift(raw, refShift_SubSpecAlign); end if raw.flags.isHERMES || raw.flags.isHERCULES @@ -526,6 +525,7 @@ raw=op_HadamardScans(raw,[-1 1 -1 1],'diff3'); end raw=op_HadamardScans(raw,[1 1 1 1],'sum'); + raw_no_subspec_aling = raw; % Correct the frequency axis so that Cr appears at 3.027 ppm [refShift_SubSpecAlign, ~] = osp_XReferencing(raw,[3.03 3.22],[1 1],[1.85 4.2]);% determine frequency shift @@ -561,7 +561,7 @@ raw=op_HadamardScans(raw,[-1 1 -1 1],'diff3'); end raw=op_HadamardScans(raw,[1 1 1 1],'sum'); - + raw = op_freqshift(raw, refShift_SubSpecAlign); end %% %%% 7. REMOVE RESIDUAL WATER %%% @@ -576,6 +576,7 @@ end % Apply iterative water filter raw = op_iterativeWaterFilter(raw, waterRemovalFreqRange, 32, fracFID*length(raw.fids), 0); + raw_no_subspec_aling = op_iterativeWaterFilter(raw_no_subspec_aling, waterRemovalFreqRange, 32, fracFID*length(raw.fids), 0); if MRSCont.flags.hasMM %re_mm raw_mm = op_iterativeWaterFilter(raw_mm, waterRemovalFreqRange, 32, fracFID*length(raw_mm.fids), 0); @@ -588,6 +589,7 @@ [refShift, ~] = osp_XReferencing(raw,2.01,1,[1.85 4.2]);% determine frequency shift end [raw] = op_freqshift(raw,-refShift); % Reference spectra by cross-correlation + [raw_no_subspec_aling] = op_freqshift(raw_no_subspec_aling,-refShift); % Reference spectra by cross-correlation if MRSCont.flags.hasMM %re_mm if raw_mm.flags.isMEGA @@ -613,6 +615,7 @@ [raw,SNR] = op_get_Multispectra_SNR(raw); FWHM = op_get_Multispectra_LW(raw); MRSCont.processed.metab{metab_ll,kk} = raw; + MRSCont.processed_no_align.metab{metab_ll,kk} = raw_no_subspec_aling; for ss = 1 : length(SubSpec) MRSCont.QM.SNR.metab(metab_ll,kk,ss) = SNR{ss}; MRSCont.QM.FWHM.metab(metab_ll,kk,ss) = FWHM(ss); % in Hz @@ -620,11 +623,11 @@ MRSCont.QM.res_water_amp.metab(metab_ll,kk,ss) = sum(MRSCont.processed.metab{kk}.watersupp{ss}.amp); if strcmp(SubSpec{ss},'diff1') ||strcmp(SubSpec{ss},'diff2') || strcmp(SubSpec{ss},'diff3') ||strcmp(SubSpec{ss},'sum') if raw.flags.isMEGA - MRSCont.QM.drift.pre.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.pre.A{kk}'; MRSCont.QM.drift.pre.B{kk}'], [], 1)'; - MRSCont.QM.drift.post.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.post.A{kk}'; MRSCont.QM.drift.post.B{kk}'], [], 1)'; + MRSCont.QM.drift.pre.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.pre.A{kk}'; MRSCont.QM.drift.pre.B{kk}'], 1, [])'; + MRSCont.QM.drift.post.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.post.A{kk}'; MRSCont.QM.drift.post.B{kk}'], 1, [])'; else - MRSCont.QM.drift.pre.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.pre.A{kk}'; MRSCont.QM.drift.pre.B{kk}'; MRSCont.QM.drift.pre.C{kk}'; MRSCont.QM.drift.pre.D{kk}'], [], 1)'; - MRSCont.QM.drift.post.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.post.A{kk}'; MRSCont.QM.drift.post.B{kk}'; MRSCont.QM.drift.pre.C{kk}'; MRSCont.QM.drift.pre.D{kk}'], [], 1)'; + MRSCont.QM.drift.pre.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.pre.A{kk}'; MRSCont.QM.drift.pre.B{kk}'; MRSCont.QM.drift.pre.C{kk}'; MRSCont.QM.drift.pre.D{kk}'], 1, [])'; + MRSCont.QM.drift.post.(SubSpec{ss}){metab_ll,kk} = reshape([MRSCont.QM.drift.post.A{kk}'; MRSCont.QM.drift.post.B{kk}'; MRSCont.QM.drift.post.C{kk}'; MRSCont.QM.drift.post.D{kk}'], 1, [])'; end end MRSCont.QM.drift.pre.AvgDeltaCr.(SubSpec{ss})(metab_ll,kk) = mean(MRSCont.QM.drift.pre.(SubSpec{ss}){kk} - 3.02);