Skip to content
Snippets Groups Projects
Commit d41c227f authored by Tamas Borbath's avatar Tamas Borbath
Browse files

Adapted plotBaselines_and_fitted_spectra.m to include also the difference... 

Adapted plotBaselines_and_fitted_spectra.m to include also the difference spectrum between baselines
parent 0e65526e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
createTestData/plotBaselines_and_fitted_spectra.m
+ 142
− 138
View file @ d41c227f
......@@ -5,7 +5,7 @@ exportFolder = {'11 Simulations'};
pathNameExport = 'ProFit test data';
pathBaseExportFiles = pathToDataFolder(pathNameExport, exportFolder);
pathBaseExportFiles = pathBaseExportFiles(1:end-1);
pathBaseExportFiles = [pathBaseExportFiles ' - R3\'];
pathBaseExportFiles = [pathBaseExportFiles ' - R3_new\'];
dataExportPathBase = strcat(pathBaseExportFiles, exportFolder{1}, '\');
plotPath = [dataExportPathBase, 'Plots\BaselinePlots\'];
......@@ -48,7 +48,8 @@ backgroundInlayColor = [0.8 0.8 0.8];
LineWidth = 1.5;
offsetBaseline = 7; % percent of max peak
offsetResidual = 20; % percent of max peak
offsetBaselineDiff = 18; % percent of max peak
offsetResidual = 30; % percent of max peak
offsetPlots = 4; % percent of max peak
subplotGridX = 2;
subplotGridY = 2;
......@@ -63,114 +64,125 @@ for indexParam = 1:length(baselines)
% figure;
figure('units','normalized','outerposition',[0.13,-0.014,0.85,0.99]); %Works on the Dell screen
end
indexParamSorted = sortedFigure(indexParam);
%% load data: Input data, LCModel fit, ProFit Fit
% import the spectra but also the used concentrations and all parameters
fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) '_truncOff'];
load([dataExportPath, fileName, '.mat'], 'summedFidWithNoise', 'current_em', 'currentConcentrationRm', ...
'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
if plot_dkntmn
fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) truncSuffix];
[~, ~, phasedData_lcm_05, fitData_05, baselineData_lcm_05, residualData_05, ...
~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
extractFits(dataExportPath_05,[fileName, '.coord']);
[~, ~, phasedData_lcm_015, fitData_015, baselineData_lcm_015, residualData_015, ...
~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
extractFits(dataExportPath_015,[fileName, '.coord']);
end
[nOfPoints, ppmVector_lcm, phasedData_lcm, fitData, baselineData_lcm, residualData, ...
metaboliteNames, metaboliteData, tableConcentrations, ppmStart_lcm, ppmEnd_lcm, ...
scanFrequency, frequencyShift, ppmGap, phase0_deg_lcm, phase1_deg_ppm_lcm] = ...
extractFits(dataExportPath,[fileName, '.coord']);
load(strcat(dataExportPath,fileName, '_profit.mat'),'fitresult', 'data');
%% process input data
nTimepoints = length(summedFidWithNoise);
ppmVector = ppmVectorFunction(scanFrequency*1e6, bandwidth, nTimepoints, '1H');
indexParamSorted = sortedFigure(indexParam);
%% load data: Input data, LCModel fit, ProFit Fit
% import the spectra but also the used concentrations and all parameters
fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) '_truncOff'];
load([dataExportPath, fileName, '.mat'], 'summedFidWithNoise', 'current_em', 'currentConcentrationRm', ...
'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
if plot_dkntmn
fileName = ['Simulated_', paramKeyword, '_', num2str(indexParamSorted) truncSuffix];
[~, ~, phasedData_lcm_05, fitData_05, baselineData_lcm_05, residualData_05, ...
~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
extractFits(dataExportPath_05,[fileName, '.coord']);
[~, ~, phasedData_lcm_015, fitData_015, baselineData_lcm_015, residualData_015, ...
~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = ...
extractFits(dataExportPath_015,[fileName, '.coord']);
end
[nOfPoints, ppmVector_lcm, phasedData_lcm, fitData, baselineData_lcm, residualData, ...
metaboliteNames, metaboliteData, tableConcentrations, ppmStart_lcm, ppmEnd_lcm, ...
scanFrequency, frequencyShift, ppmGap, phase0_deg_lcm, phase1_deg_ppm_lcm] = ...
extractFits(dataExportPath,[fileName, '.coord']);
spectrum = fftshift(fft(summedFidWithNoise));
baselineSpectrum = fftshift(fft(baselineScaled));
indexEnd = find(ppmVector > ppmEnd_lcm, 1);
spectrumInputOffset = real(spectrum(indexEnd));
baselineInputOffset = real(baselineSpectrum(indexEnd));
%% process ProFit data
currentFitresult = fitresult{1,iteration};
ppm_profit = currentFitresult.common_data.data.fitroi.ppm2;
%phasedData
pData_profit = currentFitresult.processed_spec;
%residual
rData_profit = currentFitresult.residual_spec;
%spline baseline
if ~isempty(currentFitresult.splines_baseline)
spBB = currentFitresult.splines_baseline.basis_matrix;
concSpline = currentFitresult.splines_baseline.coeff;
bData_profit = (spBB * concSpline(1:size(spBB,2)))';
splines_lambda = currentFitresult.splines_lambda;
splines_ed_vec = currentFitresult.splines_ed_vec;
splines_optim_model_vec = currentFitresult.splines_optim_model_vec;
splines_optim_idx = currentFitresult.splines_optim_idx;
else
bData_profit = zeros(size(currentFitresult.residual_spec));
splines_lambda = 0;
splines_ed_vec = [];
splines_optim_model_vec = [];
splines_optim_idx = 0;
end
load(strcat(dataExportPath,fileName, '_profit.mat'),'fitresult', 'data');
%% process input data
nTimepoints = length(summedFidWithNoise);
ppmVector = ppmVectorFunction(scanFrequency*1e6, bandwidth, nTimepoints, '1H');
spectrum = fftshift(fft(summedFidWithNoise));
baselineSpectrum = fftshift(fft(baselineScaled));
indexEnd = find(ppmVector > ppmEnd_lcm, 1);
spectrumInputOffset = real(spectrum(indexEnd));
baselineInputOffset = real(baselineSpectrum(indexEnd));
%% process ProFit data
currentFitresult = fitresult{1,iteration};
ppm_profit = currentFitresult.common_data.data.fitroi.ppm2;
%phasedData
pData_profit = currentFitresult.processed_spec;
%residual
rData_profit = currentFitresult.residual_spec;
%spline baseline
if ~isempty(currentFitresult.splines_baseline)
spBB = currentFitresult.splines_baseline.basis_matrix;
concSpline = currentFitresult.splines_baseline.coeff;
bData_profit = (spBB * concSpline(1:size(spBB,2)))';
splines_lambda = currentFitresult.splines_lambda;
splines_ed_vec = currentFitresult.splines_ed_vec;
splines_optim_model_vec = currentFitresult.splines_optim_model_vec;
splines_optim_idx = currentFitresult.splines_optim_idx;
else
bData_profit = zeros(size(currentFitresult.residual_spec));
splines_lambda = 0;
splines_ed_vec = [];
splines_optim_model_vec = [];
splines_optim_idx = 0;
end
%% Scale for plotting
% calculate scaling. Maximum has to be from the same fitted range as in LCModel/ProFit
ppmMask = (ppmVector < ppmStart_lcm) & (ppmVector > ppmEnd_lcm);
maxSpectrum = max(real(spectrum(ppmMask)));
maxSpectrum_lcm = max(phasedData_lcm);
maxSpectrum_profit = max(real(pData_profit));
scalingFactor_lcm = maxSpectrum ./ maxSpectrum_lcm;
scalingFactor_profit = maxSpectrum ./ maxSpectrum_profit;
minBaseline = min(real(baselineSpectrum(ppmMask)));
%scale the offsets
offsetBaselineScaled = offsetBaseline * maxSpectrum / 100;
offsetBaselineDiffScaled = offsetBaselineDiff * maxSpectrum / 100 - minBaseline;
offsetResidualScaled = offsetResidual * maxSpectrum / 100 - minBaseline;
offsetPlotsScaled = offsetPlots * maxSpectrum / 100;
%% create the figure
indexSubplot = mod(indexParam-1, nGrid)+1;
hs = subplot(subplotGridX,subplotGridY,indexSubplot);
coordinates = hs.Position;
hold on; %if not alligned, mind that the ppmVector has to be defined with the water resonance at 4.7 ppm
plot(ppmVector,real(spectrum),'LineWidth', LineWidth, 'Color', colorSpectrum);
% plot(ppmVector_lcm,real(phasedData_lcm)*scalingFactor_lcm,'LineWidth', LineWidth, 'Color', colorLCModel);
% plot(ppm_profit,real(pData_profit)*scalingFactor_profit,'LineWidth', LineWidth, 'Color', colorProFit);
if plot_dkntmn
offsetBaselineScaled = offsetBaselineScaled + offsetPlotsScaled;
plot(ppmVector_lcm,real(baselineData_lcm_05)*scalingFactor_lcm-offsetBaselineScaled+2.4*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_05);
plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+1.6*offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
plot(ppmVector_lcm,real(baselineData_lcm_015)*scalingFactor_lcm-offsetBaselineDiffScaled*2+0.8*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel);
else
plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
end
plot(ppmVector,real(baselineSpectrum)-offsetBaselineScaled,'-.','LineWidth', LineWidth, 'Color', colorOriginal);
plot(ppm_profit,real(bData_profit)*scalingFactor_profit-offsetBaselineScaled-offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorProFit);
fitBaselineLCM = fit(ppmVector_lcm, real(baselineData_lcm)*scalingFactor_lcm, 'linearinterp');
fitBaselineSpectrum = fit(ppmVector', real(baselineSpectrum)', 'linearinterp');
fitBaselineProFit = fit(ppm_profit', real(bData_profit)'*scalingFactor_profit, 'linearinterp');
plot(ppmVector_lcm, fitBaselineLCM(ppmVector_lcm)-fitBaselineSpectrum(ppmVector_lcm)-offsetBaselineDiffScaled,':','LineWidth', LineWidth, 'Color', colorLCModel)
plot(ppm_profit,fitBaselineProFit(ppm_profit)-fitBaselineSpectrum(ppm_profit)-offsetBaselineDiffScaled-offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorProFit);
plot(ppmVector_lcm,real(residualData)*scalingFactor_lcm-offsetResidualScaled,'LineWidth', LineWidth, 'Color', colorLCModel);
plot(ppm_profit,real(rData_profit)*scalingFactor_profit-offsetResidualScaled-offsetPlotsScaled,'LineWidth', LineWidth, 'Color', colorProFit);
set(gca,'xdir','reverse');
xlim([ppmEnd_lcm, ppmStart_lcm]);
ylim([-offsetResidualScaled-3*offsetPlotsScaled maxSpectrum+5*offsetPlotsScaled])
xlabel('\delta (ppm)');
ylabel('Signal (arb. u.)');
set(gca,'ytick',[]);
%% Scale for plotting
% calculate scaling. Maximum has to be from the same fitted range as in LCModel/ProFit
ppmMask = (ppmVector < ppmStart_lcm) & (ppmVector > ppmEnd_lcm);
maxSpectrum = max(real(spectrum(ppmMask)));
maxSpectrum_lcm = max(phasedData_lcm);
maxSpectrum_profit = max(real(pData_profit));
scalingFactor_lcm = maxSpectrum ./ maxSpectrum_lcm;
scalingFactor_profit = maxSpectrum ./ maxSpectrum_profit;
minBaseline = min(real(baselineSpectrum(ppmMask)));
%scale the offsets
offsetBaselineScaled = offsetBaseline * maxSpectrum / 100;
offsetResidualScaled = offsetResidual * maxSpectrum / 100 - minBaseline;
offsetPlotsScaled = offsetPlots * maxSpectrum / 100;
%% create the figure
indexSubplot = mod(indexParam-1, nGrid)+1;
hs = subplot(subplotGridX,subplotGridY,indexSubplot);
coordinates = hs.Position;
hold on; %if not alligned, mind that the ppmVector has to be defined with the water resonance at 4.7 ppm
plot(ppmVector,real(spectrum),'LineWidth', LineWidth, 'Color', colorSpectrum);
% plot(ppmVector_lcm,real(phasedData_lcm)*scalingFactor_lcm,'LineWidth', LineWidth, 'Color', colorLCModel);
% plot(ppm_profit,real(pData_profit)*scalingFactor_profit,'LineWidth', LineWidth, 'Color', colorProFit);
if plot_dkntmn
offsetBaselineScaled = offsetBaselineScaled + offsetPlotsScaled;
plot(ppmVector_lcm,real(baselineData_lcm_05)*scalingFactor_lcm-offsetBaselineScaled+2.4*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_05);
plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+1.6*offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
plot(ppmVector_lcm,real(baselineData_lcm_015)*scalingFactor_lcm-offsetBaselineScaled+0.8*offsetPlotsScaled,':','LineWidth', LineWidth, 'Color', colorLCModel_015);
else
plot(ppmVector_lcm,real(baselineData_lcm)*scalingFactor_lcm-offsetBaselineScaled+offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorLCModel);
end
plot(ppmVector,real(baselineSpectrum)-offsetBaselineScaled,'-.','LineWidth', LineWidth, 'Color', colorOriginal);
plot(ppm_profit,real(bData_profit)*scalingFactor_profit-offsetBaselineScaled-offsetPlotsScaled,'-.','LineWidth', LineWidth, 'Color', colorProFit);
plot(ppmVector_lcm,real(residualData)*scalingFactor_lcm-offsetResidualScaled,'LineWidth', LineWidth, 'Color', colorLCModel);
plot(ppm_profit,real(rData_profit)*scalingFactor_profit-offsetResidualScaled-offsetPlotsScaled,'LineWidth', LineWidth, 'Color', colorProFit);
set(gca,'xdir','reverse');
xlim([ppmEnd_lcm, ppmStart_lcm]);
ylim([-offsetResidualScaled-3*offsetPlotsScaled maxSpectrum+5*offsetPlotsScaled])
xlabel('\delta (ppm)');
ylabel('Signal (arb. u.)');
set(gca,'ytick',[]);
if plot_dkntmn
legend('Spectrum',sprintf('LCModel baseline\ndkntmn=0.5'),sprintf('LCModel baseline\ndkntmn=0.25'),sprintf('LCModel baseline\ndkntmn=0.15'),'Input baseline','ProFit baseline','LCModel residual','ProFit residual')
else
legend('Spectrum','LCModel baseline','Input baseline','ProFit baseline','LCModel residual','ProFit residual')
% legend('Spectrum','LCModel baseline','Input baseline','ProFit baseline','LCModel residual','ProFit residual')
if mod(indexParam,nGrid)==2
handle_legend = legend('Spectrum','LCModel baseline','Input baseline','ProFit-1D baseline',['LCModel ' char(8211) ' Input baseline'],['ProFit-1D ' char(8211) ' Input baseline'], 'LCModel residual','ProFit-1D residual','FontSize',12);
title(handle_legend,'Legend');
end
end
%arrows residual
......@@ -188,7 +200,13 @@ for indexParam = 1:length(baselines)
%arrows Input baseline - ProFit
yUpper = baselineInputOffset-offsetBaselineScaled+offsetPlotsScaled;
yLower = baselineInputOffset-offsetBaselineScaled;
xa = [ppmEnd_lcm ppmEnd_lcm]-0.1;
xa = [ppmEnd_lcm ppmEnd_lcm]-0.12;
doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
%arrows Baseline differences ProFit & LCModel
yUpper = baselineInputOffset-offsetBaselineDiffScaled;
yLower = baselineInputOffset-offsetBaselineDiffScaled-offsetPlotsScaled;
xa = [ppmEnd_lcm ppmEnd_lcm]-0.12;
doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
% text box with label "Plot offsets"
......@@ -197,12 +215,13 @@ for indexParam = 1:length(baselines)
[xaf,yaf] = axescoord2figurecoord(xa,ya);
str = {'Plot','offsets'};
annotation('textbox',[xaf(1) yaf(1)-0.02 0.1 0.1],'String',str,'FitBoxToText','On', ...
'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor);
'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor,'FontSize',13);
set(gcf,'Color','w')
set(gca,'FontSize',13)
annotation('rectangle',[coordinates(1)+coordinates(3)*0.02 coordinates(2)+coordinates(4)*0.62 ...
.35./subplotGridX .4./subplotGridY],'FaceColor', 'w', 'FaceAlpha',0, 'Color',[0 0.4 0.1])
% annotation('rectangle',[coordinates(1)+coordinates(3)*0.02 coordinates(2)+coordinates(4)*0.62 ...
% .35./subplotGridX .4./subplotGridY],'FaceColor', 'w', 'FaceAlpha',0, 'Color',[0 0.4 0.1])
axes('Position',[coordinates(1)+coordinates(3)*0.15 coordinates(2)+coordinates(4)*0.81 ...
.2./subplotGridX .2./subplotGridY])
......@@ -210,20 +229,28 @@ for indexParam = 1:length(baselines)
plot(splines_ed_vec, splines_optim_model_vec,'Color',[0. 0.5 0.2]);
hold on
if ~isempty(splines_ed_vec(splines_optim_idx))
xl = xline(splines_ed_vec(splines_optim_idx),'-.', {['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx))],...
['mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]}, ...
% xl = xline(splines_ed_vec(splines_optim_idx),'-.', {['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx))],...
% ['mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]}, ...
% 'HandleVisibility','off', 'LineWidth', 2);
% xl.LabelVerticalAlignment = 'top';
% xl.LabelHorizontalAlignment = 'center';
% xl.FontSize = 9;
xl = xline(splines_ed_vec(splines_optim_idx),'-.',...
'HandleVisibility','off', 'LineWidth', 2);
xl.LabelVerticalAlignment = 'top';
xl.LabelHorizontalAlignment = 'center';
xl.FontSize = 9;
end
set(gca,'FontSize',10)
set(gca,'FontSize',12)
set(gca,'Xscale','log');
% set(gca,'Color',backgroundInlayColor)
xlabel('Baseline ED / ppm');
ylabel('mAIC');
ylim([min(splines_optim_model_vec)*0.99 max(splines_optim_model_vec)*1.01])
title('ProFit Baseline Smoothness');
% title('ProFit Baseline Smoothness');
if isempty(splines_ed_vec(splines_optim_idx))
title('ProFit Baseline Smoothness');
else
title(['ED: ' sprintf('%.f',splines_ed_vec(splines_optim_idx)),...
' mAIC: ', sprintf('%.1f',splines_optim_model_vec(splines_optim_idx))]);
end
annotation('textbox',[coordinates(1)+coordinates(3)*0.02-0.03 coordinates(2)+coordinates(4)*0.62 ...
.35./subplotGridX .4./subplotGridY],'String',subplotNames{indexSubplot},'FontSize', 16, 'FontWeight', 'bold','LineStyle','none')
......@@ -233,26 +260,3 @@ for indexParam = 1:length(baselines)
end
end
end
function doublePointedArrowInside(xa,yUpper,yLower, arrowsColor)
set(gcf,'Units','normalized');
yaDown = [yUpper+1e-8 yUpper];
yaLine = [yLower yUpper];
yaUp = [yLower-1e-8 yLower];
[xaf,yafLine] = axescoord2figurecoord(xa,yaLine);
line = annotation('line',xaf,yafLine);
[xaf,yafDown] = axescoord2figurecoord(xa,yaDown);
arDown = annotation('arrow',xaf,yafDown);
[xaf,yafUp] = axescoord2figurecoord(xa,yaUp);
arUp = annotation('arrow',xaf,yafUp);
%styling
line.Color = arrowsColor;
arDown.HeadStyle = 'vback2';
arDown.HeadLength = 5;
arDown.HeadWidth = 7;
arDown.Color = arrowsColor;
arUp.HeadStyle = 'vback2';
arUp.HeadLength = 5;
arUp.HeadWidth = 7;
arUp.Color = arrowsColor;
end
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment