Merge branch 'master' of https://git.physi.uni-heidelberg.de/karthik/Calculations

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
@@ -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])
@@ -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);
@@ -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,200 @@ 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;
+
+%% 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
 
 function ret = getBkgOffsetFromCorners(img, x_fraction, y_fraction)
@@ -439,25 +625,48 @@ function [optrefimages] = removefringesInImage(absimages, refimages, bgmask)
 end
 
 function [M] = ImagingResponseFunction(B)
-    x = -100:100;
+    x     = -128:127;
-    y = x;
+    y     = x;
     [X,Y] = meshgrid(x,y);
-    R = sqrt(X.^2+Y.^2);
+    R     = sqrt(X.^2+Y.^2);
-    PHI = atan2(X,Y)+pi;
+    PHI   = atan2(X,Y)+pi;
-    %fit parameters
+    
-    tau = B(1);
+    % Fit parameters
+    tau   = B(1);
     alpha = B(2);
-    S0 = B(3);
+    S0    = B(3);
-    phi = B(4);
+    phi   = B(4);
-    beta = B(5);
+    beta  = B(5);
     delta = B(6);
-    A = B(7);
+    A     = B(7);
-    C = B(8);
+    C     = B(8);
-    a = B(9);
+    a     = B(9);
-    U = heaviside(1-R/a).*exp(-R.^2/a^2/tau^2);
+    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);
+    p     = U.*exp(1i.*THETA);
-    M = A*abs((ifft2(real(exp(1i*delta).*fftshift(fft2(p)))))).^2 + C;
+    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)
@@ -483,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