555 lines
20 KiB
Matlab
555 lines
20 KiB
Matlab
%% 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
|