Point Spread Function of an optical system

Hi all,

I need to estimate the PSF of an optical system. I have all the parameters and have been following an excellent code that was mentioned in the 'Spherical aberration and Chromatic aberration' question ( https://de.mathworks.com/matlabcentral/answers/36064-spherical-aberration-and-chromatic-aberration ). I am trying to model the PSF of an optical system and then apply it to an image to see how blurry the image would be.

Everything works great but the only problem is the pixel sampling. With the code I am having trouble with the psf_sampling. My psf_sampling is larger (10e-6) and every time I run the code I get the following error:

'Sampling is not sufficient to reconstruct the entire wavefront.'

Please help me with this, much appreciated.

Kind regards,

Vanessa

%%Define parameters
clear all;close all;
psf_sampling = 0.5e-6;%focal plane sampling in meters
lambda = 0.6328e-6;%wavelength in meters
N = 256;
aperture_diameter = 0.0254;%meters; 1 inch
focal_length = 5*aperture_diameter;%meters
RMS_SA = 0.25;%RMS spherical aberration content in waves
%%Calculate pupil plane sampling
delta_fx = 1/(psf_sampling*N);
x_pupil = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length;
[X_pupil,Y_pupil] = meshgrid(x_pupil);
R_pupil = sqrt(X_pupil.^2 + Y_pupil.^2);
R_norm = R_pupil/(aperture_diameter/2);%Pupil normalized to aperture diameter
assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');
%%Create wavefront
W = RMS_SA * sqrt(5) * (6*R_norm.^4 - 6*R_norm.^2 + 1);%Spherical Aberration wavefront
W(R_norm>1) = 0;
E = exp(1i*2*pi*W);%Complex amplitude
E(R_norm>1) = 0;%Impose aperture size
figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');
%%Create point-spread function
psf = abs(fftshift(fft2(ifftshift(E)))).^2;
psf = psf/sum(psf(:));%Normalize to unity energy
x_psf = (-fix(N/2):fix((N-1)/2)) * psf_sampling;
figure;imagesc(x_psf*1e6,x_psf*1e6,psf);
title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));

Answer 1

[Eric]

Sorry, I was out on vacation for a couple weeks, then I posted to

https://www.mathworks.com/matlabcentral/answers/36064-spherical-aberration-and-chromatic- aberration

Here's that response:

The assert statement is checking that the sampling allows reconstruction of the entire wavefront. If this assertion fails, then the pupil sampling is such that you cannot reconstruct the entirety of the wavefront area. This can happen if the pixel pitch is relatively large. Assuming you don’t want to adjust the aperture diameter, you need to increase the values in R_norm.

For instance, if I set the psf_sampling to 10e-6 meters, the code from the page results in max(R_norm(:)) of 0.4475.

A simple fix is to calculate the PSF at a finer resolution and then bin that PSF to the desired low-resolution. So for the case where you want a PSF sampling of 10 microns, you could calculate the PSF at a sampling of 2.0 microns. Calculate the psf as shown in the script. Then do the following:

kernel = ones(5,5)/5^2;
psf_lowres_all = conv2(psf, kernel, 'same');
psf_lowres = psf_lowres_all(3:5:end,3:5:end);

In this case you’re doing convolution to spatially sum 5x5 neighborhoods of the high-resolution PSF (5 is the correct value because it’s the ratio of the desired low-resolution pixel spacing to the calculated spacing). This convolution produces 25 possible low-resolution PSFs, each at a different registration of the optical PSF to the pixel grid. In the last line I’ve selected the one in the middle. In the indexing you skip by 5 because you want to get sums of adjacent blocks (i.e., not sliding blocks). You can change the indexing to select another possible low-resolution PSF. For instance, you could also use

psf_lowres = psf_lowres_all(1:5:end,4:5:end);

to select another PSF.

You’ll notice when this is done that the wavefront calculation area is zero-padded to be larger than the specified system aperture. You can check this by plotting abs(E).

[Eric]

By the way, the system you specified has a Q of 0.5, where Q is lambda*f_number/pixel_pitch. This is a fairly common value for imaging systems. The code I wrote will work with a Q of 1 or higher, but you still might worry about aliasing. I would still use a psf_sampling of 2 microns to get the PSF at 5X resolution and then downsample from there.

A good reference on Q, though he doesn't use that term, is the following paper: R. D. Fiete, "Image quality and lambda FN/p for remote sensing systems," Optical Engineering 38, 1229-1240 (1999). This includes a nice discussion of why most imaging systems are designed with Q values less than 2 (a system is Nyquist sampled when Q=2, super-sampled when Q>2).

