[PYTHON] Numerical calculation of lens point image distribution function and MTF curve (diffraction calculation)

I recently researched the point spread function (PSF), which often appears in the context of lens design, and the specific numerical calculation method for MTF curves, so I summarized them. I am not a specialist in lens design and I am not from the physics department, so I think there are some mistakes, but I would appreciate it if you could point out any mistakes.

About the point spread function (PSF)

** Point Spread Function (PSF) ** is a function that expresses the light intensity distribution on the image plane when a point light source is imaged on the image plane.

psf.png

The image above is an example of PSF. The closer it is to red, the higher the light intensity. In an ideal imaging, a point light source focuses on one point, so PSF is a function that has a peak at only one point, such as the delta function. However, in reality, due to the effects of aberration and diffraction, the light is not focused on one point, and an image with a certain degree of spread as shown in Fig. 1 can be obtained.

There are two types of PSF: ** geometrical optics PSF **, which does not consider the effect of diffraction, and ** wave optical PSF **, which considers the effect of diffraction. Generally speaking, PSF seems to refer to wave optical PSF. The above example is a wave-optical PSF that takes diffraction into account, with another peak around the center that appears to be wavy.

About MTF curve

The explanation of here is easy to understand, so I will quote it.

MTF (Modulation Transfer Function) is one of the scales for evaluating lens performance, and expresses how faithfully the contrast of the subject can be reproduced on the image plane as a spatial frequency characteristic.

The MTF is generally represented by a graph called the MTF Curve, which looks like this: This graph is also taken from here.

figure_diffraction.gif

The vertical axis represents the contrast, and the horizontal axis represents the image height, which is the distance from the center of the image plane. The red and green lines correspond to specific spatial frequencies, and the solid and dashed lines correspond to the sagittal (S) and meridional (M) directions. I don't really understand the meaning of the sagittal and meridional directions here, but probably the dynamic direction $ r $ and the azimuth direction $ \ theta when the image plane is represented by the polar coordinate system $ (r, \ theta) $. It seems that it corresponds to $.

In summary, the MTF curve represents the contrast in the dynamic and azimuth directions at specific spatial frequencies with varying image heights.

Like PSF, MTF has ** geometrical optics MTF ** and ** wave optical MTF ** depending on whether or not the influence of light diffraction is taken into consideration.

Numerical calculation of PSF

There were two types of PSF: geometrical optics PSF and wave optical PSF. First, let's see how to calculate the geometrical optics PSF.

Numerical calculation of geometrical optics PSF

The geometrical optics PSF is a representation of how a point light source forms an image on the image plane without considering the effects of diffraction. This can be calculated relatively easily using ** ray tracing **.

A point light source is placed at an appropriate position (generally at infinity?), And a large number of rays are emitted from there toward the front of the lens. Each time these rays hit the lens, they are refracted to calculate the next direction, and then the calculation of the collision position with the next lens is repeated. Finally, the position where the light ray intersects the image plane is calculated, and the radiance of the light is added to that position.

By repeating ray tracing for a large number of rays in this way and appropriately interpolating between the sample points, the intensity distribution of the light on the image plane can be obtained. In this way, the geometrical optics PSF can be obtained.

The graph showing the position where the ray intersects the image plane is called ** Spot Diagram **. The figure below is an example.

spot_diagram.png

Diffraction calculation

In order to calculate the wave optical PSF, it is necessary to express light as a wave and obtain the intensity of light on the image plane in consideration of diffraction when it passes through the lens. First, let's look at the basic diffraction calculation method. Finally, I will try numerical calculation using Python.

Here, we will look at the following diffraction calculation methods.

How to express as a wave of light

To think of light as a wave and express it in a mathematical formula, expressions using complex numbers are used.

For plane waves, the complex amplitude $ u (\ boldsymbol {r}, t) $ at position $ \ boldsymbol {r} $, time $ t $ is

u(\boldsymbol{r}, t) = A\exp[i(\boldsymbol{k}\cdot\boldsymbol{r} - \omega t + \delta)]

