diff --git a/createTestData/plotBaselines_and_fitted_spectra.m b/createTestData/plotBaselines_and_fitted_spectra.m index 3d519c0652607eca535698b46254f40c1b20a026..100c6a07b480efbc714824961d9dca73feedd113 100644 --- a/createTestData/plotBaselines_and_fitted_spectra.m +++ b/createTestData/plotBaselines_and_fitted_spectra.m @@ -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