problem with defocus aberration simulation code

Hey everyone,

I've been working on this code based on Chapter 5 of the book 'Numerical Simulation of Optical Wave Propagation with Examples in MATLAB'. It's about imaging systems. I wanted to represent defocus aberration, and I came up with the following code. I didn't use the 'zernike.m' and 'zernike.index' files like in the book. I kept it simple and just used the expression for the defocus mode z(0,2). < 2*r^.2-1

Here's the code:

% example_incoh_img.m


% Parameter Definitions
    N = 2000;              % number of grid points per side
    L = 0.3;               % total size of the grid [m]
    D = 0.004;             % diameter of pupil [m]
    delta = L / N;         % grid spacing [m]
    wvl = 5.5e-7;          % optical wavelength [m]
    z = 0.7;               % image distance [m]

% pupil-plane coordinates
    [x, y] = meshgrid((-N/2 : N/2-1) * delta);
    [theta, r] = cart2pol(x, y);

% wavefront aberration
    W = 0.05 * sqrt(3)*(2*r.^2-1);

% complex pupil function
    P = circ(x, y, D).* exp(-1i * 2*pi * W);

% amplitude spread function
    h = ft2(P, delta);
    U = wvl * z / (N*delta);

%PSF
    PSF = abs(h).^2;

% image-plane coordinates
    [u, v] = meshgrid((-N/2 : N/2-1) * U);

%    object    (same coordinates as h)
    obj = (rect((u-1.4e-4)/5e-5) + rect(u/5e-5)+ rect((u+1.4e-4)/5e-5)) .* rect(v/2e-4);

% convolve the object with the PSF to simulate imaging
    img = myconv2(abs(obj).^2, PSF, 1);

% Display the original object image.
    figure,imshow(obj);
    title('Original image');

% Resize the image for visualizing the PSF
    psf_img = abs(h).^2;
    psf_img = psf_img / max(psf_img(:)); % Normalize the PSF for better visualization

% Visualize the PSF
    figure;
    imshow(psf_img, []);
    title('Point Spread Function (PSF)');
    colormap jet;
    colorbar;

% Visualize the simulated result image
    figure;
    imshow(img, []);
    title('Simulated Image');

>> Code used for performing convolution

    function C = myconv2(A, B, delta)
% function C = myconv2(A, B, delta)
    N = size(A, 1);
    C = ift2(ft2(A, delta) .* ft2(B, delta), 1/(N*delta));

>> Code used for Fourier transform

    function G = ft2(g, delta)
% function G = ft2(g, delta)

    G = fftshift(fft2(fftshift(g))) * delta^2;

>> Code used for inverse Fourier transform

    function g = ift2(G, delta_f)
% function g = ift2(G, delta_f)
    N = size(G, 1);
    g = ifftshift(ifft2(ifftshift(G))) * (N * delta_f)^2;

>> Code used to create the amplitude function

    function z = circ(x, y, D)
% function z = circ(x, y, D)
    r = sqrt(x.^2+y.^2);
    z = double(r<D/2);
    z(r==D/2) = 0.5;

>> Code used for generating the object

    function y = rect(x, D)
 % function y = rect(x, D)
    if nargin == 1, D = 1; end
     x = abs(x);
     y = double(x<D/2);
    y(x == D/2) = 0.5;

So, my question is, in the original book's code, changing the constants in N, L, D, delta, wvl, z, and W should affect the blurriness of the image. But for me, changing the constant in the W equation (like, say, 0.2 or 1) doesn't seem to influence the blur and the PSF's shape. The other values do affect the result, but not this W equation's constant value.

Can anyone point out where the issue might be in my code? I've attached the relevant parts of the code I'm using.

Answer 1

[Garmit Pant]

Hello

It is my understanding that you are trying to perform defocus aberration using the method based on Chapter 5 of the book 'Numerical Simulation of Optical Wave Propagation with Examples in MATLAB'.

On inspecting your code and comparing it with the method given in the book, the difference seems to be in the calculation of wavefront aberration where you have directly implemented the defocus Zernike polynomial equation.

% wavefront aberration 
W = 0.05 * sqrt(3)*(2*r.^2-1); 

But on inspecting the code given in the book, the calculation of the wavefront aberration and Zernike polynomial has been implemented differently. On changing the code by using the implementation given in the book, the change would be the following:

% wavefront aberration 
% W = 0.05 * sqrt(3)*(2*r.^2-1) 
W = 0.05*sqrt(3)*zrf(2,0,2*r/D); 

The function ‘zrf’ has been implemented as follows:

function R = zrf(n, m, r)  
    R = 0; 
    for s = 0 : ((n-m)/2) 
        num = (-1)^s * gamma(n-s+1); 
        denom = gamma(s+1) * gamma((n+m)/2-s+1)*gamma((n-m)/2-s+1); 
        R = R + num / denom * r.^(n-2*s); 
    end 
end 

The difference in the values of variable ‘W’ calculated is quite significant.

Variable ‘W’ calculated using the method given in the book.

Variable ‘W’ calculated using your method:

The difference in the values of variable ‘W’ is leading to there being minimal change in the generated image and PSF on changing the values of the constants.

I tried changing the value of constants in the method given in the book and saw the following changes:

% wavefront aberration 
% W = 0.05*sqrt(3)*zrf(2,0,2*r/D) 
W = 0.5*sqrt(3)*zrf(2,0,2*r/D)