%% Extract Images

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

TitleString  = 'Change across transition: Droplets To Stripes';

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


 TitleString  = 'Change across transition: 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
    x = Transf.x * Params.l0 * 1e6;
    y = Transf.y * Params.l0 * 1e6;
    z = Transf.z * Params.l0 * 1e6;
    dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);

    % 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)
    
    % Calculate maximum frequencies
    % kx_max       = pi / dx;
    % ky_max       = pi / dy;
    
    % Generate reciprocal axes
    % kx           = linspace(-kx_max, kx_max * (Nx-2)/Nx, Nx);
    % ky           = linspace(-ky_max, ky_max * (Ny-2)/Ny, Ny);
    
    % Create the Wavenumber axes
    kx           = 2*pi*vx;  % Wavenumber axis in X
    ky           = 2*pi*vy;  % Wavenumber axis in Y

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

    imgs{k} = nxy;
end

%% Analyze Images


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

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

% Run Fourier analysis over images

fft_imgs              = cell(1, nimgs);
spectral_contrast     = zeros(1, nimgs);
spectral_weight       = zeros(1, nimgs);
g2_all                = cell(1, nimgs);
theta_values_all      = cell(1, nimgs);

N_bins                = 180;
Threshold             = 25;
Sigma                 = 2;

if makeMovie
    % Create VideoWriter object for movie
    videoFile = VideoWriter(movieFileName, '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

% Display the cropped image
for k = 1:nimgs
    IMG                    = imgs{k};
    [IMGFFT, IMGPR]        = computeFourierTransform(IMG, skipPreprocessing, skipMasking, skipIntensityThresholding, skipBinarization);
    [theta_vals, S_theta]  = computeNormalizedAngularSpectralDistribution(IMGFFT, 10, 35, N_bins, Threshold, Sigma);

    g2 = zeros(1, N_bins);  % Preallocate
    
    for dtheta = 0:N_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 - 1;
    end
    
    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, 'jet');
    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
    
    % Plot the power spectrum
    ax2 = nexttile;
    [rows, cols] = size(IMGFFT);
    zoom_size    = 50;             % Zoomed-in region around center
    mid_x        = floor(cols/2);
    mid_y        = floor(rows/2);
    fft_imgs{k}  = IMGFFT(mid_y-zoom_size:mid_y+zoom_size, mid_x-zoom_size:mid_x+zoom_size);
    imagesc(log(1 + abs(fft_imgs{k}).^2));
    % 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
    axis square;
    hcb = colorbar;
    colormap(ax2, 'jet');
    set(gca, 'FontSize', 14);   % For tick labels only
    set(gca,'YDir','normal')
    hXLabel = xlabel('k_x', 'Interpreter', 'tex');
    hYLabel = ylabel('k_y', 'Interpreter', 'tex');
    hTitle  = title('Power Spectrum - S(k_x,k_y)', 'Interpreter', 'tex');
    set([hXLabel, hYLabel, hText], 'FontName', font)
    set([hXLabel, hYLabel], 'FontSize', 14)
    set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title

    % Plot the angular distribution
    nexttile
    spectral_contrast(k)  = computeSpectralContrast(fft_imgs{k}, 10, 25, Threshold);
    [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(fft_imgs{k}, 10, 25, N_bins, Threshold, Sigma);
    spectral_weight(k)    = trapz(theta_vals, S_theta);
    plot(theta_vals/pi, S_theta,'Linewidth',2);
    axis square;
    set(gca, 'FontSize', 14);   % For tick labels only
    hXLabel = xlabel('\theta/\pi [rad]', 'Interpreter', 'tex');
    hYLabel = ylabel('Normalized magnitude (a.u.)', 'Interpreter', 'tex');
    hTitle = title('Angular Spectral Distribution - S(\theta)', 'Interpreter', 'tex');
    set([hXLabel, hYLabel, hText], 'FontName', font)
    set([hXLabel, hYLabel], 'FontSize', 14)
    set(hTitle, 'FontName', font, 'FontSize', 16, 'FontWeight', 'bold');  % Set font and size for title
    grid on

    nexttile
    plot(theta_vals/pi, g2, 'o-', 'LineWidth', 1.2, 'MarkerSize', 5);
    set(gca, 'FontSize', 14);
    ylim([-1.5 3.0]);  % Set y-axis limits here
    hXLabel = xlabel('$\delta\theta / \pi$', 'Interpreter', 'latex');
    hYLabel = ylabel('$g^{(2)}(\delta\theta)$', 'Interpreter', 'latex');
    hTitle  = title('Autocorrelation', '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 makeMovie
        frame = getframe(gcf);
        writeVideo(videoFile, frame);
    else
        pause(0.5);  % Only pause when not recording
    end

end

if makeMovie
    close(videoFile);
    disp(['Movie saved to ', movieFileName]);
end

%% Track across the transition

figure(2);
set(gcf,'Position',[100 100 950 750])
plot(alphas, spectral_contrast, 'o--', ...
    'LineWidth', 1.5, 'MarkerSize', 6);
set(gca, 'FontSize', 14);   % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex');
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
hTitle = title(TitleString, '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

figure(3);
set(gcf,'Position',[100 100 950 750])
plot(alphas, spectral_weight, 'o--', ...
    'LineWidth', 1.5, 'MarkerSize', 6);
set(gca, 'FontSize', 14);   % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
% hXLabel = xlabel('B_z (G)', 'Interpreter', 'tex');
hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
hTitle = title(TitleString, '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

save(datafileName, 'alphas', 'spectral_contrast', 'spectral_weight');

figure(4);
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([-1.5 3.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;

%% Track across the transition

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

format long

font                       = 'Bahnschrift';

% Load data
Data                       = load('./DropletsToStripes.mat', 'alphas', 'spectral_contrast', 'spectral_weight');
dts_alphas                 = Data.alphas;
dts_sc                     = Data.spectral_contrast;
dts_sw                     = Data.spectral_weight;

Data                       = load('./StripesToDroplets.mat', 'alphas', 'spectral_contrast', 'spectral_weight');
std_alphas                 = Data.alphas;
std_sc                     = Data.spectral_contrast;
std_sw                     = Data.spectral_weight;

% Normalize dts data
dts_min                    = min(dts_sw);
dts_max                    = max(dts_sw);
dts_range                  = dts_max - dts_min;
dts_sf_norm                = (dts_sw - dts_min) / dts_range;

% Normalize std data
std_min                    = min(std_sw);
std_max                    = max(std_sw);
std_range                  = std_max - std_min;
std_sf_norm                = (std_sw - std_min) / std_range;

figure(5);
set(gcf,'Position',[100 100 950 750])
plot(dts_alphas, dts_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes');
hold on
plot(std_alphas, std_sc, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets');
set(gca, 'FontSize', 14);   % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Spectral Contrast', 'Interpreter', 'tex');
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

figure(6);
set(gcf,'Position',[100 100 950 750])
plot(dts_alphas, dts_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Droplets to Stripes');
hold on
plot(std_alphas, std_sw, 'o--', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName' , 'Stripes to Droplets');
set(gca, 'FontSize', 14);   % For tick labels only
hXLabel = xlabel('\alpha (degrees)', 'Interpreter', 'tex');
hYLabel = ylabel('Spectral Weight', 'Interpreter', 'tex');
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 [theta_vals, S_theta] = computeNormalizedAngularSpectralDistribution(IMGFFT, r_min, r_max, num_bins, threshold, sigma)
    % Apply threshold to isolate strong peaks
    IMGFFT(IMGFFT < threshold) = 0;

    % Prepare polar coordinates
    [ny, nx]  = size(IMGFFT);
    [X, Y]    = meshgrid(1:nx, 1:ny);
    cx        = ceil(nx/2); 
    cy        = ceil(ny/2);
    R         = sqrt((X - cx).^2 + (Y - cy).^2);
    Theta     = atan2(Y - cy, X - cx); % range [-pi, pi]

    % Choose radial band
    radial_mask      = (R >= r_min) & (R <= r_max);

    % Initialize the angular structure factor array
    S_theta = zeros(1, num_bins); % Pre-allocate for 180 angle bins
    % Define the angle values for the x-axis
    theta_vals = linspace(0, pi, num_bins);

    % Loop through each angle bin
    for i = 1:num_bins
        angle_start = (i-1) * pi / num_bins;
        angle_end   = i * pi / num_bins;
        
        % Define a mask for the given angle range
        angle_mask  = (Theta >= angle_start & Theta < angle_end);
        
        bin_mask    = radial_mask & angle_mask;

        % Extract the Fourier components for the given angle
        fft_angle   = IMGFFT .* bin_mask;
        
        % Integrate the Fourier components over the radius at the angle
        S_theta(i)  = sum(sum(abs(fft_angle).^2)); % sum of squared magnitudes
    end
    
    % Create a 1D Gaussian kernel
    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); % normalize
    
    % Apply convolution (circular padding to preserve periodicity)
    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); % crop back to original size
    
    % Normalize to 1
    S_theta = S_theta / max(S_theta);
end

function contrast = computeSpectralContrast(IMGFFT, r_min, r_max, threshold)
    % Apply threshold to isolate strong peaks
    IMGFFT(IMGFFT < threshold) = 0;

    % Prepare polar coordinates
    [ny, nx] = size(IMGFFT);
    [X, Y]   = meshgrid(1:nx, 1:ny);
    cx       = ceil(nx/2); 
    cy       = ceil(ny/2);
    R        = sqrt((X - cx).^2 + (Y - cy).^2);

    % Ring region (annulus) mask
    ring_mask = (R >= r_min) & (R <= r_max);

    % Squared magnitude in the ring
    ring_power = abs(IMGFFT).^2 .* ring_mask;

    % Maximum power in the ring
    ring_max = max(ring_power(:));

    % Power at the DC component
    dc_power = abs(IMGFFT(cy, cx))^2;

    % Avoid division by zero
    if dc_power == 0
        contrast = Inf;  % or NaN or 0, depending on how you want to handle this
    else
        contrast = ring_max / dc_power;
    end
end