J'ai récemment étudié la fonction de distribution d'image ponctuelle (PSF), qui apparaît souvent dans le contexte de la conception de lentilles, et la méthode de calcul numérique spécifique de la courbe MTF, je les ai donc résumées. Je ne suis pas un spécialiste de la conception de lentilles et je ne suis pas du département de physique, donc je suis sûr qu'il y a des erreurs, mais j'apprécierais que vous puissiez signaler les erreurs.
** Fonction de diffusion ponctuelle --PSF ** est une fonction qui exprime la distribution d'intensité lumineuse sur le plan image lorsqu'une source lumineuse ponctuelle est imagée sur le plan image.
L'image ci-dessus est un exemple de PSF. Plus il est proche du rouge, plus l'intensité lumineuse est élevée. Dans l'imagerie idéale, une source de lumière ponctuelle se concentre sur un point, donc PSF est une fonction qui a un pic en un seul point, telle que la fonction delta. Cependant, en réalité, en raison des effets d'aberration et de diffraction, l'image n'est pas focalisée sur un point, et une image avec un certain degré d'étalement comme le montre la figure 1 peut être obtenue.
Il existe deux types de PSF: ** PSF optique géométrique ** qui ne prend pas en compte l'effet de diffraction, et ** PSF optique d'onde ** qui considère l'effet de diffraction. D'une manière générale, PSF semble signifier PSF optique à ondes. L'exemple ci-dessus est une PSF à onde optique qui prend en compte la diffraction, avec un autre pic autour du centre qui semble ondulé.
L'explication de ici est facile à comprendre, je vais donc la citer.
MTF (Modulation Transfer Function) est l'une des échelles d'évaluation des performances de l'objectif, et exprime avec quelle fidélité le contraste du sujet peut être reproduit sur le plan de l'image en tant que caractéristique de fréquence spatiale.
MTF est généralement représenté par un graphique appelé MTF Curve comme suit. Ce graphique est également tiré de ici.
L'axe vertical représente le contraste et l'axe horizontal représente la hauteur de l'image, qui est la distance par rapport au centre du plan d'image. Les lignes rouges et vertes correspondent à des fréquences spatiales spécifiques, et les lignes pleines et discontinues correspondent à la direction sagittale (S) et à la direction méridienne (M). Je ne comprends pas vraiment la signification des directions sagittale et méridienne ici, mais probablement la direction dynamique $ r $ et la direction de l'angle azimutal $ \ theta lorsque le plan image est représenté par le système de coordonnées polaires $ (r, \ theta) $. Il semble que cela corresponde à $.
En résumé, la courbe MTF représente le contraste dans les directions dynamique et azimutale à des fréquences spatiales spécifiques avec des hauteurs d'image variables.
Comme le PSF, le MTF a ** MTF optique géométrique ** et ** MTF optique d'onde ** selon que l'influence de la diffraction de la lumière est prise en compte ou non.
Il existe deux types de PSF: PSF optique géométrique et PSF optique ondulatoire. Voyons d'abord comment calculer le PSF optique géométrique.
La PSF géométrique et optique est une représentation de la façon dont une source de lumière ponctuelle forme une image sur le plan image sans tenir compte des effets de diffraction. Cela peut être calculé relativement facilement en utilisant le ** traçage de rayons **.
Une source lumineuse ponctuelle est placée à une position appropriée (généralement à l'infini?), Et un grand nombre de rayons sont émis à partir de là vers l'avant de la lentille. Chaque fois que ces rayons frappent la lentille, ils sont réfractés pour calculer la direction suivante, puis le calcul de la position de collision avec la lentille suivante est répété. Enfin, la position où le faisceau lumineux coupe le plan image est calculée et la luminosité de rayonnement de la lumière est ajoutée à cette position.
En répétant le laçage tardif pour un grand nombre de rayons de cette manière et en interpolant de manière appropriée entre les points d'échantillonnage, la distribution d'intensité de la lumière sur le plan image peut être obtenue. De cette manière, vous pouvez obtenir un PSF optique géométrique.
Le graphique montrant la position où le faisceau lumineux coupe le plan image est appelé ** Diagramme de points **. La figure ci-dessous est un exemple.
Pour calculer l'onde optique PSF, il est nécessaire d'exprimer la lumière sous forme d'onde et d'obtenir l'intensité de la lumière sur le plan image en tenant compte de la diffraction lors du passage à travers la lentille. Tout d'abord, regardons la méthode de base de calcul de diffraction. Enfin, je vais essayer le calcul numérique avec Python.
Ici, nous examinerons les méthodes de calcul de diffraction suivantes.
Pour considérer la lumière comme une onde et l'exprimer dans une formule mathématique, on utilise l'expression utilisant des nombres complexes.
Pour les ondes planes, l'amplitude complexe $ u (\ boldsymbol {r}, t) $ à la position $ \ boldsymbol {r} $, le temps $ t $ est
Est exprimé comme. Où $ A $ est l'amplitude, $ i $ est l'unité imaginaire, $ \ boldsymbol {k} $ est le vecteur du nombre d'onde, $ \ omega $ est la fréquence angulaire et $ \ delta $ est la phase initiale.
Pour les ondes sphériques, la distance de la source d'onde est $ r $
Est exprimé comme. Où $ k $ est le nombre de vagues.
Dans les calculs suivants, les termes $ \ omega t $ et $ \ delta $ n'ont aucun effet et sont ignorés.
Pour expliquer le calcul de la diffraction, considérez le système de coordonnées suivant.
$ \ Sigma $ représente l'ouverture et $ P (x_p, y_p, 0) $ est le point au-dessus de $ \ Sigma $. $ Q (x_q, y_q, z_q) $ est un point sur le plan image. $ r $ est la longueur de la ligne $ PQ $, et $ \ theta $ est l'angle entre la ligne $ PQ $ et l'axe $ z $.
Considérons la situation où une onde plane passant par l'ouverture $ \ Sigma $ atteint le point $ Q $ en diffractant. En utilisant le principe de Fresnel-Huihens, l'amplitude complexe du point $ Q $ peut être exprimée comme une superposition d'ondes secondaires générées sur $ \ Sigma $.
Tout d'abord, considérons l'amplitude complexe de l'onde quadratique qui atteint le point $ Q $ à partir du point $ P
Puisque l'amplitude complexe $ u (x_q, y_q, z_q) $ au point $ Q $ est la superposition de ces ondes secondaires sur $ \ 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}
Ce sera. Dans la transformation de formule, nous avons utilisé $ \ cos {\ theta} = \ frac {z_q} {r} $.
Cette intégrale est appelée ** Intégrale de diffraction de Fresnel-Kirchhoff **. Puisque l'amplitude complexe au point $ Q $ peut être calculée en l'intégrant numériquement, l'intensité sur le plan image peut être obtenue en calculant le carré de la valeur absolue. Cependant, cette méthode de calcul n'est pas pratique car le montant du calcul est $ O (N ^ 4) $ si le nombre de divisions entre la surface d'ouverture et la surface de l'image est $ N $.
Considérons une situation où $ \ theta $ est petit et $ \ cos {\ theta} \ approx 1 $ (approximation de l'axe proche).
En se concentrant sur l'équation 1, la partie $ \ cos {\ theta} $ devient 1, et $ \ frac {1} {r} $ peut être approximé comme $ \ frac {1} {z} $. Par contre, le terme $ e ^ {ikr} $ a une grande valeur de $ k = \ frac {2 \ pi} {\ lambda} $, donc si vous remplacez $ r $ par $ z $, la précision de l'approximation sera bonne. Il n'y en a pas.
Par conséquent, envisagez de développer $ r $ en une série et de l'utiliser jusqu'au deuxième terme. Quand $ r = \ sqrt {(x_q --x) ^ 2 + (y_q --y) ^ 2 + z_q ^ 2} $ est étendu au deuxième terme
En le substituant dans l'équation pour l'intégration de diffraction
Ce sera. Cette approximation est appelée ** approximation de Frenel **, et l'équation 2 est appelée ** intégration de diffraction de Frenel **.
Expansion du contenu de l'épaule $ e $ dans l'équation 2
Ce sera. Dans l'équation 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}
On voit donc qu'elle peut être calculée par la transformée de Fourier bidimensionnelle de $ u (x, y, 0) e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} $.
Pour que l'approximation de Fresnel soit maintenue, $ r $ multiplié par le troisième terme, qui est une extension de série, et $ k $ est nécessaire.
Doit être rencontré et réécrit pour $ z_q $
Ce sera.
Dans l'équation 3, $ \ frac {x ^ 2 + y ^ 2} {2z_q} $ est suffisamment petit et se rapproche de $ e ^ {ik \ frac {x ^ 2 + y ^ 2} {2z_q}} \ approx 1 $ Puis
Ce sera. Cette approximation est appelée ** approximation de Fraunhofer **, et l'équation 4 est appelée ** intégrale de diffraction de Fraunhofer **.
Dans l'équation 4, soit $ 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}
Par conséquent, nous pouvons voir qu'il peut être calculé par la transformée de Fourier bidimensionnelle de $ u (x_p, y_p, 0)
L'approximation de Fraunhofer est une approximation valide lorsque l'ouverture est suffisamment petite et que l'image est observée très loin de l'ouverture. D'autre part, l'approximation de Fresnel est une approximation valide encore plus proche que l'approximation de Fraunhofer.
La transformée de Fourier inverse de $ U (f_x, f_y, 0) $, qui est la transformée de Fourier de l'amplitude complexe $ u (x, y, 0) $ de la surface ouverte, est
est. Dans cette équation, l'amplitude complexe $ u (x, y, 0) $ de la surface d'ouverture est la somme de l'onde plane $ U (f_x, f_y, 0) \ exp [2 \ pi i (x f_x + y f_y)] $. Cela montre qu'il peut être exprimé comme.
Recherche de chaque composante du vecteur numérique d'onde $ \ boldsymbol {k} $ de chaque onde plane
\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}
Ce sera.
Lorsque chaque onde plane se déplace $ z_q $ dans la direction $ z $, son amplitude complexe
Ce sera. Alors, l'amplitude complexe $ u (x ', y', z_q) $ sur le plan image est la somme de ces ondes planes.
\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}
On voit que l'amplitude complexe $ u (x ', y, z_q) $ sur le plan image peut être calculée par la transformation de Fourier inverse de $ U (f_x, f_y, 0) \ exp [i k_z z_q] $. Cette méthode de calcul est appelée ** Méthode du spectre angulaire **.
La méthode du spectre angulaire peut calculer la diffraction dans la région plus proche de la surface d'ouverture que l'approximation de Frenel ou l'approximation de Fraunhofer. De plus, comme le calcul n'est effectué qu'une seule fois pour la transformée de Fourier et la transformée de Fourier inverse, il peut être calculé à grande vitesse.
Ouverture $ \ Sigma $ avec un rayon $ r = 1 ~ \ mathrm {[mm]} $ ouverture circulaire, $ \ lambda = 550 ~ \ mathrm {[nm]} $, $ z_q = 5 ~ \ mathrm {[m] Calculons numériquement comme} $. Le bloc-notes Jupyter utilisé pour le calcul se trouve sur Google Colab.
Tout d'abord, calculons numériquement avec l'approximation de Fraunhofer. Si vous vérifiez les conditions d'approximation, $ x $ et $ y $ seront $ r $ au maximum.
Ce sera. Il ne répond pas aux hypothèses, mais je vais procéder au calcul tel quel.
La diffraction de Fraunhofer peut être calculée en effectuant une transformée de Fourier bidimensionnelle sur l'amplitude complexe $ u (x, y, 0) $ sur l'ouverture $ \ Sigma $. Pour effectuer des calculs numériques par programme, divisez la zone de calcul avec une longueur de côté de $ D $ en une grille de $ N \ fois N $ afin de couvrir la surface d'ouverture, et stockez les amplitudes complexes sur chaque grille. Peut être calculé en effectuant une transformée de Fourier discrète bidimensionnelle.
Calculons numériquement en utilisant Python. Tout d'abord, préparez un tableau qui stocke des amplitudes complexes sur la grille.
import numpy as np
import matplotlib.pyplot as plt
l = 0.55e-6 #longueur d'onde[m]
k = 2 * np.pi / l #Nombre de vagues
zq = 5 #Distance entre la surface d'ouverture et la surface de l'image[m]
N = 200 #Nombre de divisions de grille
r = 0.001 #Rayon d'ouverture[m]
D = 0.05 #Zone de calcul sur la surface ouverte[m]
fscale = l*zq
#Renvoie l'amplitude complexe sur la surface ouverte
def P(x, y):
if x**2 + y**2 < r**2:
return 1
else:
return 0
#Génération de maillage
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)
Tracons-le.
plt.imshow(Z, extent=[-D/2, D/2, -D/2, D/2], cmap='gray')
plt.xlabel('$x_p$ [m]')
plt.ylabel('$y_p$ [m]')
La partie blanche correspond à la surface d'ouverture.
Ensuite, celle-ci est soumise à une transformation de Fourier discrète, mais avant cela, la région sur le plan image correspondant à la région de calcul est calculée.
v = np.fft.fftfreq(N, d=D/N)
v = fscale * np.fft.fftshift(v)
#Unités[mm]
extent= [1e3 * v[0], 1e3 * v[-1], 1e3 * v[0], 1e3 * v[-1]]
L'image de diffraction de Fraunhofer est obtenue par transformée de Fourier discrète d'amplitude complexe sur la surface ouverte.
#Calculer l'amplitude complexe sur le plan de l'image
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)
#Calculer la force
I_fraun = np.abs(U_fraun)**2
#Normalisation
I_fraun = I_fraun / np.max(I_fraun)
Tracons l'intensité.
plt.imshow(I_fraun, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fraunhofer Diffraction')
Il y a aussi un pic d'intensité autour du centre, ce qui montre que l'effet de diffraction peut être calculé.
Vérifiez d'abord les conditions d'approximation. L'expression conditionnelle pour l'approximation de la diffraction de Fresnel est
était. C'est un peu compliqué car la taille de la surface d'observation doit être prise en compte, contrairement au cas de la diffraction de Fraunhofer. Pour le moment, la taille d'un côté de la surface d'observation
Ce sera.
La diffraction de Frennell peut être calculée en utilisant la transformée de Fourier discrète bidimensionnelle ainsi que la diffraction de Fraunhofer.
#Calculer l'amplitude complexe sur le plan de l'image
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)
#Calculer la force
I_fresnel = np.abs(U_fresnel)**2
#Normalisation
I_fresnel = I_fresnel / np.max(I_fresnel)
Tracons l'intensité.
plt.imshow(I_fresnel, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Diffraction')
Il semble que le pic autour du centre soit plus fort que dans l'approximation de Fraunhofer.
La formule de diffraction de Fresnel Kirchhoff
Calculons en intégrant directement numériquement. Cette méthode est relativement rigoureuse car elle n'utilise pas d'approximations, mais elle n'est pas pratique en raison de l'énorme quantité de calculs. Je voulais le comparer avec la méthode d'approximation, alors je l'ai calculé.
#Intégration numérique diffractive de Fresnel-Kirchhoff
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()
#Zone de calcul sur le plan image
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))
#Intégration numérique(La mise à l'échelle est appropriée)
U_strict = np.vectorize(Uf)(X_Q, Y_Q)
#Calculer la force
I_strict = np.abs(U_strict)**2
#Normalisation
I_strict = I_strict / np.max(I_strict)
Tracons l'intensité.
plt.imshow(I_strict, cmap='jet', extent=extent)
plt.xlabel('$x_q$ [mm]')
plt.ylabel('$y_q$ [mm]')
plt.title('Fresnel Kirchhoff Diffraction')
Enfin, effectuons un calcul numérique en utilisant la méthode du spectre angulaire. La procédure consiste à trouver d'abord la transformée de Fourier $ U (f_x, f_y, 0) $ de l'amplitude complexe $ u (x, y, 0) $ de la surface ouverte, puis $ U (f_x, f_y, 0) $ et $ . Inverse Fourier transforme le produit de exp [ik_z z_q] $.
#Transforme de Fourier l'amplitude complexe de la surface ouverte
U = np.fft.fft2(Z)
#Calcul du vecteur de nombre d'onde
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)
#Calculer l'amplitude complexe du plan image
U_angular = np.fft.ifft2(U * np.exp(1.0j * k_z * zq))
#Calculer la force
I_angular = np.abs(U_angular)**2
#Normalisation
I_angular = I_angular / np.max(I_angular)
Tracons-le.
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')
Contrairement au cas de l'approximation de Fresnel et de l'approximation de Fraunhofer, la zone du plan image correspond à la zone de calcul, elle semble donc petite.
Cela fait longtemps, alors j'écrirai le reste sur une autre page.