Calculations/Dipolar-Gas-Simulator/+Scripts/extractLatticeProperties.m

function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y] = extractLatticeProperties(I, x, y)
    % Detects lattice geometry, extracts periodic spacings, and reconstructs real-space lattice vectors.
    % Handles arbitrary lattice geometries
    %
    % Inputs:
    %   I - Grayscale image with periodic structures.
    %   x - X-coordinates in micrometers.
    %   y - Y-coordinates in micrometers.
    %
    % Outputs:
    %   lattice_type - Identified lattice type (Square, Rectangular, Hexagonal, Oblique, or Unknown)
    %   dx_um - Spacing along x-axis in micrometers.
    %   dy_um - Spacing along y-axis in micrometers.
    %   real_lattice_vectors - [2x2] matrix of lattice vectors in micrometers.
    %   angle_in_reciprocal_space - Angle of the reciprocal lattice primitive vectors in degrees.

    % Compute 2D Fourier Transform
    F            = fft2(double(I));
    fftMagnitude = abs(fftshift(F))';  % Shift zero frequency to center
    
    % 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 
    
    % Find peaks in Fourier domain
    peak_mask = imregionalmax(fftMagnitude); % Detect local maxima
    threshold = max(fftMagnitude(:)) * 0.1; % Adaptive threshold
    peak_mask = peak_mask & (fftMagnitude > threshold);

    % Extract peak positions
    [rows, cols] = find(peak_mask);
    
    % Convert indices to frequency values
    freq_x = vx(cols);
    freq_y = vy(rows);
    distances = sqrt(freq_x.^2 + freq_y.^2);

    % Remove DC component (center peak)
    valid_idx = distances > 1e-3;
    freq_x = freq_x(valid_idx);
    freq_y = freq_y(valid_idx);
    distances = distances(valid_idx);

    % Sort by frequency magnitude
    [~, idx] = sort(distances);
    freq_x = freq_x(idx);
    freq_y = freq_y(idx);

    % Select first two unique peaks as lattice vectors
    unique_vectors = unique([freq_x', freq_y'], 'rows', 'stable');

    if size(unique_vectors, 1) < 4
        warning('Not enough detected peaks for lattice vector determination.');
        lattice_type   = 'Undetermined';
        dx_um          = NaN;
        dy_um          = NaN;
        freq_x         = NaN;
        freq_y         = NaN;
        angle_in_reciprocal_space = NaN;
    else

        % Select two shortest independent lattice vectors - Reciprocal lattice vectors
        % Reciprocal lattice vectors
        G1 = unique_vectors(1, :);  % Reciprocal lattice vector 1
        G2 = unique_vectors(3, :);  % Reciprocal lattice vector 2
        
        reciprocal_lattice_vectors = abs(vertcat(G1, G2));
    
        % Calculate the angle between the reciprocal lattice vectors using dot product
        dotProduct = dot(G1, G2);  % Dot product of G1 and G2
        magnitudeG1 = norm(G1);  % Magnitude of G1
        magnitudeG2 = norm(G2);  % Magnitude of G2
        
        % Angle between the reciprocal lattice vectors (in degrees)
        angle_in_reciprocal_space = rad2deg(acos(dotProduct / (magnitudeG1 * magnitudeG2)));
        
        % Convert to real-space lattice vectors
        
        % Compute the determinant of the reciprocal lattice vectors
        detG = G1(1) * G2(2) - G1(2) * G2(1); % Determinant of the matrix formed by G1 and G2
    
        if detG == 0
            % Handle the case of reciprocal lattice vector matrix being singular
            
            v1 = reciprocal_lattice_vectors(1, :);
            v2 = reciprocal_lattice_vectors(2, :);
        
            % Check if either vector is the zero vector
            if all(v1 == 0) && all(v2 == 0)
                % If both vectors are zero, return the same zero matrix
                real_lattice_vectors = [0, 0; 0, 0];
            end
            
            % Compute the norms of the vectors
            norm_v1 = norm(v1);
            norm_v2 = norm(v2);
            
            % Find the vector with the smallest norm
            if norm_v2 < norm_v1
                % If v2 has the smaller norm, set v1 to zero and keep v2
                reciprocal_lattice_vectors = [0, 0; v2];
                if v2(1) == 0
                    real_lattice_vectors = [0, 0; 0, 1/v2(2)];
                elseif v2(2) == 0
                    real_lattice_vectors = [0, 0; 1/v2(1), 0];
                else
                    real_lattice_vectors = [1/v2(1), 1/v2(2); 0, 0];
                end
            else
                % If v1 has the smaller norm, set v2 to zero and keep v1
                reciprocal_lattice_vectors = [v1; 0, 0];
                if v1(1) == 0
                    real_lattice_vectors = [0, 1/v1(2); 0, 0];
                elseif v1(2) == 0
                    real_lattice_vectors = [1/v1(1), 0; 0, 0];
                else
                    real_lattice_vectors = [1/v1(1), 1/v1(2); 0, 0];
                end
            end

        else
            % If reciprocal vectors are not collinear, compute the 2D real-space lattice
            real_lattice_vectors = inv(reciprocal_lattice_vectors'); % If vx, vy are used, no need to multiply by 2*pi (crystallographic convention); If kx, ky are used, need to multiply by 2*pi (physics convention)
        end
        
        % Ensure correct orientation
        real_lattice_vectors(isinf(real_lattice_vectors) | isnan(real_lattice_vectors)) = 0;
        
        % Separate the components of real space vectors
        real_x          = real_lattice_vectors(:, 1);
        real_y          = real_lattice_vectors(:, 2);
        
        % Compute projections
        proj_real_x          = real_lattice_vectors * [1; 0];  % Projections onto x-axis
        [smaller_value, idx] = min(proj_real_x);
        proj_real_x(idx)     = [];  % Remove the element
        proj_real_y          = real_lattice_vectors * [0; 1];  % Projections onto y-axis
        [smaller_value, idx] = min(proj_real_y);
        proj_real_y(idx)     = [];  % Remove the element

        % Separate the components of reciprocal space vectors
        reciprocal_x    = reciprocal_lattice_vectors(:, 1);
        reciprocal_y    = reciprocal_lattice_vectors(:, 2);

        % Compute projections
        proj_reciprocal_x     = reciprocal_lattice_vectors * [1; 0];  % Projections onto x-axis
        [smaller_value, idx] = min(proj_reciprocal_x);
        proj_reciprocal_x(idx)     = [];  % Remove the element
        proj_reciprocal_y     = reciprocal_lattice_vectors * [0; 1];  % Projections onto y-axis
        [smaller_value, idx] = min(proj_reciprocal_y);
        proj_reciprocal_y(idx)     = [];  % Remove the element
        
        %{
        % Create a figure for side-by-side subplots
        figure('Position', [100, 100, 1600, 800]);
        clf
        t = tiledlayout(1, 2, 'TileSpacing', 'compact', 'Padding', 'compact'); % 2x3 grid
        
        % Plot real space vectors in the first subplot (2D case)
        nexttile;  % 1 row, 2 columns, first plot
        quiver(0, 0, real_x(1) * 1E6, real_y(1) * 1E6, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);
        hold on;
        quiver(0, 0, real_x(2) * 1E6, real_y(2) * 1E6, 'b', 'LineWidth', 2, 'MaxHeadSize', 0.5);
        
        % Plot projections of the real space vectors
        quiver(0, 0, proj_real_x * 1E6, 0, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on X-axis (horizontal dotted line)
        quiver(0, 0, 0, proj_real_y * 1E6, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on Y-axis (vertical dotted line)
        
        hold off;
        xlabel('X');
        ylabel('Y');
        title('Real Space Primitive Vectors', 'FontSize', 14);
        grid on;
        axis equal;
        
        % Plot reciprocal space vectors in the second subplot (2D case)
        nexttile;  % 1 row, 2 columns, second plot
        quiver(0, 0, reciprocal_x(1) * 1E-6, reciprocal_y(1) * 1E-6, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);
        hold on;
        quiver(0, 0, reciprocal_x(2) * 1E-6, reciprocal_y(2) * 1E-6, 'b', 'LineWidth', 2, 'MaxHeadSize', 0.5);

        % Plot projections of the real space vectors
        quiver(0, 0, proj_reciprocal_x * 1E-6, 0, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on X-axis (horizontal dotted line)
        quiver(0, 0, 0, proj_reciprocal_y * 1E-6, 'k--', 'LineWidth', 2, 'MaxHeadSize', 0.5); % Projection on Y-axis (vertical dotted line)
        
        hold off;
        xlabel('X');
        ylabel('Y');
        title('Reciprocal Space Primitive Vectors', 'FontSize', 14);
        grid on;
        axis equal;
        %}

        % Compute lattice spacings
        dx_um               = proj_real_x * 1E6; 
        dy_um               = proj_real_y * 1E6;  
        
        % Classify lattice type based on angle symmetry
        angle_in_real_space = abs(atan2d(real_lattice_vectors(2,2), real_lattice_vectors(2,1)) - ...
                                atan2d(real_lattice_vectors(1,2), real_lattice_vectors(1,1)));
        angle_in_real_space_mod180 = mod(angle_in_real_space, 180);

        if all(real_lattice_vectors(1, :) == 0) || all(real_lattice_vectors(2, :) == 0)
            lattice_type              = 'Stripe';
            angle_in_real_space       = NaN;
            angle_in_reciprocal_space = NaN;
        else
            if abs(angle_in_real_space_mod180 - 60) < 5 || abs(angle_in_real_space_mod180 - 120) < 5
                if abs(dx_um - dy_um) < 1E-6
                    lattice_type = 'Triangular';
                else
                    lattice_type = 'Sheared Triangular';
                end
            elseif abs(angle_in_real_space_mod180 - 90) < 5
                lattice_type = 'Square';
            elseif abs(angle_in_real_space_mod180 - 90) > 5 && abs(angle_in_real_space_mod180 - 120) > 5
                lattice_type = 'Sheared';
            else
                lattice_type = 'Undetermined';
            end
        end

        % Display results
        fprintf('Detected Lattice Type: %s\n', lattice_type);
        fprintf('Estimated Spacing (dx): %.2f µm\n', dx_um);
        fprintf('Estimated Spacing (dy): %.2f µm\n', dy_um);
        fprintf('Angle between real space lattice primitve vectors: %.2f°\n', angle_in_real_space);
        fprintf('Angle between reciprocal lattice primitve vectors: %.2f°\n', angle_in_reciprocal_space);
        fprintf('Real Space Primitive Vectors:\n');
        disp(real_lattice_vectors);
    end
end