New analysis routine added - extracts anisotropy index from real space g2 matrices.
This commit is contained in:
parent
fd42630628
commit
c012497eef
162
Data-Analyzer/+Analyzer/analyzeG2Structures.m
Normal file
162
Data-Analyzer/+Analyzer/analyzeG2Structures.m
Normal file
@ -0,0 +1,162 @@
|
||||
function analysis_results = analyzeG2Structures(g2_results, opts)
|
||||
%% analyzeG2Structures
|
||||
% Computes peak anisotropy of g² correlation matrices using ROI and boundary
|
||||
% analysis. Anisotropy is zero if the region cannot be well approximated
|
||||
% by an ellipse or does not stretch enough within the ROI.
|
||||
%
|
||||
% Inputs:
|
||||
% g2_results : struct with fields
|
||||
% - g2_matrices : cell array of 2D matrices
|
||||
% - dx_phys : cell array of x-coordinates (µm)
|
||||
% - dy_phys : cell array of y-coordinates (µm)
|
||||
% opts : struct with fields
|
||||
% - roi.center : [x0, y0] center of ROI box (µm)
|
||||
% - roi.size : [w, h] width and height (µm)
|
||||
% - roi.angle : rotation angle (radians, CCW)
|
||||
% - threshold : fraction of peak maximum to define core (default 0.5)
|
||||
% - deviationThreshold : RMS deviation threshold to accept ellipse (default 0.23)
|
||||
% - minEllipseFraction : fraction of ROI diagonal that major axis must span (default 0.8)
|
||||
% - angleLimit : required ellipse angle (radians)
|
||||
% - angleTolerance : tolerance around angleLimit (radians)
|
||||
% - skipLivePlot : true/false (default false)
|
||||
%
|
||||
% Outputs:
|
||||
% analysis_results : struct with fields
|
||||
% - anisotropy_vals : anisotropy per image (0 = irregular or too small, 0–1 = elliptical)
|
||||
% - peak_centroid : [x, y] per image
|
||||
% - boundary_coords : [x, y] coordinates of peak boundary
|
||||
% - ellipse_params : [x0, y0, a, b, theta] per image
|
||||
|
||||
if ~isfield(opts, "skipLivePlot"), opts.skipLivePlot = false; end
|
||||
if ~isfield(opts, "threshold"), opts.threshold = 0.5; end
|
||||
if ~isfield(opts, "deviationThreshold"), opts.deviationThreshold = 0.23; end
|
||||
if ~isfield(opts, "minEllipseFraction"), opts.minEllipseFraction = 0.8; end
|
||||
|
||||
N_images = numel(g2_results.g2_matrices);
|
||||
|
||||
analysis_results = struct();
|
||||
analysis_results.anisotropy_vals = zeros(1, N_images);
|
||||
analysis_results.peak_centroid = cell(1, N_images);
|
||||
analysis_results.boundary_coords = cell(1, N_images);
|
||||
analysis_results.ellipse_params = cell(1, N_images);
|
||||
|
||||
% ROI definition
|
||||
x0 = opts.roi.center(1);
|
||||
y0 = opts.roi.center(2);
|
||||
w = opts.roi.size(1);
|
||||
h = opts.roi.size(2);
|
||||
roi_theta = opts.roi.angle;
|
||||
|
||||
% Precompute ROI diagonal for minimum major axis
|
||||
roi_diag = sqrt(w^2 + h^2);
|
||||
minMajorAxis_auto = opts.minEllipseFraction * roi_diag;
|
||||
|
||||
for k = 1:N_images
|
||||
g2_matrix = g2_results.g2_matrices{k};
|
||||
dx_phys = g2_results.dx_phys{k};
|
||||
dy_phys = g2_results.dy_phys{k};
|
||||
|
||||
[X, Y] = meshgrid(dx_phys, dy_phys);
|
||||
|
||||
% ---- Rotated ROI mask ----
|
||||
Xc = X - x0; Yc = Y - y0;
|
||||
Xr = cos(roi_theta)*Xc + sin(roi_theta)*Yc;
|
||||
Yr = -sin(roi_theta)*Xc + cos(roi_theta)*Yc;
|
||||
roi_mask = (abs(Xr) <= w/2) & (abs(Yr) <= h/2);
|
||||
|
||||
% ---- Threshold ROI to define core region ----
|
||||
g2_roi = zeros(size(g2_matrix));
|
||||
g2_roi(roi_mask) = g2_matrix(roi_mask);
|
||||
thresh_val = opts.threshold * max(g2_roi(:));
|
||||
BW = g2_roi >= thresh_val;
|
||||
|
||||
if ~any(BW(:))
|
||||
warning('No peak found in ROI for image %d', k);
|
||||
analysis_results.peak_centroid{k} = [NaN, NaN];
|
||||
analysis_results.anisotropy_vals(k) = NaN;
|
||||
analysis_results.ellipse_params{k} = [NaN, NaN, NaN, NaN, NaN];
|
||||
analysis_results.boundary_coords{k} = [NaN, NaN];
|
||||
continue;
|
||||
end
|
||||
|
||||
%% ---- Boundary coordinates ----
|
||||
B = bwboundaries(BW);
|
||||
boundary = B{1};
|
||||
x_bound = dx_phys(boundary(:,2));
|
||||
y_bound = dy_phys(boundary(:,1));
|
||||
analysis_results.boundary_coords{k} = [x_bound, y_bound];
|
||||
|
||||
%% ---- Ellipse fitting and anisotropy using boundary region ----
|
||||
stats = regionprops(BW, 'Centroid','MajorAxisLength','MinorAxisLength','Orientation');
|
||||
if ~isempty(stats)
|
||||
a = stats(1).MajorAxisLength / 2 * mean(diff(dx_phys));
|
||||
b = stats(1).MinorAxisLength / 2 * mean(diff(dy_phys));
|
||||
theta = -stats(1).Orientation * pi/180;
|
||||
|
||||
x_c = interp1(1:numel(dx_phys), dx_phys, stats(1).Centroid(1), 'linear', 'extrap');
|
||||
y_c = interp1(1:numel(dy_phys), dy_phys, stats(1).Centroid(2), 'linear', 'extrap');
|
||||
|
||||
Xb = x_bound - x_c;
|
||||
Yb = y_bound - y_c;
|
||||
Xrot = cos(theta)*Xb + sin(theta)*Yb;
|
||||
Yrot = -sin(theta)*Xb + cos(theta)*Yb;
|
||||
|
||||
r_norm = (Xrot/a).^2 + (Yrot/b).^2;
|
||||
deviation = sqrt(mean((r_norm - 1).^2));
|
||||
|
||||
% ---- Angle criterion (± allowed) ----
|
||||
allowed = opts.angleLimit;
|
||||
tol = opts.angleTolerance;
|
||||
angle_ok = (abs(theta - allowed) <= tol) || (abs(theta + allowed) <= tol);
|
||||
|
||||
% Only accept ellipse if deviation small, major axis long enough, and angle ok
|
||||
if (deviation > opts.deviationThreshold) || (a < minMajorAxis_auto) || (~angle_ok)
|
||||
anisotropy = 0;
|
||||
style = '--'; % dashed if rejected
|
||||
else
|
||||
anisotropy = 1 - b/a;
|
||||
style = '-'; % solid if accepted
|
||||
end
|
||||
|
||||
analysis_results.peak_centroid{k} = [x_c, y_c];
|
||||
analysis_results.anisotropy_vals(k) = anisotropy;
|
||||
analysis_results.ellipse_params{k} = [x_c, y_c, a, b, theta];
|
||||
else
|
||||
x_c = NaN; y_c = NaN; a=NaN; b=NaN; theta=NaN; anisotropy=NaN; style='--';
|
||||
analysis_results.peak_centroid{k} = [NaN, NaN];
|
||||
analysis_results.anisotropy_vals(k) = NaN;
|
||||
analysis_results.ellipse_params{k} = [NaN, NaN, NaN, NaN, NaN];
|
||||
end
|
||||
|
||||
%% ---- Visualization ----
|
||||
if ~opts.skipLivePlot
|
||||
fig=figure(100); clf;
|
||||
set(fig, 'Color', 'w', 'Position',[100 100 950 750]);
|
||||
imagesc(dx_phys, dy_phys, g2_matrix);
|
||||
axis image;
|
||||
set(gca, 'YDir', 'normal', 'FontName', opts.font, 'FontSize', 14);
|
||||
colormap(Colormaps.coolwarm()); colorbar;
|
||||
hold on;
|
||||
|
||||
corners = [ -w/2, -h/2; w/2, -h/2; w/2, h/2; -w/2, h/2];
|
||||
R = [cos(roi_theta), -sin(roi_theta); sin(roi_theta), cos(roi_theta)];
|
||||
corners_rot = (R*corners')' + [x0, y0];
|
||||
corners_rot = [corners_rot; corners_rot(1,:)];
|
||||
plot(corners_rot(:,1), corners_rot(:,2), 'r--', 'LineWidth',1.5);
|
||||
|
||||
plot(x_bound, y_bound, 'g-', 'LineWidth',2);
|
||||
plot(x_c, y_c, 'mo', 'MarkerSize',8, 'LineWidth',2);
|
||||
|
||||
t = linspace(0, 2*pi, 200);
|
||||
ellipse_x = x_c + a*cos(t)*cos(theta) - b*sin(t)*sin(theta);
|
||||
ellipse_y = y_c + a*cos(t)*sin(theta) + b*sin(t)*cos(theta);
|
||||
plot(ellipse_x, ellipse_y, ['y' style], 'LineWidth',2);
|
||||
|
||||
title(sprintf('Image %d | Anisotropy = %.2f | Deviation = %.2f | θ = %.2f°', ...
|
||||
k, anisotropy, deviation, rad2deg(theta)));
|
||||
xlabel('x (\mum)'); ylabel('y (\mum)');
|
||||
drawnow; pause(1.0);
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
@ -61,6 +61,9 @@ function results = conductCorrelationAnalysis(od_imgs, scan_parameter_values, op
|
||||
g2_angular = cell(1, N_shots);
|
||||
r_vals = cell(1, N_shots);
|
||||
theta_vals = cell(1, N_shots);
|
||||
dx_phys_ind = cell(1, N_shots);
|
||||
dy_phys_ind = cell(1, N_shots);
|
||||
|
||||
|
||||
dx = pixel_size / magnification;
|
||||
shifts = -maximumShift:maximumShift;
|
||||
@ -79,8 +82,11 @@ function results = conductCorrelationAnalysis(od_imgs, scan_parameter_values, op
|
||||
% Compute g² in physical units
|
||||
[g2_matrix, dx_phys, dy_phys] = Calculator.compute2DAutocorrelation( ...
|
||||
IMG, maximumShift, pixel_size, magnification);
|
||||
|
||||
g2_matrices{k} = g2_matrix;
|
||||
|
||||
dx_phys_ind{k} = dx_phys;
|
||||
dy_phys_ind{k} = dy_phys;
|
||||
|
||||
% Extract radial profile (within angular window)
|
||||
[r_vals{k}, g2_radial{k}] = Calculator.computeRadialCorrelation( ...
|
||||
g2_matrix, dx_phys, dy_phys, radial_theta);
|
||||
@ -182,6 +188,8 @@ function results = conductCorrelationAnalysis(od_imgs, scan_parameter_values, op
|
||||
results.g2_angular = g2_angular;
|
||||
results.r_vals = r_vals;
|
||||
results.theta_vals = theta_vals;
|
||||
results.dx_phys = dx_phys_ind;
|
||||
results.dy_phys = dy_phys_ind;
|
||||
results.scan_parameter_values = unique(scan_parameter_values);
|
||||
|
||||
end
|
||||
|
||||
@ -69,8 +69,8 @@ function plotMeanWithSE(scan_values, data_values, varargin)
|
||||
if ~isempty(opts.YLim)
|
||||
ylim(opts.YLim);
|
||||
end
|
||||
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
|
||||
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
|
||||
xlabel(opts.XLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
|
||||
ylabel(opts.YLabel, 'Interpreter', 'tex', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
|
||||
title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold');
|
||||
grid on;
|
||||
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
function plotPDF(dataCell, referenceValues, varargin)
|
||||
function plotPDF(data_values, scan_reference_values, varargin)
|
||||
%% plotPDF
|
||||
% Author: Karthik
|
||||
% Date: 2025-09-12
|
||||
@ -31,11 +31,11 @@ function plotPDF(dataCell, referenceValues, varargin)
|
||||
parse(p, varargin{:});
|
||||
opts = p.Results;
|
||||
|
||||
N_params = numel(referenceValues);
|
||||
N_params = numel(scan_reference_values);
|
||||
|
||||
% --- Determine y-range ---
|
||||
if isempty(opts.DataRange)
|
||||
allData = cell2mat(dataCell(:));
|
||||
allData = cell2mat(data_values(:));
|
||||
y_min = min(allData);
|
||||
y_max = max(allData);
|
||||
else
|
||||
@ -54,7 +54,7 @@ function plotPDF(dataCell, referenceValues, varargin)
|
||||
|
||||
% --- Compute PDFs ---
|
||||
for i = 1:N_params
|
||||
data = dataCell{i};
|
||||
data = data_values{i};
|
||||
data = data(~isnan(data));
|
||||
if isempty(data), continue; end
|
||||
|
||||
@ -76,14 +76,14 @@ function plotPDF(dataCell, referenceValues, varargin)
|
||||
set(fig, 'Color', 'w', 'Position',[100 100 950 750]);
|
||||
|
||||
if strcmpi(opts.PlotType,'kde')
|
||||
imagesc(referenceValues, y_grid, pdf_matrix);
|
||||
imagesc(scan_reference_values, y_grid, pdf_matrix);
|
||||
else
|
||||
imagesc(referenceValues, binCenters, pdf_matrix);
|
||||
imagesc(scan_reference_values, binCenters, pdf_matrix);
|
||||
end
|
||||
|
||||
set(gca, 'YDir', 'normal', 'FontName', opts.FontName, 'FontSize', opts.FontSize);
|
||||
xlabel(opts.XLabel, 'Interpreter', 'latex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
|
||||
ylabel(opts.YLabel, 'Interpreter', 'latex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
|
||||
xlabel(opts.XLabel, 'Interpreter', 'tex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
|
||||
ylabel(opts.YLabel, 'Interpreter', 'tex', 'FontSize', opts.FontSize, 'FontName', opts.FontName);
|
||||
title(opts.Title, 'FontName', opts.FontName, 'FontSize', opts.FontSize + 2, 'FontWeight', 'bold');
|
||||
colormap(feval(opts.Colormap));
|
||||
c = colorbar;
|
||||
|
||||
@ -99,4 +99,82 @@ end
|
||||
|
||||
%% Conduct correlation analysis
|
||||
|
||||
results = Analyzer.conductCorrelationAnalysis(od_imgs, scan_parameter_values, options);
|
||||
g2_analysis_results = Analyzer.conductCorrelationAnalysis(od_imgs, scan_parameter_values, options);
|
||||
|
||||
%% Analyze G2 matrices
|
||||
|
||||
% ROI definition
|
||||
options.roi.center = [3, 3]; % center of ROI in µm (x0, y0)
|
||||
options.roi.size = [3, 11]; % width and height in µm
|
||||
options.roi.angle = pi/4; % rotation angle (radians, CCW)
|
||||
options.threshold = 0.85;
|
||||
|
||||
options.deviationThreshold = 0.30;
|
||||
options.minEllipseFraction = 0.30;
|
||||
options.angleLimit = deg2rad(45);
|
||||
options.angleTolerance = deg2rad(5);
|
||||
|
||||
|
||||
% Plot control
|
||||
options.skipLivePlot = false; % set true if you don't want per-image figures
|
||||
|
||||
analysis_results = Analyzer.analyzeG2Structures(g2_analysis_results, options);
|
||||
|
||||
%% Plot mean and standard error of anisotropy
|
||||
|
||||
Plotter.plotMeanWithSE(scan_parameter_values, analysis_results.anisotropy_vals, ...
|
||||
'Title', options.titleString, ...
|
||||
'YLim', [0,1], ...
|
||||
'XLabel', '\alpha (degrees)', ...
|
||||
'YLabel', 'Anisotropy of correlation peaks', ...
|
||||
'FigNum', 1, ...
|
||||
'FontName', options.font, ...
|
||||
'SaveFileName', 'RadialSpectralContrast.fig', ...
|
||||
'SaveDirectory', pwd, ... % save figures inside dataset-specific folder
|
||||
'SkipSaveFigures', options.skipSaveFigures);
|
||||
|
||||
%% Plot distribution of anisotropy
|
||||
grouped_data = groupDataByScan(scan_parameter_values, analysis_results.anisotropy_vals);
|
||||
|
||||
% call plotPDF
|
||||
Plotter.plotPDF(grouped_data, ...
|
||||
scan_reference_values, ...
|
||||
'Title', options.titleString, ...
|
||||
'XLabel', '\alpha (degrees)', ...
|
||||
'YLabel', 'Anisotropy of correlation peaks', ...
|
||||
'FigNum', 2, ...
|
||||
'FontName', options.font, ...
|
||||
'SkipSaveFigures', options.skipSaveFigures, ...
|
||||
'SaveFileName', 'PDF_MaxG2AcrossTransition.fig', ...
|
||||
'SaveDirectory', pwd, ...
|
||||
'NumberOfBins', 10, ...
|
||||
'NormalizeHistogram', true, ...
|
||||
'DataRange', [0 1.0], ...
|
||||
'Colormap', @Colormaps.coolwarm, ...
|
||||
'XLim', [min(scan_reference_values) max(scan_reference_values)]);
|
||||
|
||||
function groupedData = groupDataByScan(scan_values, data_values)
|
||||
%% groupByScanValues
|
||||
% Groups data according to unique scan parameter values.
|
||||
%
|
||||
% Inputs:
|
||||
% scan_values : array of scan parameters (length = N_reps * N_scan)
|
||||
% data_values : numeric array or cell array of measured values
|
||||
% (same length as scan_values)
|
||||
%
|
||||
% Output:
|
||||
% groupedData : 1 x N_unique cell array, each containing all repetitions
|
||||
% corresponding to a unique scan value
|
||||
|
||||
[unique_vals, ~, idx] = unique(scan_values, 'stable'); % preserve order
|
||||
groupedData = cell(1, numel(unique_vals));
|
||||
|
||||
for i = 1:numel(unique_vals)
|
||||
if iscell(data_values)
|
||||
group = data_values(idx == i);
|
||||
groupedData{i} = [group{:}]; % concatenate if nested cells
|
||||
else
|
||||
groupedData{i} = data_values(idx == i);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user