# spherical aberration and Chromatic aberration

172 views (last 30 days)
Alan AF on 20 Apr 2012
Commented: VK on 26 Jul 2022
Hi
How can we simulate spherical aberration in matlab? It might be by using the below code, but what are the parameters ??
PSF = fspecial(?); % create PSF
Blurred = imfilter(I,PSF,?,'conv');
Thank you
##### 2 CommentsShow 1 older commentHide 1 older comment
VK on 26 Jul 2022
I am trying to do convnfft of psf and image for convolution, but unfortunately it shows error. Have anyone face such problem?

Eric on 20 Apr 2012
The following code lets you simulate the PSF associated with spherical aberration (as well as diffraction). You need to specify the PSF sampling pitch, the wavelength, the aperture diameter, the system focal length, the amount of spherical aberration, and the PSF array size.
The idea is to calculate a wavefront with spherical aberration and then use the Fourier transform to calculate the PSF. You need to be careful with sampling, hence the assert() statement.
You can easily include other aberrations as well, if desired. Simply modify how the wavefront is calculated.
You can check that this code is working by setting RMS_SA to zero (i.e., simulate a diffraction-limited system, in which case the PSF is an Airy disk). Then the first null of the PSF should have a radius of 1.22*lambda*focal_length/aperture_diameter.
You'll have to do some work to simulate a spectrally broad system. In that case you could simulate a number of monochromatic PSFs and then sum them together. You would need to know how the spherical aberration varies as a function of wavelength (variation as the reciprocal of wavelength would be a good estimate - i.e., longer wavelengths have less aberration) as well as the spectral output of the source, the spectral transmittance of the optical system, and the spectral responsivity of the detector. The product of these spectral terms would be used as the weights in a weighted average to calculate the broadband PSF from the monochromatic PSFs.
Once you get the PSF you can convolve this with your scene to simulate imaging performance. imfilter() may work, but I like convnfft ( http://www.mathworks.com/matlabcentral/fileexchange/24504-fft-based-convolution) from the File Exchange. You could use conv2() as well.
Good luck,
Eric
%%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)');
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)'));
Eric on 26 Nov 2013
That's right. Add other aberrations to the variable W and you should be set.
-Eric

