From bfe828d1a9aedf6a551b3a2b41a8689477a4cfc1 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 28 Feb 2025 20:40:43 +0100 Subject: [PATCH] Added functionality to extract Zernike Coefficients from decomposing the Imaging Response Function in terms of the polynomials. --- .../extractIRF.m | 149 +++++++++++++++++- 1 file changed, 144 insertions(+), 5 deletions(-) diff --git a/Imaging-Response-Function-Extractor/extractIRF.m b/Imaging-Response-Function-Extractor/extractIRF.m index 2506f4d..f697838 100644 --- a/Imaging-Response-Function-Extractor/extractIRF.m +++ b/Imaging-Response-Function-Extractor/extractIRF.m @@ -4,14 +4,15 @@ groupList = ["/images/MOT_3D_Camera/in_situ_absorption", "/images/ODT_1_Axis folderPath = "C:/Users/Karthik/Documents/GitRepositories/Calculations/Imaging-Response-Function-Extractor/"; -run = '0096'; +run = '0072'; folderPath = strcat(folderPath, run); cam = 5; angle = 0; -center = [1137, 2023]; +% center = [1137, 2023]; +center = [1141, 2049]; span = [255, 255]; fraction = [0.1, 0.1]; @@ -109,7 +110,6 @@ R = 32; % aperture radius A = (X.^2 + Y.^2 <= R^2); % circular aperture of radius R mask(A) = 1; % set mask elements inside aperture to 1 - % Calculate Power Spectrum and plot figure(1) clf @@ -198,7 +198,7 @@ end % Compute the average power spectrum. averagePowerSpectrum = mean(cat(3, density_noise_spectrum{:}), 3, 'double'); -%% Plot the average power spectrum and the 1-D Radial Imaging Response Function +% Plot the average power spectrum and the 1-D Radial Imaging Response Function figure(2) clf set(gcf,'Position',[100, 100, 1500, 700]) @@ -223,7 +223,7 @@ yc = centers(:,2); [yDim, xDim] = size(averagePowerSpectrum); [xx,yy] = meshgrid(1:yDim,1:xDim); mask = false(xDim,yDim); -for ii = 1:length(centers) +for ii = 1:size(centers, 1) mask = mask | hypot(xx - xc(ii), yy - yc(ii)) <= radius; end mask = not(mask); @@ -365,6 +365,119 @@ colormap(flip(jet)); title('Fit residuals', 'FontSize', 16); grid on; +%% 3-D Plot of Density Noise Spectrum + +figure(4) +clf +set(gcf,'Position',[100, 100, 1500, 700]) +surf(10 * log10(averagePowerSpectrum)); +shading interp; % Creates a 3D surface plot +xlabel('X-axis', 'FontSize', 16); +ylabel('Y-axis', 'FontSize', 16); +zlabel('Rescaled amplitude (a.u.)', 'FontSize', 16); +title('Imaging Response Function', 'FontSize', 16); +colorbar; % Shows the color scale + +% Set axis limits based on the size of the averagePowerSpectrum +[xSize, ySize] = size(averagePowerSpectrum); +xlim([1, xSize]); +ylim([1, ySize]); +zlim([min(10 * log10(averagePowerSpectrum(:))), max(10 * log10(averagePowerSpectrum(:)))]); % Optional for Z-axis limits + + +%% Decompose in Zernike Polynomial basis + +N = size(averagePowerSpectrum, 1); +[X, Y] = meshgrid(linspace(-1, 1, N)); + +max_n = 6; % Adjust based on your needs +basis = []; +orders = []; + +for n = 0:max_n + for m = -n:2:n + % Generate Zernike polynomial for (n, m) + Z = zernike_polynomial(n, m, X, Y); + % Flatten and store valid points + basis = [basis, Z(mask)]; + orders = [orders; [n, m]]; + end +end + +data = 10 * log10(averagePowerSpectrum); +valid_data = data(mask); + +% Solve Ax = b (A = basis matrix, b = data) +coeffs = basis \ valid_data(:); + +% Reconstruct the surface using the coefficients +reconstructed = basis * coeffs; +reconstructed_surface = zeros(size(X)); +reconstructed_surface(mask) = reconstructed; + +figure(5) +clf +set(gcf,'Position',[100, 100, 1500, 700]) + +% Create tiled layout with 2 rows and 3 columns +t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); + +nexttile(1); +imagesc(data); title('Imaging Response Function', 'FontSize', 16); +axis square; +colorbar +colormap(jet); +grid on; + +nexttile(2); +imagesc(reconstructed_surface); title('Reconstructed with Zernike', 'FontSize', 16); +axis square; +colorbar +colormap(jet); +grid on; + +nexttile(3); +imagesc(data - reconstructed_surface); title('Residuals', 'FontSize', 16); +axis square; +colorbar +colormap(jet); +grid on; + +disp('Zernike Coefficients:'); +disp('---------------------'); +for i = 1:length(coeffs) + fprintf('Order (n=%d, m=%d): Coefficient = %.4f\n', orders(i,1), orders(i,2), coeffs(i)); +end + +% Plot Zernike Coeffecients + +% Find the index of the (n=0, m=0) term +idx_remove = find(orders(:,1) == 0 & orders(:,2) == 0); + +% Remove the Z_0^0 term from coefficients and orders +coeffs_filtered = coeffs; +coeffs_filtered(idx_remove) = []; +orders_filtered = orders; +orders_filtered(idx_remove, :) = []; + +% Generate labels for filtered modes (n, m) +labels_filtered = cell(length(coeffs_filtered), 1); +for i = 1:length(coeffs_filtered) + labels_filtered{i} = sprintf('(%d, %d)', orders_filtered(i,1), orders_filtered(i,2)); +end + +figure(6) +clf +set(gcf,'Position',[100, 100, 1500, 700]) +bar(coeffs_filtered, 'FaceColor', [0.2, 0.6, 0.8]); % Customize bar color +ylim([-1.0, 1.0]) +title('Zernike Coefficients', 'FontSize', 16); +xlabel('Zernike Mode (n, m)', 'FontSize', 16); +ylabel('Coefficient Value', 'FontSize', 16); +xticks(1:length(coeffs_filtered)); % Set x-ticks for all coefficients +xticklabels(labels_filtered); % Assign (n, m) labels +xtickangle(45); % Rotate labels for readability +grid on; %% Helper Functions @@ -579,3 +692,29 @@ function [RadialResponseFunc] = RadialImagingResponseFunction(C, k, kmax) end RadialResponseFunc = C(6) * 1/2 * A .* RadialResponseFunc; end + +function R = zernike_radial(n, m, rho) + % Compute radial part of Zernike polynomial for radial order n and azimuthal order m + R = zeros(size(rho)); + for k = 0:(n - abs(m))/2 + coeff = (-1)^k * factorial(n - k) / ... + (factorial(k) * factorial((n + abs(m))/2 - k) * factorial((n - abs(m))/2 - k)); + R = R + coeff * rho.^(n - 2*k); + end +end + +function Z = zernike_polynomial(n, m, X, Y) + % Convert Cartesian coordinates to polar coordinates + [theta, rho] = cart2pol(X, Y); + rho(rho > 1) = 0; % Restrict to unit disk + + % Compute radial component + R = zernike_radial(n, m, rho); + + % Compute azimuthal component (sine/cosine terms) + if m >= 0 + Z = R .* cos(m * theta); + else + Z = R .* sin(-m * theta); + end +end \ No newline at end of file