From fbd808fc8bc19d460f98b25ad68ed0f8d7b02143 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 28 Feb 2025 17:30:35 +0100 Subject: [PATCH 1/2] Latest analysis routine. --- .../extractIRF.m | 136 +++++++++++++++--- 1 file changed, 116 insertions(+), 20 deletions(-) diff --git a/Imaging-Response-Function-Extractor/extractIRF.m b/Imaging-Response-Function-Extractor/extractIRF.m index 2950cbe..2506f4d 100644 --- a/Imaging-Response-Function-Extractor/extractIRF.m +++ b/Imaging-Response-Function-Extractor/extractIRF.m @@ -193,12 +193,12 @@ for k = 1 : length(od_imgs) drawnow; end -%% Compute the average 2D spectrum and do radial averaging to get the 1D spectrum +%% Compute the average 2D spectrum % Compute the average power spectrum. averagePowerSpectrum = mean(cat(3, density_noise_spectrum{:}), 3, 'double'); -% Plot the average power spectrum. +%% Plot the average power spectrum and the 1-D Radial Imaging Response Function figure(2) clf set(gcf,'Position',[100, 100, 1500, 700]) @@ -272,7 +272,7 @@ initialGuess = [2E6, 1E-6, 1E-6, 1E-6, 1E-6, 0.8]; % Set upper and lower bounds for the parameters (optional) lb = [-Inf, -Inf, -Inf, -Inf, -Inf, 0]; % Lower bounds -ub = [Inf, Inf, Inf, Inf, Inf, 1]; % Upper bounds +ub = [Inf, Inf, Inf, Inf, Inf, 1]; % Upper bounds % Perform the non-linear least squares fitting using lsqcurvefit options = optimoptions('lsqcurvefit', 'Display', 'iter'); % Display iterations during fitting @@ -285,14 +285,87 @@ fittedProfile = RadialImagingResponseFunction(C_fit, k_new, kmax); nexttile(2) plot(kvec, NormalisedProfile, 'o-', 'MarkerSize', 4, 'MarkerFaceColor', 'none', 'DisplayName', 'Radial (average) profile'); hold on -plot(k_new, fittedProfile, 'r-', 'DisplayName', 'Fitted Curve'); +plot(k_new, fittedProfile, 'r-', 'DisplayName', 'Fit'); set(gca, 'XScale', 'log'); % Setting x-axis to log scale xlabel('k_\rho (\mum^{-1})', 'Interpreter', 'tex', 'FontSize', 16) ylabel('Normalised amplitude', 'FontSize', 16) -title('Modulation Transfer Function', 'FontSize', 16); +title('Imaging Response Function', 'FontSize', 16); legend('FontSize', 16); grid on; +%% 2-D Imaging Response Function fit + +% Get the size of the matrix +[nRows, nCols] = size(averagePowerSpectrum); + +% Compute the center index +centerRow = round(nRows / 2); +centerCol = round(nCols / 2); + +% Remove the central DC component by setting it to zero +averagePowerSpectrum(centerRow, centerCol) = 0; + +% Define the objective function for fitting +function residuals = fitImagingResponse(params, averagePowerSpectrum) + % Calculate the ImagingResponseFunction with the given parameters + M = ImagingResponseFunction(params); + + % Compute the residuals as the difference between the matrices + residuals = M(:) - averagePowerSpectrum(:); % Flatten both matrices +end + +% Initial guess for parameters (from your example) +initialParams = [0.65, 0.1, 0.1, 0.1, 0.1, 0.1, 1E-6, 1E-8, 80]; + +% Set options for the optimizer (optional, helps with debugging) +options = optimoptions('lsqnonlin', 'Display', 'iter', 'MaxFunctionEvaluations', 5000); + +% Perform the fitting using lsqnonlin +optimalParams = lsqnonlin(@(params) fitImagingResponse(params, averagePowerSpectrum), initialParams, [], [], options); +AstigmatismCoefficient = optimalParams(2); +SphericalAberrationCoefficient = optimalParams(3); +DefocussingCoefficient = optimalParams(5); + +% Compute the best-fit M using the optimized parameters +bestFitM = ImagingResponseFunction(optimalParams); + +residuals = fitImagingResponse(optimalParams, averagePowerSpectrum); + +% Plot the fitted result + +figure(3) +clf +set(gcf,'Position',[100, 100, 1500, 500]) + +% Create tiled layout with 2 rows and 3 columns +t = tiledlayout(1, 3, 'TileSpacing', 'compact', 'Padding', 'compact'); + +nexttile(1); +imagesc(abs(10*log10(averagePowerSpectrum))) +axis equal tight; +colorbar +colormap(flip(jet)); +title('Average Density Noise Spectrum', 'FontSize', 16); +grid on; + +nexttile(2); +imagesc(abs(10*log10(bestFitM))) +axis equal tight; +colorbar +colormap(flip(jet)); +title('Fitted Imaging Response Function', 'FontSize', 16); +grid on; + +nexttile(3); +resmat = reshape(residuals, size(bestFitM)); +imagesc(abs(10*log10(resmat))) +axis equal tight; +colorbar +colormap(flip(jet)); +title('Fit residuals', 'FontSize', 16); +grid on; + + %% Helper Functions function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction) @@ -439,25 +512,48 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask) end function [M] = ImagingResponseFunction(B) - x = -100:100; - y = x; + x = -128:127; + y = x; [X,Y] = meshgrid(x,y); - R = sqrt(X.^2+Y.^2); - PHI = atan2(X,Y)+pi; - %fit parameters - tau = B(1); + R = sqrt(X.^2+Y.^2); + PHI = atan2(X,Y)+pi; + + % Fit parameters + tau = B(1); alpha = B(2); - S0 = B(3); - phi = B(4); - beta = B(5); + S0 = B(3); + phi = B(4); + beta = B(5); delta = B(6); - A = B(7); - C = B(8); - a = B(9); - U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2); + A = B(7); + C = B(8); + a = B(9); + U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2); THETA = S0*(R/a).^4 + alpha*(R/a).^2.*cos(2*PHI-2*phi) + beta*(R/a).^2; - p = U.*exp(1i.*THETA); - M = A*abs((ifft2(real(exp(1i*delta).*fftshift(fft2(p)))))).^2 + C; + p = U.*exp(1i.*THETA); + M = A*abs((ifft2(real(exp(1i*delta).*fftshift(fft2(p)))))).^2 + C; +end + +function [p] = exitPupilFunction(B) + x = -128:127; + y = x; + [X,Y] = meshgrid(x,y); + R = sqrt(X.^2+Y.^2); + PHI = atan2(X,Y)+pi; + + % Fit parameters + tau = B(1); + alpha = B(2); + S0 = B(3); + phi = B(4); + beta = B(5); + delta = B(6); + A = B(7); + C = B(8); + a = B(9); + U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2); + THETA = S0*(R/a).^4 + alpha*(R/a).^2.*cos(2*PHI-2*phi) + beta*(R/a).^2; + p = U.*exp(1i.*THETA); end function [R, Zr] = getRadialProfile(data, radialStep) From bfe828d1a9aedf6a551b3a2b41a8689477a4cfc1 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 28 Feb 2025 20:40:43 +0100 Subject: [PATCH 2/2] 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