Eric, hi again; The all problem is, I want to simulate aberrations by Zernike polynomials to 15th sentence, for a 16" telescope, but i think, somewhere there are some defeats on code, could you please help me to solve the problem? this is the code:
clear all;
close all;
psf_sampling = 0.5e-6;%focal plane sampling in meters
lambda = 0.6328e-6;%wavelength in meters
N=256;
teta=45;% the angle to compute
aperture_diameter = 0.407;%meters; 16 inch
focal_length = 5*aperture_diameter;%meters
%RMS_SA = 0.15;%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 & Aberations
%z0 is equal to 1
z1 = (2 * R_norm) *cos(teta);%Tip Aberration
z1(R_norm>1)=0;
z2 = (2*R_norm) *sin(teta);%Tilt Aberration
z2(R_norm>1)=0;
z3 = sqrt(3) * (2* R_norm.^2 - 1);% De-focus Aberration
z3(R_norm>1)=0;
z4 = sqrt(6) * (R_norm.^2) * sin(2*teta);% Astigmatism Aberration
z4(R_norm>1)=0;
z5 = sqrt(6) * (R_norm.^2) * cos(2*teta);% Astigmatism Aberration
z5(R_norm>1)=0;
z6 = sqrt(8) * (3 * R_norm.^3 - 2 * R_norm) * sin(teta);% Coma Aberration
z6(R_norm>1)=0;
z7 = sqrt(8) * (3 * R_norm.^3 - 2 * R_norm) * cos(teta);% Coma Aberration
z7(R_norm>1)=0;
z8 = sqrt(8) * (R_norm.^3) * sin(3*teta);% Trefoil aberration
z8(R_norm>1)=0;
z9 = sqrt(8) * (R_norm.^3) * cos(3*teta);% Trefoil aberration
z9(R_norm>1)=0;
z10 = sqrt(5) * (6 * R_norm.^4 - 6 * R_norm.^2 + 1);%Spherical Aberration wavefront
z10(R_norm>1) = 0;
z11 = sqrt(10) * (4 * R_norm.^4 - 3 * R_norm.^2) * cos(2 * teta);% ----- Aberration
z11(R_norm>1) = 0;
z12 = sqrt(10) * (4 * R_norm.^4 - 3 * R_norm.^2) * sin(2 * teta);% ----- Aberration
z12(R_norm>1) = 0;
z13 = sqrt(10) * (R_norm.^4) * cos(4 * teta)% ----- Aberration
z13(R_norm>1) = 0;
z14 = sqrt(10) * (R_norm.^4) * sin(4 * teta)% ----- Aberration
z14(R_norm>1) = 0;
E= exp(1i*2*pi*z1) + exp(1i*2*pi*z2) + exp(1i*2*pi*z3) + exp(1i*2*pi*z4) + exp(1i*2*pi*z5) + exp(1i*2*pi*z6) + exp(1i*2*pi*z7) + exp(1i*2*pi*z8) + exp(1i*2*pi*z9) + exp(1i*2*pi*z10) + exp(1i*2*pi*z11) + exp(1i*2*pi*z12) + exp(1i*2*pi*z13) + exp(1i*2*pi*z14);%Complex amplitude for first 15 zernike polynomials
E(R_norm>1) = 0;%Impose aperture size
figure;imagesc(angle(E)/(2*pi));colorbar;title('Wavefront Phase (waves)');
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 Aberration'));
xlabel(sprintf('Position (\\mum)'));
ylabel(sprintf('Position (\\mum)'));
I write only for specific angle(45), and the weights for all aberrations are equal to 1. Could u tell me this code is OK or not?
##### 2 CommentsShow 1 older commentHide 1 older comment
Hi,
DEAR ERIC
First, Excuse me. I only check the mathworks page.
You are right the coefficients extremely effect on aberrations, but for the first step i used 1 for all of them. Yesterday i found a function, Created by Dick Brunson, on 12/1/97. This code generates all Zernike polynomial, and size of row or col of the zern output array.
The code is here:
%
% function to generate a matrix of phase values from a single
% Zernike Polynomial coefficient. Enter genzern with no parameters
% to view the definitions of first 15 Zernike Polynomial terms.
%
% Created by:
% Dick Brunson
% 12/1/97
%
% Inputs:
%
% nz number of Zernike Coefficient desired(if <0 print out
% definition of Zernike Polynomial
% arraysize size of row or col of the zern output array.
%
% mask[optional] array of 1's & 0's to define boundaries of Zernike
% output array. If undefined, output array is max circle
% that can be inscribed inside the array.
%
% Outputs:
% zern output array of phase values
% mask pupil function or definition of good phase point area
% (defined by mask input or default)
%
% Zernike Polynomials definiton from "Optical Shop Testing", Daniel Malacara, 1992
% pg 465
%
if nargin < 1
disp('Zernike # n m Zernike Definition')
disp(' 1 0 0 1 Piston')
disp(' 2 1 0 r*sin(phi) Tilt Y')
disp(' 3 1 1 r*cos(phi) Tilt X')
disp(' 4 2 0 r^2*sin(2*phi) Astig 1st ord. 45 deg')
disp(' 5 2 1 2*r^2 - 1 Defocus')
disp(' 6 2 2 r^2*cos(2*phi) Astig 1st ord. 0 deg')
disp(' 7 3 0 r^3*sin(3*phi) Trifoil 30 deg')
disp(' 8 3 1 (3*r^3 - 2*r)sin(phi) Coma Y')
disp(' 9 3 2 (3*r^3 - 2*r)cos(phi) Coma X')
disp(' 10 3 3 r^3*cos(3*phi) Trifoil 0 deg')
disp(' 11 4 0 r^4*sin(4*phi) Tetrafoil 22.5 deg')
disp(' 12 4 1 (4*r^4 - 3*r^2)sin(2*phi) Astig 2nd ord 45 deg')
disp(' 13 4 2 6*r^4 - 2*r^2 - 1 Spherical Aberration')
disp(' 14 4 3 (4*r^4 - 3*r^2)cos(2*phi) Astig 2nd ord 0 deg')
disp(' 15 4 4 r^4*cos(4*phi) Tetrafoil 0 deg')
return
end
xv = [1:arraysize]';
yv = [1:arraysize];
x = [];
y = [];
for i = 1:arraysize
x = [x, xv];
y = [y; yv];
end
center = 0.5 + arraysize/2;
r = sqrt((x - center).^2 + (y - center).^2);
r = r/(center - 1);
if nargin < 3
% define inscribed circle for mask
router = r(arraysize/2,1);
rmax_new = max(max(r))/router;
r = r * rmax_new/max(max(r));
else
r = r/router;
end
phi = atan2((x - center),(y - center));
zern =zeros(arraysize);
% compute indices m & n for Zernikes per definition
csum = cumsum(0:nz);
nindx = sum(csum < nz);
mx = csum(nindx);
n = nindx - 1;
m = nz - mx - 1;
%
% Zernike generating function from "Optical Shop Testing", Daniel Malacara, 1992
%
% ___m
% (n-2m) \ s (n-2s)
% R = /__ (-1) [(n-s)!/(s!(m-s)!(n-m-s)!)] * r
% n s=0
% ___m
% \ q(s)
% = /__ p(s) * r
% s=0
%
% Where the Zernike Polynomial is defined by:
%
% ___k ___n
% \ \ n-2m sin(n-2m)phi
% W(r,phi) = /__ /__ Anm Rn or
% n=0 m=0 cos(n-2m)phi
%
% where sine function is used for n-2m>0 & cosine for n-2m<=0
%
for s = 0:m
if (n - s) < 0 | (m - s) < 0 | (n - m - s) < 0
break
end
p = (-1)^s * prod(1:(n - s))/(prod(1:s) * prod(1:(m - s)) * prod(1:(n - m - s)));
q = n - 2*s;
zern = zern + p * r .^ q;
% disp(num2str([nz,n,m,n-2*m,p,q]))
end
if n - 2*m > 0
zern = zern .* sin((n - 2*m)*phi);
else
zern = zern .* cos((n - 2*m)*phi);
end
return
I tried to use this function on simulation code, so we can select the number of Zernike Coefficient (and the aberrations) desired in code, and increase the precision of Zernike polynomials, on wavefront. But there are some issues, the first is how to relate aperture diameter to arraysize, the second is how to apply the coefficient weights? For the second issue, i think we can insert the weights by a "For Loop" in the function.
thanks again

Eric on 3 Dec 2013
I'll have to take a look at this code later. I am reasonably sure that Malacara's Zernikes use a different ordering and scaling coefficient than the "Noll" Zernikes I use. That's one thing to be careful of. I'll have to look through it carefully to see how to apply the aperture size.
You might also be interested in the article at:
-Eric
##### 2 CommentsShow 1 older commentHide 1 older comment
Walter Roberson on 7 Apr 2017
dong,
Eric will not be notified about new comments, so I recommend that you contact him by way of his profile https://www.mathworks.com/matlabcentral/profile/authors/210783-eric

Eric on 3 Dec 2013
I looked over this function quickly and I think you want to set the "mask" input value of genzern() to R_norm<=1 from my code.
To implement this code you'll need to implement a for loop. The input argument nz is a scalar for the Zernike index and so this function returns only a single wavefront. You'll need to loop over the indices for the polynomials of interest and multiply the zern output of the function by the appropriate scalar. Sum these scaled wavefronts together to get the total wavefront.
Again, be careful about scaling. The normalization of these Zernikes is different than the ones you have in your earlier post. See the caveat in my original answer about being very specific about how Zernike coefficients are defined. There are lots of definitions. Malacara's are as valid as others, but just always make sure you know which definition you're working with. For what it's worth, I prefer Noll's definition from R. J. Noll, "Zernike Polynomials and Atmospheric-Turbulence," J Opt Soc Am 66, 207-211 (1976). The nice thing about this definition is that the root-sum-square of the Zernike coefficients is equal to the RMS wavefront error for an unobscured system. This makes comparing Zernike magnitudes of different terms considerably easier.
-Eric
##### 2 CommentsShow 1 older commentHide 1 older comment
zzzj1208 on 3 Mar 2022
Hi, anyone?
I still have questions about how to solve the RMS. What's the meaning of the 'Let obj.Map be the phase of the surface in waves. ' and how to get ogj.Map?
Or how to find RMS_SA for a system with only one lens? Any help would be greatly appreciated!