Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[FEATURE REQUEST] - Add pre/post alignment - Eric Porges #417

Merged
merged 1 commit into from
May 23, 2022
Merged

[FEATURE REQUEST] - Add pre/post alignment - Eric Porges #417

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion GUI/osp_updateProWindow.m
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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