はじめに
散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので波動光学に基づいてPythonで計算してみたシリーズの第2回です。前回FFTを使って波動光学的に結像計算をしましたが、ピラミッドみたいな形になってしまいました。
今回は回折や収差による像ボケをFFT結像に取り込みながら、なぜそうなったのかを明らかにします。
頻出のPSFやゼルニケ多項式についても腹落ちできるように説明します。
前回記事
回折や収差の取り込み方
瞳面の電場をFFTすると像面での電場分布が得られます。
この関係を線形システムとして一般化すると、瞳面上の電場分布と振幅分布関数(以下ASF)の畳み込み計算により、像面での電場分布が計算できます。
このときASFの絶対値の2乗がよく知られる点像分布関数(PSF)です。
ASFは光学系におけるの電場振幅のインパルス応答と言えます。
ASFは瞳関数をFFTすれば求められます。
瞳関数は瞳面形状に位相を取り込んだものです。
位相の項により収差を表現でき、ゼルニケ関数により収差と位相が結び付きます。
実際に1個ずつ実装していきます。
回折効果の計算
まずは瞳面を設定します。
前回記事で説明したように瞳面は対物レンズの焦点距離単位で規格化されていたのでした。
ここでは、最初はピクセル数で設定して後で物理量を与えます。
import numpy as np
import matplotlib.pyplot as plt
import math
N = 512
# 座標（中心を0にする）
x = np.arange(N) - N//2
y = np.arange(N) - N//2
X, Y = np.meshgrid(x, y)
# 瞳面の半径(ピクセル単位)
R = N // 4
# 円マスク
mask = (X**2 + Y**2) <= R**2
# 0/1配列
arr = mask.astype(float)
plt.figure(figsize=(6,5))
plt.contourf(X,Y,arr,levels=50)
plt.xlabel("X")
plt.ylabel("Y")
plt.gca().set_aspect('equal')
plt.show()瞳面は以下のように円形です。
内側は1で光を通し、外側は0で光を通しません。
これは単純な瞳関数ですが、複素数に拡張するこで収差を表現できます。
瞳関数をFFTしてASFとPSFを求めます。
ASF = np.fft.fft2(np.fft.fftshift(arr))
# FFTは左上を原点として扱うため、
# 中心を原点とした瞳関数をfftshiftで並べ替える
ASF = np.fft.fftshift(ASF)
PSF = np.abs(ASF) ** 2実スケールに変換します。
瞳面の半径はレンズの焦点距離にNAを掛ければ求まります。
半径はRピクセルで瞳面を表わしています。
#光学系の倍率は1
f_obj = 100#f :対物レンズ焦点距離
f_img =100#f :結像レンズ焦点距離
NA = 0.5
lam = 0.532 *10**-3 #: 波長
#瞳面での1ピクセルの物理的大きさ
dx = f_obj * NA / R
k = np.fft.fftfreq(N, d=dx)
x_real = lam * f_img * k
x_real = np.fft.fftshift(x_real)
y_real = np.copy(x_real)
X_real,Y_real=np.meshgrid(x_real,y_real)
#PSFをグラフ化
plt.figure(figsize=(6,5))
pix =50
plt.contourf(X_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],Y_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],PSF,levels=50)
plt.xlabel("X")
plt.ylabel("Y")
plt.gca().set_aspect('equal')
plt.show()PSFが大体2um程度に結像出来ていて、エアリーディスクに対応する回折限界を表します。
これで回折による像のボケを表せました。
PSFは理想的な点光源を集光した時の形状を表しています。
そしてPSFよりも小さく集光することは基本的にはできません。
光学系収差の取り込み方
次に収差を取り込みます。
ここまで瞳関数は実数でしたが、位相ファクタexp(i2π*W)により収差を表現できます。
ここでWは波長単位での位相ずれを表現していて、ゼルニケ多項式を使うのが便利です。
自作でゼルニケ多項式を組みます。
#ゼルニケ関数の定義
#ゼルニケ関数の半径成分
def zernike_radial(n, m, rho):
m = abs(m)
if (n - m) % 2 != 0:
return np.zeros_like(rho)
R = np.zeros_like(rho, dtype=float)
for k in range((n - m)//2 + 1):
c = ((-1)**k * math.factorial(n - k) /
(math.factorial(k) *
math.factorial((n + m)//2 - k) *
math.factorial((n - m)//2 - k)))
R += c * rho**(n - 2*k)
return R
def zernike(n, m, rho, theta):
R = zernike_radial(n, m, rho)
if m > 0:
return R * np.cos(m*theta)
elif m < 0:
return R * np.sin((-m)*theta)
else:
return Rゼルニケ多項式はnとmの値を指定すれば、収差と対応します。
nは自然数、mは整数で$-n \leq m \leq n$です。
例えば(n,m)=(2,0)とデフォーカスに対応する等、扱いやすいです。
ゼルニケ多項式の形はWikipediaの図が分かりやすいです。ゼルニケ多項式（Wikipedia）
#規格化瞳面半径
rho = np.sqrt(X**2+ Y**2) / R
#瞳面の角度
theta = np.arctan2(Y,X)
# (n,m) 指定で適当に足す(係数はwaves単位)
phi_waves = (
0.25 * zernike(2, 0, rho, theta) + # defocus
0.15 * zernike(2, 2, rho, theta) + #astigmatism
0.2 * zernike(2, 1, rho, theta) #tilt
)
phi = 2*np.pi * phi_waves
arr = mask.astype(float)
#位相ずれ成分を瞳関数にかける
arr = arr * np.exp(1j*phi)
#FFTでASFとPSFを求める
ASF = np.fft.fft2(np.fft.ifftshift(arr))
ASF = np.fft.fftshift(ASF)
PSF = np.abs(ASF)**2
#グラフ化
plt.figure(figsize=(6,5))
pix =20
plt.contourf(X_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],Y_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],PSF[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],levels=50)
plt.xlabel("X (mm)")
plt.ylabel("Y (mm)")
plt.gca().set_aspect('equal')
plt.show()収差を入れたことでPSFが左右に広がって見えます。
収差の入れ方を少しずつ変えてPSFを見ていくと収差の特徴が見えてきます。
位相の分布も確認してみましょう。
plt.figure(figsize=(6,5))
plt.contourf(X,Y,phi*mask,levels=50)
plt.xlabel("X (pixel)")
plt.ylabel("Y (pixel)")
plt.gca().set_aspect('equal')
plt.show()x方向の位相の振動が大きいことが分かります。
PSFもx方向に広がっていたので対応が見えます。
ピラミッドみたいになった理由
今回設定ではNA=0.5としたため瞳関数は512 x 512の約半分になっていますが、実はこれがうまく計算出来た理由です。
FFTは周期境界条件を仮定するため、瞳を配列いっぱいに広げると折り返しの影響がでます。
前回記事でピラミッドみたいに見えたのはこれが原因でした。
NA=0.9として計算したため瞳関数が大きくなりすぎて折り返しの影響が見えていたのです。
もしNA=0.9で計算したいときは瞳面の座標の取り方を工夫する必要があったのです。
参考文献
フーリエ光学の理論面についてさらに詳しく学びたい方は、Goodmanの『Introduction to Fourier Optics』がおすすめです。フーリエ光学の定番教科書で、今回扱った内容の背景も体系的に学べます。邦訳もでているのは有難い。
本文中でフーリエ変換が線形システムと述べましたが、畳み込み積分とフーリエ積分の数学的な関係の詳細について押さえたい方は谷田貝先生の光とフーリエ変換をおすすめします。
次回
次回は結像計算から少し離れて、散乱光の計算方法について解説してみたいと思います。
コメント