Is expressed as. Where $ A $ is the amplitude, $ i $ is the imaginary unit, $ \ boldsymbol {k} $ is the wave vector, $ \ omega $ is the angular frequency, and $ \ delta $ is the initial phase.

For spherical waves, the distance from the wave source is $ r $

u(r, t) = \frac{A}{r}\exp[i(kr - \omega t + \delta)]

Is expressed as. Where $ k $ is the wavenumber.

In the calculations below, the $ \ omega t $ and $ \ delta $ terms have no effect and are ignored.

Diffraction integration

In explaining the calculation of diffraction, consider the following coordinate system.

diagram-20200411.png

$ \ Sigma $ represents the opening and $ P (x_p, y_p, 0) $ is the point above $ \ Sigma $. $ Q (x_q, y_q, z_q) $ is a point on the image plane. $ r $ is the length of the line segment $ PQ $, and $ \ theta $ is the angle between the line segment $ PQ $ and the $ z $ axis.

Consider a situation where a plane wave passing through the opening $ \ Sigma $ reaches the point $ Q $ while diffracting. Using the Fresnel-Huygens principle, the complex amplitude of the point $ Q $ can be expressed as a superposition of the secondary waves generated on $ \ Sigma $.

First, consider the complex amplitude of the quadratic wave that reaches the point $ Q $ from the point $ P . $ \frac{1}{i\lambda}\frac{e^{ikr}}{r}u(x_p, y_p, 0)\cos{\theta} $$ It will be. Where $ \ lambda $ is the wavelength of light, $ k $ is the wavenumber, and $ u (x_p, y_p, 0) $ is the complex amplitude at point $ P $.

Since the complex amplitude $ u (x_q, y_q, z_q) $ at the point $ Q $ is the superposition of these secondary waves over $ \ Sigma $.

\begin{align}
u(x_q, y_q, z_q) &= \frac{1}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r}\cos{\theta}dxdy \tag{1} \\
&= \frac{z_q}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r^2}dxdy
\end{align}

It will be. In the expression transformation, we used $ \ cos {\ theta} = \ frac {z_q} {r} $.

This integral is called the ** Fresnel-Kirchhoff diffractive integral **. Since the complex amplitude at the point $ Q $ can be calculated by numerically integrating this, the intensity on the image plane can be obtained by calculating the square of the absolute value. However, this calculation method is not practical because the amount of calculation is $ O (N ^ 4) $ if the number of divisions between the aperture plane and the image plane is $ N $.

Fresnel approximation

Consider a situation where $ \ theta $ is small and $ \ cos {\ theta} \ approx 1 $ (paraxial approximation).

Focusing on Equation 1, the $ \ cos {\ theta} $ part becomes 1, and $ \ frac {1} {r} $ can be approximated to $ \ frac {1} {z} $. On the other hand, the term $ e ^ {ikr} $ has a large value of $ k = \ frac {2 \ pi} {\ lambda} $, so if you replace $ r $ with $ z $, the approximation accuracy will be good. There is none.

Therefore, consider expanding $ r $ into a series and using it up to the second term. When $ r = \ sqrt {(x_q --x) ^ 2 + (y_q --y) ^ 2 + z_q ^ 2} $ is expanded to the second term

r = z_q + \frac{(x_q - x)^2 + (y_q - y)^2}{2z_q} + \cdots

By substituting this into the equation for diffraction integration

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q} \int\int_{\Sigma} u(x, y, 0)e^{ik\frac{(x_q - x)^2 + (y_q - y)^2}{2z_q}}dxdy \tag{2}

It will be. This approximation is called the ** Fresnel approximation **, and Equation 2 is called the ** Fresnel diffraction integral **.

Expanding the contents of the $ e $ shoulder in Equation 2

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x_p, y_p, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}e^{ik\frac{x^2 + y^2}{2z_q}}dx dy \tag{3}

It will be. In Equation 3, $ f_x = \ frac {x_q} {\ lambda z_q} $, $ f_y = \ frac {y_q} {\ lambda z_q} $

