Skip to content

Commit

Permalink
Merge pull request #417 from HJZollner/develop
Browse files Browse the repository at this point in the history 
[FEATURE REQUEST] - Add pre/post alignment - Eric Porges
  • Loading branch information
@HJZollner
HJZollner authored May 23, 2022
2 parents f96f24d + f64a568 commit ddc906d
Show file tree
Hide file tree
Showing 3 changed files with 50 additions and 19 deletions.
2 changes: 1 addition & 1 deletion GUI/osp_updateProWindow.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
48 changes: 38 additions & 10 deletions plot/osp_plotProcess.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)));
Expand Down Expand Up @@ -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<ppmmax))), 'Color',MRSCont.colormap.LightAccent, 'LineWidth', 1);
end
text(ax_proc, ppmmin+1.5, 1, 'pre alignment', 'Color', colormap.LightAccent);
text(ax_proc, ppmmin+1.5, 0.9, 'post alignment', 'Color', colormap.Foreground);
end
if isfield(MRSCont,'plot') && isfield(MRSCont.plot, 'processed') && (MRSCont.plot.processed.match == 1)
if MRSCont.flags.hasRef || MRSCont.flags.hasWater
if MRSCont.flags.hasRef
Expand Down Expand Up @@ -948,8 +976,8 @@
if (isfield(MRSCont.flags,'isPRIAM') && (MRSCont.flags.isPRIAM == 1))
if isfield(MRSCont.QM{VoxelIndex}.drift.pre, which_sub_spec)
if length(MRSCont.QM{VoxelIndex}.drift.pre.(which_sub_spec){kk}) > 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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -1018,20 +1046,20 @@
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)
colors(dots,1) = colors(dots,1) + (1 - colors(dots,1)) * (1-weights(dots,1));
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'});
Expand Down
19 changes: 11 additions & 8 deletions process/OspreyProcess.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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 %%%
Expand All @@ -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);
Expand All @@ -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
Expand All @@ -613,18 +615,19 @@
[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
MRSCont.QM.freqShift.metab(metab_ll,kk,ss) = refShift;
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);
Expand Down

0 comments on commit ddc906d

Please sign in to comment.