%% Extract Images

baseDir       = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/DTS/Point2';
JobNumber     = 0;
runFolder     = sprintf('Run_%03d', JobNumber);
savefileName  = 'DropletsToStripes';  % Output file name
datafileName  = './DropletsToStripes.mat';
reverseOrder  = false;  % Set this to true to reverse the theta ordering

TitleString  = 'Droplets To Stripes';

%%
baseDir       = 'D:/Results - Numerics/Data_Full3D/PhaseTransition/STD/Point2';
JobNumber     = 0;
runFolder     = sprintf('Run_%03d', JobNumber);
savefileName  = 'StripesToDroplets';  % Output file name
datafileName  = './StripesToDroplets.mat';
reverseOrder  = true;  % Set this to true to reverse the theta ordering


TitleString  = 'Stripes To Droplets';

%%
folderList = dir(baseDir);
isValid = [folderList.isdir] & ~ismember({folderList.name}, {'.', '..'});
folderNames = {folderList(isValid).name};
nimgs = numel(folderNames);

% Extract theta values from folder names
PolarAngleVals = zeros(1, nimgs);
for k = 1:nimgs
    tokens = regexp(folderNames{k}, 'theta_(\d{3})', 'tokens');
    if isempty(tokens)
        warning('No theta found in folder name: %s', folderNames{k});
        PolarAngleVals(k) = NaN;
    else
        PolarAngleVals(k) = str2double(tokens{1}{1});
    end
end

% Choose sort direction
sortDirection = 'ascend';
if reverseOrder
    sortDirection = 'descend';
end

% Sort folderNames based on polar angle
[~, sortIdx]     = sort(PolarAngleVals, sortDirection);
folderNames      = folderNames(sortIdx);
PolarAngleVals   = PolarAngleVals(sortIdx);  % Optional: if you still want sorted list

imgs             = cell(1, nimgs);
alphas           = zeros(1, nimgs);

for k = 1:nimgs
    folderName = folderNames{k};
    SaveDirectory = fullfile(baseDir, folderName, runFolder);

    % Extract alpha (theta) again from folder name
    tokens = regexp(folderName, 'theta_(\d{3})', 'tokens');
    alpha_val = str2double(tokens{1}{1});
    alphas(k) = alpha_val;

    matPath = fullfile(SaveDirectory, 'psi_gs.mat');
    if ~isfile(matPath)
        warning('Missing psi_gs.mat in %s', SaveDirectory);
        continue;
    end

    try
        Data = load(matPath, 'psi', 'Params', 'Transf', 'Observ');
    catch ME
        warning('Failed to load %s: %s', matPath, ME.message);
        continue;
    end

    Params = Data.Params;
    Transf = Data.Transf;
    Observ = Data.Observ;

    psi = Data.psi;
    if isgpuarray(psi)
        psi = gather(psi);
    end
    if isgpuarray(Observ.residual)
        Observ.residual = gather(Observ.residual);
    end
    
    % Axes and projection
    z = Transf.z * Params.l0 * 1e6;
    dz = z(2)-z(1);
    
    if k == 1
        x = Transf.x * Params.l0 * 1e6;
        y = Transf.y * Params.l0 * 1e6;
        % Calculate frequency increment (frequency axes)
        Nx           = length(x);       % grid size along X 
        Ny           = length(y);       % grid size along Y 
        dx           = mean(diff(x));   % real space increment in the X direction (in micrometers)
        dy           = mean(diff(y));   % real space increment in the Y direction (in micrometers)
        dvx          = 1 / (Nx * dx);  % reciprocal space increment in the X direction (in micrometers^-1)
        dvy          = 1 / (Ny * dy);  % reciprocal space increment in the Y direction (in micrometers^-1)
        
        % Create the frequency axes
        vx           = (-Nx/2:Nx/2-1) * dvx;  % Frequency axis in X (micrometers^-1)
        vy           = (-Ny/2:Ny/2-1) * dvy;  % Frequency axis in Y (micrometers^-1)
           
        % Create the Wavenumber axes
        kx_full          = 2*pi*vx;  % Wavenumber axis in X
        ky_full          = 2*pi*vy;  % Wavenumber axis in Y
    end

    n = abs(psi).^2;
    nxy = squeeze(trapz(n * dz, 3));

    imgs{k} = nxy;
