Modified g2 feature extraction to extract contrast and peak angles as intended with plotting modifications.
This commit is contained in:
parent
0401dd3f7b
commit
f41663c801
@ -2,8 +2,9 @@ function results = analyzeAutocorrelation(autocorrresults)
|
||||
%% Computes features from g2(theta) curves and aggregates results.
|
||||
% Assumes theta_values are given in radians.
|
||||
|
||||
feature_list = table();
|
||||
per_group = cell(numel(autocorrresults.scan_parameter_values),1);
|
||||
feature_list = table();
|
||||
per_group = cell(numel(autocorrresults.scan_parameter_values),1);
|
||||
per_group_mean = cell(numel(autocorrresults.scan_parameter_values),1); % mean-curve features
|
||||
|
||||
% --- Loop over groups ---
|
||||
for i = 1:numel(autocorrresults.scan_parameter_values)
|
||||
@ -18,6 +19,10 @@ function results = analyzeAutocorrelation(autocorrresults)
|
||||
feature_list = [feature_list; feats]; %#ok<AGROW>
|
||||
end
|
||||
per_group{i} = group_tbl;
|
||||
|
||||
% --- Compute mean-curve features for this group ---
|
||||
meanFeatures = extractMeanG2Features(thetas, curves);
|
||||
per_group_mean{i} = meanFeatures;
|
||||
end
|
||||
|
||||
%% Compute group summaries (mean & SEM of features)
|
||||
@ -28,22 +33,8 @@ function results = analyzeAutocorrelation(autocorrresults)
|
||||
FT = per_group{i};
|
||||
|
||||
% Scalars
|
||||
vars = {'A2','A4','A6','S2','Q4','H6','Contrast','NumberOfPeaks','MaxProm','MaxHeight'};
|
||||
vars = {'A2','A4','A6','S2','Q4','H6','Contrast'};
|
||||
|
||||
% --- Handle spacing ---
|
||||
spacing_means = cellfun(@mean, FT.Spacing, 'UniformOutput', false);
|
||||
spacing_stds = cellfun(@std, FT.Spacing, 'UniformOutput', false);
|
||||
FT.SpacingMean = cell2mat(spacing_means);
|
||||
FT.SpacingStd = cell2mat(spacing_stds);
|
||||
|
||||
% --- Handle widths ---
|
||||
width_means = cellfun(@mean, FT.Widths, 'UniformOutput', false);
|
||||
width_stds = cellfun(@std, FT.Widths, 'UniformOutput', false);
|
||||
FT.WidthsMean = cell2mat(width_means);
|
||||
FT.WidthsStd = cell2mat(width_stds);
|
||||
|
||||
vars = [vars, {'SpacingMean','SpacingStd','WidthsMean','WidthsStd'}];
|
||||
|
||||
% Group averages
|
||||
meanVals = varfun(@mean, FT, 'InputVariables', vars);
|
||||
semVals = varfun(@(x) std(x,0,1,'omitnan')/sqrt(size(FT,1)), ...
|
||||
@ -70,11 +61,12 @@ function results = analyzeAutocorrelation(autocorrresults)
|
||||
|
||||
%% Package results
|
||||
results = struct();
|
||||
results.features_all = feature_list; % all repetitions pooled
|
||||
results.features_group = per_group; % per scan parameter
|
||||
results.summaries = summary_table; % group-level stats
|
||||
results.g2_cumulants = g2_cums; % cumulants from full distribution
|
||||
results.theta_values = autocorrresults.theta_values;
|
||||
results.features_all = feature_list; % all repetitions pooled
|
||||
results.features_group = per_group; % per scan parameter
|
||||
results.meanCurve = per_group_mean; % mean-curve contrast & peaks
|
||||
results.summaries = summary_table; % group-level stats
|
||||
results.g2_cumulants = g2_cums; % cumulants from full distribution
|
||||
results.theta_values = autocorrresults.theta_values;
|
||||
end
|
||||
|
||||
function features = extractFeatures(thetas, g2curve)
|
||||
@ -104,36 +96,74 @@ function features = extractFeatures(thetas, g2curve)
|
||||
S2 = 0; Q4 = 0; H6 = 0;
|
||||
end
|
||||
|
||||
%% --- Step 4: RMS fluctuation (contrast) ---
|
||||
RMSval = sqrt(mean(g2cWithinPi.^2));
|
||||
|
||||
%% --- Step 5: Peak finding with dynamic threshold ---
|
||||
alpha = 0.25;
|
||||
dynamicProm = alpha * range(g2WithinPi);
|
||||
minDist = round(numel(thetasWithinPi)/12);
|
||||
|
||||
[pks,locs,w,prom] = findpeaks(g2WithinPi, ...
|
||||
'MinPeakProminence', dynamicProm, ...
|
||||
'MinPeakDistance', minDist);
|
||||
|
||||
peakAngles = thetasWithinPi(locs);
|
||||
|
||||
%% --- Step 6: Spacing calculation (wrap within π) ---
|
||||
if numel(peakAngles) > 1
|
||||
d = diff([peakAngles, peakAngles(1)+pi]);
|
||||
%% --- Step 4: Contrast ---
|
||||
g2_restricted = g2curve((thetas >= pi/18 & thetas <= 17*pi/18));
|
||||
|
||||
if ~isempty(g2_restricted)
|
||||
ContrastVal = (max(g2_restricted) - min(g2_restricted)) / ...
|
||||
(max(g2_restricted) + min(g2_restricted));
|
||||
else
|
||||
d = NaN;
|
||||
ContrastVal = NaN;
|
||||
end
|
||||
|
||||
%% --- Step 7: Safely handle empty peaks ---
|
||||
if isempty(prom), MaxProm = NaN; else, MaxProm = max(prom(:),[],'omitnan'); end
|
||||
if isempty(pks), MaxHeight = NaN; else, MaxHeight = max(pks(:),[],'omitnan'); end
|
||||
|
||||
%% --- Step 8: Package as one row table ---
|
||||
features = table( ...
|
||||
An(1), An(2), An(3), RMSval, numel(pks), ...
|
||||
MaxProm, MaxHeight, ...
|
||||
{d(:)'}, {w(:)'}, ...
|
||||
An(1), An(2), An(3), ContrastVal, ...
|
||||
S2, Q4, H6, ...
|
||||
'VariableNames', {'A2','A4','A6','Contrast','NumberOfPeaks','MaxProm','MaxHeight','Spacing','Widths','S2','Q4','H6'});
|
||||
end
|
||||
'VariableNames', {'A2','A4','A6','Contrast','S2','Q4','H6'});
|
||||
end
|
||||
|
||||
function meanFeatures = extractMeanG2Features(thetas, g2_curves)
|
||||
% Computes contrast and peak angles from the mean g2 curve.
|
||||
% g2_curves: [N_reps × N_theta] array of individual g2(theta) curves
|
||||
% thetas: vector of theta values in radians
|
||||
|
||||
% --- Step 1: Compute mean curve ---
|
||||
meanG2 = mean(g2_curves, 1);
|
||||
|
||||
% --- Step 2: Contrast calculation ---
|
||||
maskContrast = (thetas >= pi/18 & thetas <= 17*pi/18);
|
||||
meanG2_contrast = meanG2(maskContrast);
|
||||
|
||||
if ~isempty(meanG2_contrast)
|
||||
ContrastVal = (max(meanG2_contrast) - min(meanG2_contrast)) / ...
|
||||
(max(meanG2_contrast) + min(meanG2_contrast));
|
||||
else
|
||||
ContrastVal = NaN;
|
||||
end
|
||||
|
||||
% --- Step 3: Peak finding restricted to (pi/18, pi/2) ---
|
||||
maskPeaks = (thetas >= pi/18 & thetas <= pi/2);
|
||||
meanG2_peaks = meanG2(maskPeaks);
|
||||
thetas_peaks = thetas(maskPeaks);
|
||||
|
||||
% Dynamic prominence: 10% of range or minimum 0.05 absolute
|
||||
alpha = 0.1;
|
||||
minAbsProm = 0.005;
|
||||
dynamicProm = max(alpha * range(meanG2_peaks), minAbsProm);
|
||||
minDist = max(round(numel(thetas_peaks)/12),1); % at least 1
|
||||
|
||||
[~, locs] = findpeaks(meanG2_peaks, ...
|
||||
'MinPeakProminence', dynamicProm, ...
|
||||
'MinPeakDistance', minDist);
|
||||
|
||||
PeakAngles = thetas_peaks(locs); % angles of detected peaks
|
||||
|
||||
% --- Step 4: Fallback if no peaks detected ---
|
||||
if isempty(PeakAngles)
|
||||
% Use global maximum of restricted curve
|
||||
[~, idxMax] = max(meanG2_peaks);
|
||||
MaxPeakAngle = thetas_peaks(idxMax);
|
||||
else
|
||||
% Otherwise, take max of detected peaks
|
||||
[~, idxMax] = max(meanG2_peaks(locs));
|
||||
MaxPeakAngle = PeakAngles(idxMax);
|
||||
end
|
||||
|
||||
% --- Step 5: Package as struct ---
|
||||
meanFeatures = struct();
|
||||
meanFeatures.MeanContrast = ContrastVal;
|
||||
meanFeatures.MeanPeakAngles = PeakAngles;
|
||||
meanFeatures.MaxPeakAngle = MaxPeakAngle;
|
||||
|
||||
end
|
||||
|
||||
@ -39,11 +39,9 @@ function plotG2Features(results, varargin)
|
||||
if isempty(opts.FigNum), fig = figure; else, fig = figure(opts.FigNum); end
|
||||
clf(fig);
|
||||
set(fig,'Color','w','Position',[100 100 950 750]);
|
||||
t = tiledlayout(3,1,'TileSpacing','Compact','Padding','Compact');
|
||||
t = tiledlayout(2,2,'TileSpacing','Compact','Padding','Compact');
|
||||
|
||||
t = tiledlayout(3,1,'TileSpacing','Compact','Padding','Compact');
|
||||
|
||||
% --- Top tile ---
|
||||
% --- Top-left tile: Fourier amplitudes ---
|
||||
ax1 = nexttile; hold on;
|
||||
errorbar(x, A2m, A2e,'-o','LineWidth',1.5,'Color',colors(1,:),'MarkerFaceColor',colors(1,:));
|
||||
errorbar(x, A4m, A4e,'-s','LineWidth',1.5,'Color',colors(2,:),'MarkerFaceColor',colors(2,:));
|
||||
@ -53,7 +51,7 @@ function plotG2Features(results, varargin)
|
||||
legend(ax1,'|A2|','|A4|','|A6|','Location','best'); grid(ax1,'on');
|
||||
set(ax1,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
|
||||
% --- Middle tile ---
|
||||
% --- Top-right tile: Symmetry fractions ---
|
||||
ax2 = nexttile; hold on;
|
||||
errorbar(x, S2m, S2e,'-o','LineWidth',1.5,'Color',colors(1,:),'MarkerFaceColor',colors(1,:));
|
||||
errorbar(x, Q4m, Q4e,'-s','LineWidth',1.5,'Color',colors(2,:),'MarkerFaceColor',colors(2,:));
|
||||
@ -64,19 +62,45 @@ function plotG2Features(results, varargin)
|
||||
legend(ax2,'S2','Q4','H6','Location','best');
|
||||
set(ax2,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
|
||||
% --- Bottom tile ---
|
||||
ax3 = nexttile; hold on;
|
||||
% --- Bottom-left tile: Contrast (mean of individual curves + mean curve) ---
|
||||
ax3 = nexttile;
|
||||
% Compute mean ± SEM of individual curve contrasts
|
||||
N_params = numel(results.features_group);
|
||||
meanContrast_ind = zeros(1,N_params);
|
||||
semContrast_ind = zeros(1,N_params);
|
||||
for i = 1:N_params
|
||||
FT = results.features_group{i};
|
||||
contrasts = FT.Contrast;
|
||||
meanContrast_ind(i) = mean(contrasts,'omitnan');
|
||||
semContrast_ind(i) = std(contrasts,0,1,'omitnan')/sqrt(numel(contrasts));
|
||||
end
|
||||
errorbar(x, meanContrast_ind, semContrast_ind,'-d','LineWidth',1.5,'Color',[0.2 0.6 0.2],'MarkerFaceColor',[0.2 0.6 0.2]);
|
||||
hold on;
|
||||
% Plot mean-curve contrast on top
|
||||
errorbar(x, Nrm, Nrme,'-d','LineWidth',1.5,'Color',[0.6 0.2 0.8],'MarkerFaceColor',[0.6 0.2 0.8]);
|
||||
ylabel(ax3,'Contrast (RMS)','FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
ylabel(ax3,'Contrast','FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
title(ax3,'g^2 Curve Contrast','FontName',opts.FontName,'FontSize',opts.FontSize+2);
|
||||
grid(ax3,'on'); set(ax3,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
legend(ax3,'Individual Curves Mean ± SEM','Mean Curve','Location','southwest');
|
||||
|
||||
% --- Bottom-right tile: Max-peak location vs scan parameter ---
|
||||
ax4 = nexttile; hold on;
|
||||
% Extract max peak for each scan parameter and convert to degrees
|
||||
maxPeaks = rad2deg(cellfun(@(s) s.MaxPeakAngle, results.meanCurve));
|
||||
plot(x, maxPeaks,'-o','LineWidth',1.5,'Color',[0.8 0.3 0.0],'MarkerFaceColor',[0.8 0.3 0.0]);
|
||||
ylim(ax4, [10, 90])
|
||||
ylabel(ax4,'\theta (degrees)','Interpreter','tex','FontSize',opts.FontSize);
|
||||
xlabel(ax4,opts.XLabel,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
title(ax4,'Max peak in mean g^2 curve','FontName',opts.FontName,'FontSize',opts.FontSize+2);
|
||||
grid(ax4,'on'); set(ax4,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
|
||||
|
||||
% --- Global x-label ---
|
||||
% --- Global x-label for entire figure ---
|
||||
xlabel(t,opts.XLabel,'FontName',opts.FontName,'FontSize',opts.FontSize);
|
||||
|
||||
|
||||
%% --- Save figure ---
|
||||
if ~opts.SkipSaveFigures
|
||||
if ~exist(opts.SaveDirectory,'dir'), mkdir(opts.SaveDirectory); end
|
||||
savefig(fig, fullfile(opts.SaveDirectory, opts.SaveFileName));
|
||||
end
|
||||
end
|
||||
end
|
||||
@ -12,7 +12,7 @@ if useLocalBaseDir
|
||||
baseDir = fullfile(thisScriptDir, 'Results');
|
||||
else
|
||||
% Use manually specified folder
|
||||
baseDir = 'E:\Results - Experiment\202507\BECToDroplets\';
|
||||
baseDir = 'E:\Results - Experiment\202508\BECToDroplets\';
|
||||
end
|
||||
|
||||
% --- Build paths ---
|
||||
|
||||
@ -11,7 +11,7 @@ options = struct();
|
||||
|
||||
% File paths
|
||||
options.baseDataFolder = '//DyLabNAS/Data';
|
||||
options.FullODImagesFolder = 'E:/Data - Experiment/FullODImages/202507';
|
||||
options.FullODImagesFolder = 'E:/Data - Experiment/FullODImages/202508';
|
||||
options.measurementName = 'BECToDroplets';
|
||||
scriptFullPath = mfilename('fullpath');
|
||||
options.saveDirectory = fileparts(scriptFullPath);
|
||||
|
||||
@ -12,7 +12,7 @@ if useLocalBaseDir
|
||||
baseDir = fullfile(thisScriptDir, 'Results');
|
||||
else
|
||||
% Use manually specified folder
|
||||
baseDir = 'E:\Results - Experiment\202507\BECToStripes\';
|
||||
baseDir = 'E:\Results - Experiment\202508\BECToStripes\';
|
||||
end
|
||||
|
||||
% --- Build paths ---
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user