[Eric]

One more suggestion:

If possible, I would do the entire image simulation in the high-resolution space and then downsample the blurred image.

That is, if you're working with simulated imagery, simulate the imagery at the 2 micron sampling as well (or any value less than 2.5 such that the Q of your system is 2 or higher). Convolve the high-resolution image with the high-resolution PSF. Then do the convolution with the kernel function and extraction of a representative image. This isn't always possible (e.g., the image you're blurring is measured), but if it is possible it most accurately represents information at aliased spatial frequencies.

[Vanessa]

Thank you so much Eric!!! I really appreciate this, it will definetly help with my thesis! I also hope you had a good vocation!

Answer 2

[Image Analyst]

Decrease psf_sampling back to the range where it is valid. If you sample the PSF too coarsely, you're not going to be able to get it's shape at all.

[Vanessa]

Thank you your reply, however I cannot change the parameters as they are set for the optical system, is there any way to estimate the PSF with the parameters I have?

[Image Analyst]

According to the program, no. Why do you want to decrease the resolution so much?

[Vanessa]

Thank you for your answer Image Analyst, much appreciated! However, it's not that I want to decrease the resolution but these are the parameters for the optical system that I have been given and cannot be changed:

• wavelength = 2.5e-6 m
• image size = 512 x 512
• f number = 2
• aperture diameter = 0.03 m
• focal length = 0.06 m
• pixel sampling = 10 e-6 m

[Image Analyst]

I don't know without delving into the program and parameters much more. Just try commenting out the assert statement and see how it goes.

[Vanessa]

Hi Image Analyst,

I have tried everything and no luck, but thank you for your replies, I am truly grateful and if by chance you figure something out, please let me know.

[Eric]

If you comment out the assert statement, you'll get poor results. You will reconstruct a wavefront that is smaller than the aperture of your system.

Answer 3

[Tilkesh]

If we just take thin lens function to find the PSF then this code does not work. (No aberrations simple lens equation)

For any simple lens with diameter=25mm, focal lengh=100mm. How to determine sampling for this code. What is the criteria? I choose Q>2 as mentioned by author still this code does not work.

If we want to find the Airy disc diameter from this code then it should be independent of PSF sampling and should only depend on wavelength and aperture diameter but it does not work why? I am unable to find it. Is it the physical problem of computentional or both still not clear to me. Please help for the benifit of many many reserachers.

%%Define parameters

clear all;close all;

psf_sampling = 0.5e-6;%focal plane sampling in meters

lambda = 0.6328e-6;%wavelength in meters

N = 256;

aperture_diameter = 0.0254;%meters; 1 inch

focal_lengt%%Define parameters

clear all;close all;clc;

psf_sampling = .01e-6;%focal plane sampling in meters

lambda = 0.6328e-6;%wavelength in meters

N = 256;

aperture_diameter = 0.0254;%meters; 1 inch

focal_length = 5*aperture_diameter;%meters

RMS_SA = 0.25;%RMS spherical aberration content in waves

%%Calculate pupil plane sampling

delta_fx = 1/(psf_sampling*N);

x_pupil = (-fix(N/2):fix((N-1)/2)) * delta_fx * lambda * focal_length;

[X_pupil,Y_pupil] = meshgrid(x_pupil);

R_pupil = sqrt(X_pupil.^2 + Y_pupil.^2);

R_norm = R_pupil/(aperture_diameter/2);%Pupil normalized to aperture diameter

assert(max(R_norm(:))>=sqrt(2),'Sampling is not sufficient to reconstruct the entire wavefront.');

%%Create wavefront

W =(-(2*pi/(lambda*focal_length))*((X_pupil.^2 + Y_pupil.^2)));

E = exp(1i*W);%Complex amplitude

figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');

%%Create point-spread function

psf = abs(fftshift(fft2(ifftshift(E)))).^2;

psf = psf/sum(psf(:));%Normalize to unity energy

x_psf = (-fix(N/2):fix((N-1)/2)) * psf_sampling;

figure;imagesc(x_psf*1e6,x_psf*1e6,psf);

title(sprintf('PSF with %.4f waves RMS of Spherical Aberration', RMS_SA));

xlabel(sprintf('Position (\\mum)'));

ylabel(sprintf('Position (\\mum)'));

Q=(lambda*focal_length)./(aperture_diameter*psf_sampling)