From fbd808fc8bc19d460f98b25ad68ed0f8d7b02143 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Fri, 28 Feb 2025 17:30:35 +0100 Subject: [PATCH] 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)