end

%% Analyze Images

skipMovieRender             = true; % Set to true to disable movie creation
skipSaveFigures             = false;
font                        = 'Bahnschrift';

skipPreprocessing           = true;
skipMasking                 = true;
skipIntensityThresholding   = true;
skipBinarization            = true;

% Run Fourier analysis over images

fft_imgs                       = cell(1, nimgs);
angular_spectral_distribution  = cell(1, nimgs);
theta_values              = cell(1, nimgs);
radial_spectral_contrast  = zeros(1, nimgs);
angular_spectral_weight   = zeros(1, nimgs);
g2_all                    = cell(1, nimgs);
theta_values_all          = cell(1, nimgs);

N_shots                   = length(imgs);
ps_list                   = cell(1, N_shots);  % 2D power spectrum
s_k_list                  = cell(1, N_shots);  % Radial spectrum
s_theta_list              = cell(1, N_shots);  % Angular spectrum

% Fourier analysis settings

% Radial Spectral Distribution
theta_min          = deg2rad(0);
theta_max          = deg2rad(180);
N_radial_bins      = 500;
Radial_Sigma       = 2;
Radial_WindowSize  = 5; % Choose an odd number for a centered moving average

% Angular Spectral Distribution
k_min              = 1.2771;    % in μm⁻¹
k_max              = 2.5541;    % in μm⁻¹
N_angular_bins     = 180;
Angular_Threshold  = 25;
Angular_Sigma      = 2;
Angular_WindowSize = 5;

zoom_size          = 50;             % Zoomed-in region around center

if ~skipMovieRender
    % Create VideoWriter object for movie
    videoFile           = VideoWriter([savefileName '.mp4'], 'MPEG-4');
    videoFile.Quality   = 100; % Set quality to maximum (0–100)
    videoFile.FrameRate = 2;  % Set the frame rate (frames per second)
    open(videoFile);  % Open the video file to write
end

if ~skipSaveFigures
    % Define folder for saving images
    saveFolder = [savefileName '_SavedFigures'];
    if ~exist(saveFolder, 'dir')
        mkdir(saveFolder);
    end