\begin{align} 
u(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}e^{ik\frac{x^2 + y^2}{2z_q}}dx dy \\ 
&= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{ik\frac{x^2 + y^2}{2z_q}}e^{-2\pi i(f_x x + f_y y)}dx dy \\ 
&= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{u(x, y, 0)e^{ik\frac{x^2 + y^2}{2z_q}}\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align}

Therefore, we can see that it can be calculated by the two-dimensional Fourier transform of $ u (x, y, 0) e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} $.

In order for the Fresnel approximation to hold, the third term of the series expansion of $ r $ multiplied by $ k $ is required.

\frac{k}{8z_q^3}[(x_q - x)^2 + (y_q - y)^2]^2 \ll 1

Must be met and rewritten for $ z_q $

z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x)^2 + (y_q - y)^2]^2

It will be.

Fraunhofer approximation

In Equation 3, $ \ frac {x ^ 2 + y ^ 2} {2z_q} $ is small enough and approximates $ e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} \ approx 1 $ Then

u(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-ik\frac{xx_q + yy_q}{z_q}}dxdy \tag{4}

It will be. This approximation is called the ** Fraunhofer approximation **, and Equation 4 is called the ** Fraunhofer diffraction integral **.

In Equation 4, let $ f_x = \ frac {x_q} {\ lambda z_q} $, $ f_y = \ frac {y_q} {\ lambda z_q} $

\begin{align} u(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} u(x, y, 0)e^{-2\pi i(f_x x + f_y y)}dx dy \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{u(x, y, 0)\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}}
\end{align}

Therefore, we can see that it can be calculated by the two-dimensional Fourier transform of $ u (x_p, y_p, 0) . For the Fraunhofer approximation to hold, it is necessary to satisfy $ \ frac {k} {2z_q} (x ^ 2 + y ^ 2) \ ll 1 $$, which can be rewritten for $ z_q $ to get $ z_q \ It becomes gg \ frac {\ pi} {\ lambda} (x ^ 2 + y ^ 2) $.

The Fraunhofer approximation is a valid approximation when the aperture plane is small enough and the image is observed very far from the aperture plane. On the other hand, the Fresnel approximation is a valid approximation even closer than the Fraunhofer approximation.

Angle spectral method

The inverse Fourier transform of $ U (f_x, f_y, 0) $, which is the Fourier transform of the complex amplitude $ u (x, y, 0) $ of the aperture surface, is

u(x, y, 0) = \int\int U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)]df_x df_y

is. In this equation, the complex amplitude $ u (x, y, 0) $ of the opening surface is the sum of the plane wave $ U (f_x, f_y, 0) \ exp [2 \ pi i (x f_x + y f_y)] $. It shows that it can be expressed as.

Finding each component of the wave vector $ \ boldsymbol {k} $ of each plane wave

\begin{align}
k_x &= 2\pi f_x \\
k_y &= 2\pi f_y \\
k_z &= \sqrt{k^2 - k_x^2 - k_y^2}
\end{align}

It will be.k_zIs\|\boldsymbol{k}\| = kI will obey from.

When each plane wave travels $ z_q $ in the $ z $ direction, its complex amplitude is

U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)] \exp[i k_z z_q]

It will be. Then the complex amplitude $ u (x', y', z_q) $ on the image plane is the sum of those plane waves.

\begin{align}
u(x', y', z_q) &= \int\int U(f_x, f_y, 0)\exp[2\pi i(x f_x + y f_y)] \exp[i k_z z_q] df_x df_y \\
&= \mathcal{F}^{-1}\{ U(f_x, f_y, 0) \exp[i k_z z_q] \}
\end{align}

