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

function [kx, ky, fftMagnitude, lattice_type, dx_um, dy_um, freq_x, freq_y, rotation_angle] = extractLatticeProperties(I, x, y)
    % Detects lattice geometry, extracts periodic spacings, and reconstructs real-space lattice vectors.
    % Handles arbitrary lattice geometries with rotation correction.
    %
    % 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.
    %   rotation_angle - Angle of rotation of the lattice in degrees.

    % Compute 2D Fourier Transform
    F            = fft2(double(I));
    fftMagnitude = abs(fftshift(F))';  % Shift zero frequency to center
    
    % Calculate wavenumber increment (frequency axes)
    Nx = length(x);
    Ny = length(y);
    dx = mean(diff(x)); 
    dy = mean(diff(y)); 
    dkx = 1 / (Nx * dx);  % Increment in the X direction (in micrometers^-1)
    dky = 1 / (Ny * dy);  % Increment in the Y direction (in micrometers^-1)

    % Create the wavenumber axes
    kx = (-Nx/2:Nx/2-1) * dkx;  % Wavenumber axis in X (micrometers^-1)
    ky = (-Ny/2:Ny/2-1) * dky;  % Wavenumber axis in Y (micrometers^-1)

    % 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 = kx(cols);
    freq_y = ky(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;
        rotation_angle = 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 = 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)
        rotation_angle = 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 collinear reciprocal lattice vectors
            disp('Reciprocal lattice vectors are collinear.');
            
            % If both G1, G2 or just G1 are along x-direction
            if G1(1) ~= 0 && G1(2) == 0 
                a1 = [1 / norm(G1), 0];  % Real-space lattice vector along x
                a2 = [0, 0];             % No periodicity in the y-direction (1D lattice)
            % If G2 is along x-direction or if G1 is zero.
            elseif G2(1) ~= 0 && G2(2) == 0
                a1 = [1 / norm(G2), 0];             % Real-space lattice vector along x
                a2 = [0, 0];  % No periodicity in the y-direction
            % If G1 and G2 are both in the y-direction (i.e., [0, non-zero])
            elseif G1(1) == 0 && G2(1) == 0
                a1 = [0, 1 / norm(G1)];  % Real-space lattice vector along y
                a2 = [0, 0];  % No periodicity in the x-direction
            % If both vectors are zero (no periodicity)
            else
                a1 = [0, 0];
                a2 = [0, 0];
            end
        else
            % If reciprocal vectors are not collinear, compute the 2D real-space lattice
            a1 = [G2(2), -G1(2)] / detG;
            a2 = [-G2(1), G1(1)] / detG;
        end
        
        real_lattice_vectors = vertcat(a1, a2);
    
        % Ensure correct orientation
        real_lattice_vectors(isinf(real_lattice_vectors) | isnan(real_lattice_vectors)) = 0;
        
        % Compute lattice spacings
        dx_um = (1/norm(G2)) * 1E6; 
        dy_um = (1/norm(G1)) * 1E6;  
    
        % Correct for rotation by applying inverse rotation matrix
        theta_rad = deg2rad(-rotation_angle); % Convert to radians
        R = [cos(theta_rad), -sin(theta_rad); sin(theta_rad), cos(theta_rad)]; % Rotation matrix
        real_lattice_vectors = (R * real_lattice_vectors')'; % Apply rotation correction
    
        % Classify lattice type based on angle symmetry
        angle_between_vectors = atan2d(real_lattice_vectors(2,2), real_lattice_vectors(2,1)) - ...
                                atan2d(real_lattice_vectors(1,2), real_lattice_vectors(1,1));
        angle_between_vectors = mod(angle_between_vectors, 180);
    
        if abs(angle_between_vectors - 60) < 5 || abs(angle_between_vectors - 120) < 5
            lattice_type = 'Hexagonal';
        elseif abs(angle_between_vectors - 90) < 5
            lattice_type = 'Square';
        elseif abs(angle_between_vectors - 90) > 5 && abs(angle_between_vectors - 120) > 5
            lattice_type = 'Oblique';
        else
            lattice_type = 'Unknown';
        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('Rotation Angle: %.2f°\n', rotation_angle);
        % fprintf('Lattice Vectors:\n');
        % disp(real_lattice_vectors);
    end
end