function analysis_results = analyzeG2Structures(g2_results, opts)
%% analyzeG2Structures
% Computes peak anisotropy of g² correlation matrices using ROI and boundary
% analysis, and fits an oriented two-mode Gaussian aligned with the
% structure's principal axis 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)
%       - skipLivePlot : true/false (default false)
%       - fitDeviationThreshold : max normalized RMS deviation for accepting fit
%       - boundaryPad : extra pixels around boundary to allow Gaussian peaks
%
% Outputs:
%   analysis_results : struct with fields
%       - anisotropy_vals : numeric array
%       - peak_centroid   : {N×1} fitted centroids [x, y]
%       - ellipse_params  : {N×1} fitted parameters + angle
%       - boundary_coords : {N×1} [x, y] boundary points

if ~isfield(opts, "skipLivePlot"), opts.skipLivePlot = false; end
if ~isfield(opts, "threshold"), opts.threshold = 0.5; end
if ~isfield(opts, "font"), opts.font = 'Arial'; end
if ~isfield(opts, "fitDeviationThreshold"), opts.fitDeviationThreshold = 0.2; end
if ~isfield(opts, "boundaryPad"), opts.boundaryPad = 2; 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;

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);

    %% ---- Core ROI extraction (rotated rectangle) ----
    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);

    x_roi = Xr(roi_mask);
    y_roi = Yr(roi_mask);
    z_roi = g2_matrix(roi_mask);

    %% ---- 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(1,7);
        analysis_results.boundary_coords{k} = [NaN, NaN];
        continue;
    end

    %% ---- Boundary coordinates (edge detection inside ROI) ----
    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];

    %% ---- Compute centroid and orientation (use PCA on boundary in ROI-local frame) ----
    % Transform boundary points to ROI-local coordinates (same local coords as Xr,Yr)
    Xc_b = x_bound - x0;
    Yc_b = y_bound - y0;
    Xr_b =  cos(roi_theta)*Xc_b + sin(roi_theta)*Yc_b;
    Yr_b = -sin(roi_theta)*Xc_b + cos(roi_theta)*Yc_b;

    % PCA on boundary points in ROI-local frame (robust and simple)
    boundary_mean = [mean(Xr_b), mean(Yr_b)];
    XYb = [Xr_b - boundary_mean(1), Yr_b - boundary_mean(2)];
    try
        coeff = pca(XYb);
        major = coeff(:,1);
        angle_region = atan2(major(2), major(1)); % angle in ROI-local coords (radians)
    catch
        % fallback: use roi_theta as orientation if PCA fails
        angle_region = 0;
    end

    % We'll define longitudinal axis in ROI-local coordinates as angle_region,
    % transverse axis is angle_region + pi/2.

    %% ---- Use points inside boundary for 1D profile (use ROI-local coords) ----
    % We already have axis-aligned ROI samples x_roi, y_roi, z_roi in ROI-local coords.
    xData_local = x_roi;
    yData_local = y_roi;
    zData       = z_roi;
    
    %% ---- Interpolate onto regular axis-aligned grid in ROI-local coords ----
    Nx = ceil(w*10); Ny = ceil(h*10);
    xq = linspace(-w/2, w/2, Nx);
    yq = linspace(-h/2, h/2, Ny);
    [Xq, Yq] = meshgrid(xq, yq);
    Zq = griddata(xData_local, yData_local, zData, Xq, Yq, 'cubic');

    % If griddata produced NaNs in large patches, fill with nearest to avoid empty bins
    if any(isnan(Zq(:)))
        Zq = fillmissing(Zq,'nearest');
    end

    % --- Build rotated coordinates so major axis aligns with new X'
    theta = angle_region; % ROI-local angle of major axis
    Rrot = [cos(theta), sin(theta); -sin(theta), cos(theta)]; % rotation that maps [x;y] -> [x';y'] where x' along major
    XY = [Xq(:)'; Yq(:)'];
    XY_rot = Rrot * XY;
    Xrot = reshape(XY_rot(1,:), size(Xq));
    Yrot = reshape(XY_rot(2,:), size(Yq));

    % --- Compute 1D longitudinal profile by binning along X' (major axis)
    nbins_long = size(Xq,1); % choose number of bins ~ Ny (vertical resolution)
    x_edges_long = linspace(min(Xrot(:)), max(Xrot(:)), nbins_long+1);
    x_centers_long = 0.5*(x_edges_long(1:end-1)+x_edges_long(2:end));
    prof_long = nan(1, nbins_long);
    for ii = 1:nbins_long
        mask_bin = Xrot >= x_edges_long(ii) & Xrot < x_edges_long(ii+1);
        vals = Zq(mask_bin);
        if isempty(vals)
            prof_long(ii) = NaN;
        else
            prof_long(ii) = mean(vals,'omitnan');
        end
    end

    % --- Compute 1D transverse profile by binning along Y' (minor axis)
    nbins_trans = size(Xq,2);
    y_edges_trans = linspace(min(Yrot(:)), max(Yrot(:)), nbins_trans+1);
    y_centers_trans = 0.5*(y_edges_trans(1:end-1)+y_edges_trans(2:end));
    prof_trans = nan(1, nbins_trans);
    for ii = 1:nbins_trans
        mask_bin = Yrot >= y_edges_trans(ii) & Yrot < y_edges_trans(ii+1);
        vals = Zq(mask_bin);
        if isempty(vals)
            prof_trans(ii) = NaN;
        else
            prof_trans(ii) = mean(vals,'omitnan');
        end
    end

    % --- Prepare x-axis positions for profile centers (in ROI-local units)
    xcent_long = x_centers_long;      % positions along major axis (µm)
    xcent_trans = y_centers_trans;    % positions along minor axis (µm)

    % --- Remove NaNs from profiles
    valid_long = ~isnan(prof_long);
    x_long = xcent_long(valid_long);
    y_long = prof_long(valid_long);

    valid_trans = ~isnan(prof_trans);
    x_trans = xcent_trans(valid_trans);
    y_trans = prof_trans(valid_trans);

    % --- If not enough points, skip
    if numel(x_long) < 8 || numel(x_trans) < 8
        warning('Too few points for longitudinal/transverse fits for image %d', k);
        analysis_results.peak_centroid{k} = [NaN, NaN];
        analysis_results.anisotropy_vals(k) = NaN;
        analysis_results.ellipse_params{k} = nan(1,7);
        continue;
    end

    %% ---- Fit two-Gaussian to longitudinal profile (use your working 1D fit) ----
    twoGauss1D = @(p,x) ...
        p(1) * exp(-0.5*((x - p(2))/max(p(3),1e-6)).^2) + ... % Gaussian 1
        p(4) * exp(-0.5*((x - p(5))/max(p(6),1e-6)).^2) + ... % Gaussian 2
        p(7);                                                % offset
    % initial guesses (longitudinal)
    [~, mainIdxL] = max(y_long);
    A1_guessL     = y_long(mainIdxL);
    mu1_guessL    = x_long(mainIdxL);
    sigma1_guessL = range(x_long)/10;
    A2_guessL     = 0.8 * A1_guessL;
    mu2_guessL    = mu1_guessL + range(x_long)/5;
    sigma2_guessL = sigma1_guessL;
    offset_guessL = min(y_long);
    p0L = [A1_guessL, mu1_guessL, sigma1_guessL, A2_guessL, mu2_guessL, sigma2_guessL, offset_guessL];
    lbL = [0, min(x_long), 1e-6, 0, min(x_long), 1e-6, -Inf];
    ubL = [2, max(x_long), range(x_long), 2, max(x_long), range(x_long), Inf];
    opts_lsq = optimoptions('lsqcurvefit','Display','off','MaxFunctionEvaluations',1e4,'MaxIterations',1e4);
    try
        pFitL = lsqcurvefit(twoGauss1D, p0L, x_long, y_long, lbL, ubL, opts_lsq);
    catch
        warning('Longitudinal 1D fit failed for image %d, using initial guess.', k);
        pFitL = p0L;
    end

    %% ---- Fit two-Gaussian to transverse profile (use same logic) ----
    [~, mainIdxT] = max(y_trans);
    A1_guessT     = y_trans(mainIdxT);
    mu1_guessT    = x_trans(mainIdxT);
    sigma1_guessT = range(x_trans)/10;
    A2_guessT     = 0.8 * A1_guessT;
    mu2_guessT    = mu1_guessT + range(x_trans)/5;
    sigma2_guessT = sigma1_guessT;
    offset_guessT = min(y_trans);
    p0T = [A1_guessT, mu1_guessT, sigma1_guessT, A2_guessT, mu2_guessT, sigma2_guessT, offset_guessT];
    lbT = [0, min(x_trans), 1e-6, 0, min(x_trans), 1e-6, -Inf];
    ubT = [2, max(x_trans), range(x_trans), 2, max(x_trans), range(x_trans), Inf];
    try
        pFitT = lsqcurvefit(twoGauss1D, p0T, x_trans, y_trans, lbT, ubT, opts_lsq);
    catch
        warning('Transverse 1D fit failed for image %d, using initial guess.', k);
        pFitT = p0T;
    end

    %% ---- Determine anisotropy from the two fits (ratio of widths) ----
    sigma_long = mean([pFitL(3), pFitL(6)]);
    sigma_trans = mean([pFitT(3), pFitT(6)]);
    anisotropy = sigma_trans / sigma_long; % transverse / longitudinal (ratio)
    % store centroid roughly as ROI-local centroid rotated back to lab coords
    centroid_local_roi = boundary_mean'; % mean of boundary in ROI-local coords
    % rotate back to lab frame
    rotBack = [cos(roi_theta), -sin(roi_theta); sin(roi_theta), cos(roi_theta)];
    centroid_lab = rotBack * centroid_local_roi + [x0; y0];
    centroid_lab = centroid_lab(:)';

    analysis_results.peak_centroid{k} = centroid_lab;
    analysis_results.anisotropy_vals(k) = anisotropy;
    analysis_results.ellipse_params{k} = [pFitL, pFitT, angle_region];

    %% ---- Visualization ----
    if ~opts.skipLivePlot
        fig = figure(100); clf;
        set(fig,'Color','w','Position',[100 100 1600 450]);
        tiledlayout(1,3,'Padding','compact','TileSpacing','compact');

        % --- Panel 1: 2D g² map + ROI + boundary
        nexttile;
        imagesc(dx_phys, dy_phys, g2_matrix);
        axis image; set(gca,'YDir','normal'); colormap(Colormaps.coolwarm()); colorbar; hold on;

        % ROI bounding box (red dashed)
        rect = [-w/2, -h/2; 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)];
        rect_rot = (R * rect')';
        plot(rect_rot(:,1)+x0, rect_rot(:,2)+y0, 'r--', 'LineWidth',2);

        % Detected peak boundary (green solid)
        if ~isempty(x_bound)
            plot(x_bound, y_bound, 'g-', 'LineWidth',2);
        end

        xlabel('\Deltax (\mum)'); ylabel('\Deltay (\mum)');
        title(sprintf('Image %d | 2D g^2', k));

        % --- Panel 2: Longitudinal 1D profile and fit
        nexttile;
        plot(x_long, y_long, 'k.-', 'DisplayName','Data'); hold on;
        xFine = linspace(min(x_long), max(x_long), 500);
        yFitL = twoGauss1D(pFitL, xFine);
        plot(xFine, yFitL, 'r-', 'LineWidth',1.5,'DisplayName','Fit');
        title('Longitudinal profile','FontName',opts.font);
        xlabel('x_{longitudinal} (\mum)');
        ylabel('g^2');
        legend('Location','northeast'); grid on;

        % --- Panel 3: Transverse 1D profile and fit
        nexttile;
        plot(x_trans, y_trans, 'b.-', 'DisplayName','Data'); hold on;
        xFineT = linspace(min(x_trans), max(x_trans), 500);
        yFitT = twoGauss1D(pFitT, xFineT);
        plot(xFineT, yFitT, 'r-', 'LineWidth',1.5,'DisplayName','Fit');
        title('Transverse profile','FontName',opts.font);
        xlabel('x_{transverse} (\mum)');
        ylabel('g^2');
        legend('Location','northeast'); grid on;

        drawnow;
    end
end

end