It can be seen that the complex amplitude $ u (x', y, z_q) $ on the image plane can be calculated by the inverse Fourier transform of $ U (f_x, f_y, 0) \ exp [i k_z z_q] $. This calculation method is called ** Angular Spectrum Method **.

The angular spectral method can calculate diffraction in a region closer to the aperture surface than the Fresnel approximation or Fraunhofer approximation. Moreover, since the calculation is performed only once for the Fourier transform and the inverse Fourier transform, it can be calculated at high speed.

Numerical calculation example of diffraction

Aperture $ \ Sigma $ with radius $ r = 1 ~ \ mathrm {[mm]} $ circular opening, $ \ lambda = 550 ~ \ mathrm {[nm]} $, $ z_q = 5 ~ \ mathrm {[m] Let's calculate numerically as} $. You can see the Jupyter Notebook used for the calculation at Google Colab.

Numerical calculation of Fraunhofer diffraction

First, let's calculate numerically with Fraunhofer approximation. If you check the approximation conditions, $ x $ and $ y $ will be $ r $ at the maximum.

z_q \gg \frac{\pi}{\lambda}r^2 \approx 5.71

It will be. It does not meet the hypothetical conditions, but I will proceed with the calculation as it is.

Fraunhofer diffraction can be calculated by performing a two-dimensional Fourier transform on the complex amplitude $ u (x, y, 0) $ on the aperture $ \ Sigma $. To perform numerical calculations programmatically, divide the calculation area with a side length of $ D $ into a grid of $ N \ times N $ so as to cover the opening surface, and store the complex amplitude on each grid. Can be calculated by performing a two-dimensional discrete Fourier transform.

Let's calculate numerically using Python. First, prepare an array that stores the complex amplitude on the grid.

import numpy as np
import matplotlib.pyplot as plt

l = 0.55e-6 #wavelength[m]
k = 2 * np.pi / l #Wave number
zq = 5 #Distance from the aperture plane to the image plane[m]
N = 200 #Number of grid divisions
r = 0.001 #Aperture radius[m]
D = 0.05 #Computation area on the open surface[m]

fscale = l*zq

#Returns the complex amplitude on the open plane
def P(x, y):
    if x**2 + y**2 < r**2:
        return 1
    else:
        return 0

#Grid generation
X, Y = np.meshgrid(np.linspace(-D/2, D/2, num=N), np.linspace(-D/2, D/2, num=N))
Pv = np.vectorize(P)
Z = Pv(X, Y)

Let's plot it.

plt.imshow(Z, extent=[-D/2, D/2, -D/2, D/2], cmap='gray')
plt.xlabel('$x_p$ [m]')
plt.ylabel('$y_p$ [m]')

download.png

The white part corresponds to the opening surface.

Next, this is subjected to the discrete Fourier transform, but before that, the area on the image plane corresponding to the calculation area is calculated.

v = np.fft.fftfreq(N, d=D/N)
v = fscale * np.fft.fftshift(v)
#Units[mm]
extent= [1e3 * v[0], 1e3 * v[-1], 1e3 * v[0], 1e3 * v[-1]]

The Fraunhofer diffraction image is obtained by discrete Fourier transforming the complex amplitude on the aperture plane.

#Calculate complex amplitude on image plane
U_fraun = np.fft.fft2(Z)
U_fraun = 1/(1.0j * l * zq) * np.exp(1.0j * k * zq) * np.exp(1.0j * k * (v**2 + v**2)/(2*zq)) * np.fft.fftshift(U_fraun)

#Calculate strength
I_fraun = np.abs(U_fraun)**2
#Normalization
I_fraun = I_fraun / np.max(I_fraun)

Let's plot the intensity.

plt.imshow(I_fraun, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fraunhofer Diffraction')

download (1).png

There is also an intensity peak around the center, which shows that the effect of diffraction can be calculated.

Numerical calculation of Fresnel diffraction

First, check the conditions of approximation. The conditional expression for the approximation of Fresnel diffraction is

z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x)^2 + (y_q - y)^2]^2

was. Unlike Fraunhofer diffraction, it is a little complicated because the size of the observation surface must be taken into consideration. For the time being, the size of one side of the observation surfacelI will leave it as. Then|x_q - x| \le r + \frac{l}{2}、yIs the same, so if you substitute this

z_q \gg \sqrt[3]{\frac{\pi}{\lambda}(r + \frac{l}{2})^4} \approx 0.21

It will be.

Fresnel diffraction can be calculated using the two-dimensional discrete Fourier transform as well as Fraunhofer diffraction.

#Calculate complex amplitude on image plane
U_fresnel = np.fft.fft2(Z * np.exp(1.0j * k * (X**2 + Y**2)/(2*zq)))
U_fresnel = 1/(1.0j * l * zq) * np.exp(1.0j * k * zq) * np.exp(1.0j * k * (v**2 + v**2)/(2*zq)) * np.fft.fftshift(U_fresnel)

#Calculate strength
I_fresnel = np.abs(U_fresnel)**2
#Normalization
I_fresnel = I_fresnel / np.max(I_fresnel)

Let's plot the intensity.

plt.imshow(I_fresnel, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Diffraction')

download (2).png

It seems that the peak around the center is stronger than in the Fraunhofer approximation.

Direct calculation of Fresnel Kirchhoff's diffraction formula

Fresnel Kirchhoff's diffractive formula

\frac{z_q}{i\lambda}\int\int_{\Sigma} u(x, y, 0)\frac{e^{ikr}}{r^2}dx dy

Let's calculate by directly integrating numerical values. This method is relatively rigorous because it does not use approximations, but it is not practical due to the overwhelming amount of computation. I wanted to compare it with the approximation method, so I calculated it.

#Fresnel Kirchhoff's diffractive numerical integration
def Uf(x_q, y_q):
    r = np.sqrt((x_q - X)**2 + (y_q - Y)**2 + zq**2)
    return zq/(1.0j * l) * (Z * np.exp(1.0j * k * r) / r**2).sum()

#Computational area on the image plane
X_Q, Y_Q = np.meshgrid(np.linspace(extent[0]*1e-3, extent[1]*1e-3, num=N), np.linspace(extent[2]*1e-3, extent[3]*1e-3, num=N))
#Numerical integration(Scaling is appropriate)
U_strict = np.vectorize(Uf)(X_Q, Y_Q)

#Calculate strength
I_strict = np.abs(U_strict)**2
#Normalization
I_strict = I_strict / np.max(I_strict)

Let's plot the intensity.

plt.imshow(I_strict, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Kirchhoff Diffraction')

download (3).png

Numerical calculation by the angular spectrum method

Finally, let's perform a numerical calculation using the angular spectral method. The procedure is to first find the Fourier transform $ U (f_x, f_y, 0) $ of the complex amplitude $ u (x, y, 0) $ of the opening surface, and then $ U (f_x, f_y, 0) $ and $ . Inverse Fourier transform is performed by multiplying exp [ik_z z_q] $.

#Fourier transform the complex amplitude of the open surface
U = np.fft.fft2(Z)

#Wave vector calculation
k_x = np.fft.fftfreq(N, d=D/N) * 2 * np.pi
k_y = np.fft.fftfreq(N, d=D/N) * 2 * np.pi
K_X, K_Y = np.meshgrid(k_x, k_y)
k_z = np.sqrt(k**2 - K_X**2 - K_Y**2)

#Calculate the complex amplitude of the image plane
U_angular = np.fft.ifft2(U * np.exp(1.0j * k_z * zq))

#Calculate strength
I_angular = np.abs(U_angular)**2
#Normalization
I_angular = I_angular / np.max(I_angular)

Let's plot it.

plt.imshow(I_angular, cmap='jet', extent=[-D/2, D/2, -D/2, D/2])
plt.xlabel('$x_q$ [m]')
plt.ylabel('$y_q$ [m]')
plt.title('Angular Spectrum Diffraction')

download (4).png

Unlike the case of Fresnel approximation and Fraunhofer approximation, the area of the image plane coincides with the calculation area, so it looks small.

It's been long, so I'll write the rest on another page.

References etc.

Recommended Posts

Numerical calculation of lens point image distribution function and MTF curve (diffraction calculation)
[Introduction to Scipy] Calculation of Lorenz curve and Gini coefficient ♬