end
% Display the cropped image
for k = 1:nimgs
    IMG                   = imgs{k};
    if ~(max(IMG(:)) > 1)
        IMGFFT            = NaN(size(IMG));
    else
        [IMGFFT, IMGPR]   = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
    end
    
    % Crop FFT image around center
    mid_x                 = floor(Nx/2);
    mid_y                 = floor(Ny/2);
    fft_imgs{k}           = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
    
    % Crop wavenumber axes to match fft_imgs{k}
    kx                    = kx_full(mid_x - zoom_size : mid_x + zoom_size);
    ky                    = ky_full(mid_y - zoom_size : mid_y + zoom_size);
    
    [theta_vals, S_theta]        = computeAngularSpectralDistribution(fft_imgs{k}, kx, ky, k_min, k_max, N_angular_bins, Angular_Threshold, Angular_Sigma, []);
    [k_rho_vals, S_k]            = computeRadialSpectralDistribution(fft_imgs{k}, kx, ky, theta_min, theta_max, N_radial_bins);
    S_k_smoothed                 = movmean(S_k, Radial_WindowSize);  % Compute moving average (use convolution) or use conv for more control
    angular_spectral_distribution{k}  = S_theta;
    theta_values{k}                   = theta_vals;
    radial_spectral_contrast(k)  = computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max);
    S_theta_norm                 = S_theta / max(S_theta); % Normalize to 1
    angular_spectral_weight(k)   = trapz(theta_vals, S_theta_norm);

    g2 = zeros(1, N_angular_bins);  % Preallocate
    
    for dtheta = 0:N_angular_bins-1
        profile         = S_theta;
        profile_shifted = circshift(profile, -dtheta, 2);
    
        num   = mean(profile .* profile_shifted);
        denom = mean(profile.^2);
    
        g2(dtheta+1) = num / denom;
    end
    
    ps_list{k}              = abs(fft_imgs{k}).^2;  % store the power spectrum
    s_k_list{k}             = S_k_smoothed;         % store smoothed radial spectrum
    s_theta_list{k}         = S_theta;              % store angular spectrum
    g2_all{k}               = g2;
    theta_values_all{k}     = theta_vals;

    figure(1);
    clf
    set(gcf,'Position',[500 100 1000 800])
    t = tiledlayout(2, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 1x4 grid

    y_min = min(y);
    y_max = max(y);
    x_min = min(x);
    x_max = max(x);
    

    % Display the cropped OD image
    ax1 = nexttile;
    imagesc(x, y, IMG')
    % Define normalized positions (relative to axis limits)
    x_offset = 0.025; % 5% offset from the edges
    y_offset = 0.025; % 5% offset from the edges
    % Top-right corner (normalized axis coordinates)
    hText = text(1 - x_offset, 1 - y_offset, ['Angle = ', num2str(alphas(k), '%.1f')], ...
        'Color', 'white', 'FontWeight', 'bold', 'Interpreter', 'tex', 'FontSize', 20, 'Units', 'normalized', 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top');
    axis square;
    hcb = colorbar;
    colormap(ax1, Colormaps.plasma());
    set(gca, 'FontSize', 14);   % For tick labels only
    hL = ylabel(hcb, 'Optical Density');
    set(hL,'Rotation',-90);
    set(gca,'YDir','normal')
    % set(gca, 'YTick', linspace(y_min, y_max, 5));  % Define y ticks
    % set(gca, 'YTickLabel', flip(linspace(y_min, y_max, 5)));  % Flip only the labels
    hXLabel = xlabel('x (pixels)', 'Interpreter', 'tex');
    hYLabel = ylabel('y (pixels)', 'Interpreter', 'tex');
    hTitle = title('Density', 'Interpreter', 'tex');
    set([hXLabel, hYLabel, hL, hText], 'FontName', font)
    set([hXLabel, hYLabel, hL], 'FontSize', 14)
    set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
    
    % ======= FFT POWER SPECTRUM (reciprocal space) =======
    ax2 = nexttile;
    imagesc(kx, ky, log(1 + ps_list{k}));
    axis tight;
    set(gca, 'FontSize', 14, 'YDir', 'normal')
    xlabel('k_x [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
    ylabel('k_y [\mum^{-1}]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
    title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex', ...
          'FontSize', 16, 'FontWeight', 'bold', 'FontName', font);
    colorbar;
    colormap(ax2, Colormaps.coolwarm());
    
    % ======= ANGULAR DISTRIBUTION (S(θ)) =======
    nexttile;
    plot(theta_vals/pi, S_theta_norm, 'LineWidth', 2);
    axis image;
    set(gca, 'FontSize', 14, 'YLim', [0, 1]);
    xlabel('\theta/\pi [rad]', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
    ylabel('Magnitude (a.u.)', 'Interpreter', 'tex', 'FontSize', 14, 'FontName', font);
    title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex', ...
        'FontSize', 16, 'FontWeight', 'bold', 'FontName', font);
    grid on;            % Enable major grid
    ax = gca;
    ax.MinorGridLineStyle = ':';  % Optional: make minor grid dotted
    ax.MinorGridColor = [0.7 0.7 0.7]; % Optional: light gray minor grid color
    ax.MinorGridAlpha = 0.5;      % Optional: transparency for minor grid
    
    ax.XMinorGrid = 'on';         % Enable minor grid for x-axis
    ax.YMinorGrid = 'on';         % Enable minor grid for y-axis (if desired)
    
    % ======= ANGULAR CORRELATION g2(θ) =======
    nexttile
    plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5);
    set(gca, 'FontSize', 14);
    ylim([0.0 1.0]);  % Set y-axis limits here
    hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
    hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
    hTitle  = title('Angular Correlation', 'Interpreter', 'tex');
    set([hXLabel, hYLabel], 'FontName', font)
    set([hXLabel, hYLabel], 'FontSize', 14)
    set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
    grid on;
    
    if ~skipMovieRender
        % Capture the current frame and write it to the video
        frame = getframe(gcf);  % Capture the current figure as a frame
        writeVideo(videoFile, frame);  % Write the frame to the video
    end
    if ~skipSaveFigures
        % Construct a filename for each image
        fileNamePNG = fullfile(saveFolder, sprintf('fft_analysis_img_%03d.png', k));
        
        % Save current figure as PNG with high resolution
        print(gcf, fileNamePNG, '-dpng', '-r100');  % 300 dpi for high quality
    end
    if skipMovieRender & skipSaveFigures
        pause(0.5);
    end
end

if ~skipMovieRender
    % Close the video file
    close(videoFile);
    disp(['Movie saved to ', savefileName]);
end

%% Track across the transition

% Convert spectral distribution to matrix (N_shots x N_angular_bins)
delta_nkr_all = zeros(N_shots, N_angular_bins);
for k = 1:N_shots
    delta_nkr_all(k, :) = angular_spectral_distribution{k};
end

N_params = length(alphas);
% Define angular range and conversion
angle_range     = 180;
angle_per_bin   = angle_range / N_angular_bins;
max_peak_angle  = 180;
max_peak_bin    = round(max_peak_angle / angle_per_bin);

% Parameters for search
window_size     = 10;
angle_threshold = 100;

% Initialize containers for final results
mean_max_g2_values       = zeros(1, N_params);
mean_max_g2_angle_values = zeros(1, N_params);
var_max_g2_values        = zeros(1, N_params);
var_max_g2_angle_values  = zeros(1, N_params);

% Also store raw data per group
g2_all_per_group         = cell(1, N_params);
angle_all_per_group      = cell(1, N_params);

for i = 1:N_params
    group_data     = delta_nkr_all(i, :);
    N_reps         = size(group_data, 1);

    g2_values         = zeros(1, N_reps);
    angle_at_max_g2   = zeros(1, N_reps);

    for j = 1:N_reps
        profile                = group_data(j, :);

        % Restrict search to 0–60° for highest peak
        restricted_profile     = profile(1:max_peak_bin);
        [~, peak_idx_rel]      = max(restricted_profile);
        peak_idx               = peak_idx_rel;
        peak_angle             = (peak_idx - 1) * angle_per_bin;

        if peak_angle < angle_threshold
            offsets            = round(50 / angle_per_bin) : round(70 / angle_per_bin);
        else
            offsets            = -round(70 / angle_per_bin) : -round(50 / angle_per_bin);
        end

        ref_window             = mod((peak_idx - window_size):(peak_idx + window_size) - 1, N_angular_bins) + 1;
        ref                    = profile(ref_window);

        correlations           = zeros(size(offsets));
        angles                 = zeros(size(offsets));

        for k = 1:length(offsets)
            shifted_idx        = mod(peak_idx + offsets(k) - 1, N_angular_bins) + 1;
            sec_window         = mod((shifted_idx - window_size):(shifted_idx + window_size) - 1, N_angular_bins) + 1;
            sec                = profile(sec_window);

            num                = mean(ref .* sec);
            denom              = mean(ref.^2);
            g2                 = num / denom;

            correlations(k)    = g2;
            angles(k)          = mod((peak_idx - 1 + offsets(k)) * angle_per_bin, angle_range);
        end

        [max_corr, max_idx]    = max(correlations);
        g2_values(j)           = max_corr;
        angle_at_max_g2(j)     = angles(max_idx);
    end

    % Store raw values
    g2_all_per_group{i}      = g2_values;
    angle_all_per_group{i}   = angle_at_max_g2;

    % Final stats
    mean_max_g2_values(i)      = mean(g2_values, 'omitnan');
    var_max_g2_values(i)       = var(g2_values, 0, 'omitnan');
    mean_max_g2_angle_values(i)= mean(angle_at_max_g2, 'omitnan');
    var_max_g2_angle_values(i) = var(angle_at_max_g2, 0, 'omitnan');
end

% Plot peak offset angular correlation
figure(2);
set(gcf,'Position',[100 100 950 750])
set(gca, 'FontSize', 14);    % For tick labels only
plot(alphas, ...             % x-axis
     mean_max_g2_values, ... % y-axis (mean)
     '--o', 'LineWidth', 1.8, 'MarkerSize', 6 );

set(gca, 'FontSize', 14, 'YLim', [0, 1]);
hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
hTitle  = title(TitleString, 'Interpreter', 'tex');
set([hXLabel, hYLabel], 'FontName', font);
set([hXLabel, hYLabel], 'FontSize', 14);
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');
grid on;

% Plot full angular correlation
figure(3);
clf;
set(gcf,'Position',[100 100 950 750])
hold on;

% Reconstruct theta axis from any one of the stored values
theta_vals = theta_values_all{1};  % assuming it's in radians

legend_entries = cell(nimgs, 1);

% Generate a colormap with enough unique colors
cmap = sky(nimgs);  % You can also try 'jet', 'turbo', 'hot', etc.

for i = 1:nimgs
    plot(theta_vals/pi, g2_all{i}, ...
        'o-', 'Color', cmap(i,:), 'LineWidth', 1.2, ...
        'MarkerSize', 5);
    legend_entries{i} = sprintf('$\\alpha = %g^\\circ$', alphas(i));
end

ylim([0.0 1.0]);  % Set y-axis limits here
set(gca, 'FontSize', 14);
hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
hTitle  = title(TitleString, 'Interpreter', 'tex');
legend(legend_entries, 'Interpreter', 'latex', 'Location', 'bestoutside');
set([hXLabel, hYLabel], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
grid on;

save(datafileName, 'alphas', 'mean_max_g2_values', 'theta_vals', 'g2_all');

%% Compare between transitions

set(0,'defaulttextInterpreter','latex')
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');

format long

font                       = 'Bahnschrift';

% Load data
Data                       = load('./DropletsToStripes.mat', 'alphas', 'mean_max_g2_values');
dts_alphas                 = Data.alphas;
dts_max_g2                 = Data.mean_max_g2_values;

Data                       = load('./StripesToDroplets.mat', 'alphas', 'mean_max_g2_values');
std_alphas                 = Data.alphas;
std_max_g2                 = Data.mean_max_g2_values;

figure(5);
set(gcf,'Position',[100 100 950 750])
plot(dts_alphas, dts_max_g2, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes');
hold on
plot(std_alphas, std_max_g2, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets');
set(gca, 'FontSize', 14);   % For tick labels only
hXLabel = xlabel('$\alpha$ (degrees)', 'Interpreter', 'latex');
hYLabel = ylabel('$\mathrm{max}[g^{(2)}_{[50,70]}(\delta\theta)]$', 'Interpreter', 'latex');
hTitle = title('Change across transition', 'Interpreter', 'tex');
legend
set([hXLabel, hYLabel], 'FontName', font)
set([hXLabel, hYLabel], 'FontSize', 14)
set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
grid on

%%
function [IMGFFT, IMGPR] = computeFourierTransform(I, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization)
    % computeFourierSpectrum - Computes the 2D Fourier power spectrum 
    % of binarized and enhanced lattice image features, with optional central mask.
    %
    % Inputs:
    %   I           - Grayscale or RGB image matrix
    %
    % Output:
    %   F_mag       - 2D Fourier power spectrum (shifted)
    
    if ~skipPreprocessing
        % Preprocessing: Denoise
        filtered             = imgaussfilt(I, 10);
        IMGPR                = I - filtered; % adjust sigma as needed
    else
        IMGPR                = I;
    end
    
    if ~skipMasking
        [rows, cols]     = size(IMGPR);
        [X, Y]           = meshgrid(1:cols, 1:rows);
        % Elliptical mask parameters
        cx               = cols / 2;
        cy               = rows / 2;
        
        % Shifted coordinates
        x                = X - cx;
        y                = Y - cy;
        
        % Ellipse semi-axes
        rx               = 0.4 * cols;
        ry               = 0.2 * rows;
        
        % Rotation angle in degrees -> radians
        theta_deg        = 30;        % Adjust as needed
        theta            = deg2rad(theta_deg);
        
        % Rotated ellipse equation
        cos_t            = cos(theta);
        sin_t            = sin(theta);
        
        x_rot            = (x * cos_t + y * sin_t);
        y_rot            = (-x * sin_t + y * cos_t);
        
        ellipseMask      = (x_rot.^2) / rx^2 + (y_rot.^2) / ry^2 <= 1;
        
        % Apply cutout mask
        IMGPR            = IMGPR .* ellipseMask;
    end

    if ~skipIntensityThresholding
        % Apply global intensity threshold mask
        intensity_thresh = 0.20; 
        intensity_mask   = IMGPR > intensity_thresh;
        IMGPR            = IMGPR .* intensity_mask;
    end

    if ~skipBinarization
        % Adaptive binarization and cleanup
        IMGPR            = imbinarize(IMGPR, 'adaptive', 'Sensitivity', 0.0);
        IMGPR            = imdilate(IMGPR, strel('disk', 2));
        IMGPR            = imerode(IMGPR, strel('disk', 1));
        IMGPR            = imfill(IMGPR, 'holes');
        F                = fft2(double(IMGPR)); % Compute 2D Fourier Transform
        IMGFFT           = abs(fftshift(F))';  % Shift zero frequency to center
    else
        F                = fft2(double(IMGPR)); % Compute 2D Fourier Transform
        IMGFFT           = abs(fftshift(F))';  % Shift zero frequency to center
    end
end

function [k_rho_vals, S_radial] = computeRadialSpectralDistribution(IMGFFT, kx, ky, thetamin, thetamax, num_bins)
    % IMGFFT    : 2D FFT image (fftshifted and cropped)
    % kx, ky    : 1D physical wavenumber axes [μm⁻¹] matching FFT size
    % thetamin  : Minimum angle (in radians)
    % thetamax  : Maximum angle (in radians)
    % num_bins  : Number of radial bins

    [KX, KY] = meshgrid(kx, ky);
    K_rho = sqrt(KX.^2 + KY.^2);
    Theta = atan2(KY, KX);

    if thetamin < thetamax
        angle_mask = (Theta >= thetamin) & (Theta <= thetamax);
    else
        angle_mask = (Theta >= thetamin) | (Theta <= thetamax);
    end

    power_spectrum = abs(IMGFFT).^2;

    r_min = min(K_rho(angle_mask));
    r_max = max(K_rho(angle_mask));
    r_edges = linspace(r_min, r_max, num_bins + 1);
    k_rho_vals = 0.5 * (r_edges(1:end-1) + r_edges(2:end));
    S_radial = zeros(1, num_bins);

    for i = 1:num_bins
        r_low = r_edges(i);
        r_high = r_edges(i + 1);
        radial_mask = (K_rho >= r_low) & (K_rho < r_high);
        full_mask = radial_mask & angle_mask;
        S_radial(i) = sum(power_spectrum(full_mask));
    end
end

function [theta_vals, S_theta] = computeAngularSpectralDistribution(IMGFFT, kx, ky, k_min, k_max, num_bins, threshold, sigma, windowSize)
    % Apply threshold to isolate strong peaks
    IMGFFT(IMGFFT < threshold) = 0;

    % Create wavenumber meshgrid
    [KX, KY] = meshgrid(kx, ky);
    Kmag     = sqrt(KX.^2 + KY.^2);       % radial wavenumber magnitude
    Theta    = atan2(KY, KX);             % range [-pi, pi]

    % Restrict to radial band in wavenumber space
    radial_mask = (Kmag >= k_min) & (Kmag <= k_max);

    % Initialize angular structure factor
    S_theta    = zeros(1, num_bins);
    theta_vals = linspace(0, pi, num_bins);  % only 0 to pi due to symmetry

    % Loop over angular bins
    for i = 1:num_bins
        angle_start = (i - 1) * pi / num_bins;
        angle_end   = i * pi / num_bins;
        angle_mask  = (Theta >= angle_start) & (Theta < angle_end);
        bin_mask    = radial_mask & angle_mask;
        fft_angle   = IMGFFT .* bin_mask;
        S_theta(i)  = sum(sum(abs(fft_angle).^2));
    end

    % Optional smoothing
    if exist('sigma', 'var') && ~isempty(sigma)
        % Gaussian smoothing
        half_width   = ceil(3 * sigma);
        x            = -half_width:half_width;
        gauss_kernel = exp(-x.^2 / (2 * sigma^2));
        gauss_kernel = gauss_kernel / sum(gauss_kernel);

        % Circular convolution
        S_theta = conv([S_theta(end - half_width + 1:end), S_theta, S_theta(1:half_width)], ...
                       gauss_kernel, 'same');
        S_theta = S_theta(half_width + 1:end - half_width);
    elseif exist('windowSize', 'var') && ~isempty(windowSize)
        % Moving average smoothing
        pad    = floor(windowSize / 2);
        kernel = ones(1, windowSize) / windowSize;
        S_theta = conv([S_theta(end - pad + 1:end), S_theta, S_theta(1:pad)], kernel, 'same');
        S_theta = S_theta(pad + 1:end - pad);
    end
end

function contrast = computeRadialSpectralContrast(k_rho_vals, S_k_smoothed, k_min, k_max)
%   Computes the ratio of the peak in S_k_smoothed within [k_min, k_max]
%   to the value at (or near) k = 0.

    % Ensure inputs are column vectors
    k_rho_vals     = k_rho_vals(:);
    S_k_smoothed   = S_k_smoothed(:);

    % Step 1: Find index of k ≈ 0
    [~, idx_k0] = min(abs(k_rho_vals));  % Closest to zero
    S_k0 = S_k_smoothed(idx_k0);

    % Step 2: Find indices in specified k-range
    in_range = (k_rho_vals >= k_min) & (k_rho_vals <= k_max);
    
    if ~any(in_range)
        warning('No values found in the specified k-range. Returning NaN.');
        contrast = NaN;
        return;
    end

    % Step 3: Find peak value in the specified k-range
    S_k_peak = max(S_k_smoothed(in_range));

    % Step 4: Compute contrast
    contrast = S_k_peak / S_k0;

end

function drawPSOverlays(kx, ky, r_min, r_max)
% drawFFTOverlays - Draw overlays on existing FFT plot:
%   - Radial lines every 30°
%   - Annular highlight with white (upper half) and gray (lower half) circles between r_min and r_max
%   - Horizontal white bands at ky=0 in annulus region
%   - Scale ticks and labels every 1 μm⁻¹ along each radial line
%
% Inputs:
%   kx, ky  - reciprocal space vectors (μm⁻¹)
%   r_min   - inner annulus radius offset index (integer)
%   r_max   - outer annulus radius offset index (integer)
%
% Example:
%   hold on;
%   drawFFTOverlays(kx, ky, 10, 30);
    
    hold on

    % === Overlay Radial Lines + Scales ===
    [kx_grid, ky_grid] = meshgrid(kx, ky);
    [~, kr_grid] = cart2pol(kx_grid, ky_grid);  % kr_grid in μm⁻¹

    max_kx = max(kx);
    max_ky = max(ky);

    for angle = 0 : pi/6 : pi
        x_line = [0, max_kx] * cos(angle);
        y_line = [0, max_ky] * sin(angle);

        % Plot radial lines
        plot(x_line, y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);
        plot(x_line, -y_line, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.2);

        % Draw scale ticks along positive radial line
        drawTicksAlongLine(0, 0, x_line(2), y_line(2));

        % Draw scale ticks along negative radial line (reflect y)
        drawTicksAlongLine(0, 0, x_line(2), -y_line(2));
    end

    % === Overlay Annular Highlight: White (r_min to r_max), Gray elsewhere ===
    theta_full = linspace(0, 2*pi, 500);

    center_x = ceil(size(kr_grid, 2) / 2);
    center_y = ceil(size(kr_grid, 1) / 2);

    k_min = kr_grid(center_y, center_x + r_min);
    k_max = kr_grid(center_y, center_x + r_max);

    % Upper half: white dashed circles
    x1_upper = k_min * cos(theta_full(theta_full <= pi));
    y1_upper = k_min * sin(theta_full(theta_full <= pi));
    x2_upper = k_max * cos(theta_full(theta_full <= pi));
    y2_upper = k_max * sin(theta_full(theta_full <= pi));
    plot(x1_upper, y1_upper, 'k--', 'LineWidth', 1.2);
    plot(x2_upper, y2_upper, 'k--', 'LineWidth', 1.2);

    % Lower half: gray dashed circles
    x1_lower = k_min * cos(theta_full(theta_full > pi));
    y1_lower = k_min * sin(theta_full(theta_full > pi));
    x2_lower = k_max * cos(theta_full(theta_full > pi));
    y2_lower = k_max * sin(theta_full(theta_full > pi));
    plot(x1_lower, y1_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);
    plot(x2_lower, y2_lower, '--', 'Color', [0.5 0.5 0.5], 'LineWidth', 1.0);

    % === Highlight horizontal band across k_y = 0 ===
    x_vals = kx;
    xW1 = x_vals((x_vals >= -k_max) & (x_vals < -k_min));
    xW2 = x_vals((x_vals > k_min) & (x_vals <= k_max));

    plot(xW1, zeros(size(xW1)), 'k--', 'LineWidth', 1.2);
    plot(xW2, zeros(size(xW2)), 'k--', 'LineWidth', 1.2);

    hold off


    % --- Nested helper function to draw ticks along a radial line ---
    function drawTicksAlongLine(x_start, y_start, x_end, y_end)
        % Tick parameters
        tick_spacing = 1;      % spacing between ticks in μm⁻¹
        tick_length = 0.05 * sqrt((x_end - x_start)^2 + (y_end - y_start)^2); % relative tick length
        line_color = [0.5 0.5 0.5];
        tick_color = [0.5 0.5 0.5];
        font_size = 8;

        % Vector along the line
        dx = x_end - x_start;
        dy = y_end - y_start;
        L = sqrt(dx^2 + dy^2);
        ux = dx / L;
        uy = dy / L;

        % Perpendicular vector for ticks
        perp_ux = -uy;
        perp_uy = ux;

        % Number of ticks (from 0 up to max length)
        n_ticks = floor(L / tick_spacing);

        for i = 1:n_ticks
            % Position of tick along the line
            xt = x_start + i * tick_spacing * ux;
            yt = y_start + i * tick_spacing * uy;

            % Tick endpoints
            xt1 = xt - 0.5 * tick_length * perp_ux;
            yt1 = yt - 0.5 * tick_length * perp_uy;
            xt2 = xt + 0.5 * tick_length * perp_ux;
            yt2 = yt + 0.5 * tick_length * perp_uy;

            % Draw tick
            plot([xt1 xt2], [yt1 yt2], '-', 'Color', tick_color, 'LineWidth', 1);

            % Label with distance (integer)
            text(xt, yt, sprintf('%d', i), ...
                'Color', tick_color, ...
                'FontSize', font_size, ...
                'HorizontalAlignment', 'center', ...
                'VerticalAlignment', 'bottom', ...
                'Rotation', atan2d(dy, dx));
        end
    end

end