FFT結像に回折と収差を取り込む 〜PSFと瞳関数の話〜

導入

散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので計算してみたシリーズの第4回。

今回はついに回折による像のボケについて書きたいと思います。
PSF(Point Spread Function)やゼルニケ多項式を何となく使っている人向けです。

点光源からの光でも結像すると、回折の効果で有限の大きさに結像されます。
またレンズ含めた光学系の収差に由来してボケが発生します。

これらの効果をこれまでのFFT結像に取り込んでいきます。

FFTによる結像計算

結像計算したけど、像の大きさってどう決まるの？

散乱光の偏光のS/Pってどう決まるの？

FFT結像で回折や収差の効果の取り込み方

瞳面での電場をFFTすると像面での電場分布が得られます。
この関係を線形システムとして一般化すると、瞳面上の電場分布と振幅分布関数(以下ASF)の畳み込み計算により、像面での電場分布が計算できます。

注:ここでは照明光はレーザ光としていて、コヒーレント光学系と呼ばれます。

このときASFの絶対値の2乗がよく知られるPSFです。
ASFは光学系におけるの電場振幅のインパルス応答と言えます。

ASFが分かれば回折や収差の効果を取り込むことができます。
ASFは瞳関数をFFTすれば求められます。
瞳関数は瞳面形状に位相を取り込んだものです。
位相の項により収差を表現でき、ゼルニケ関数により収差と位相が結び付きます。

次節から実際に1個ずつ実装していきます。

理論をちゃんと抑えたい人はJoseph GoodmanのIntroduction to Fourier Opticsを参照ください。
この辺りをちゃんと書いている和書は殆どありません。

回折効果の計算

まずは瞳面を設定します。
瞳面の物理的な大きさを決める必要がありますが、後で考えます。

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.ifftshift(arr))
# FFTは左上を原点として扱うため、
# 中心を原点とした瞳関数をifftshiftで並べ替える
ASF = np.fft.fftshift(ASF)
PSF = np.abs(ASF) ** 2

実スケールに変換します。
瞳面の半径はレンズの焦点距離にNAを掛ければ求まります。
半径はRピクセルで瞳面を表わしています。

# dx : 瞳面のサンプリング間隔の実サイズ
#lam: 波長
#f :焦点距離
#単位はmm
f = 100
NA = 0.5
lam = 0.532 *10**-3

#瞳面での1ピクセルの物理的大きさ
dx = f * NA / R

k = np.fft.fftfreq(N, d=dx)
x_real  = lam * f * 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は以下のように計算出来て、縦軸と横軸の単位はmmです。
大体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<=m<=nです。
例えば(n,m)=(2,0)はデフォーカスに対応する、などです。

wikipediaの図は一度ご覧ください。
https://ja.wikipedia.org/wiki/%E3%82%BC%E3%83%AB%E3%83%8B%E3%82%B1%E5%A4%9A%E9%A0%85%E5%BC%8F

        
#規格化瞳面半径
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

適当に収差を組み合わせてPSFを計算しました。
可視化してみましょう。

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方向に広がっているので対応しています。

瞳関数の半径(ピクセル単位)について

瞳関数は512 x 512の内接円にはせずちょっと小さくしています。
小さくすると0のピクセルが増えるだけで、情報増えないのに無駄な気がします。

# 瞳面の半径(ピクセル単位)
R = N // 4
# R = N // 2とするとNG

実際に内接円にしてみると分かりますが、PSFが十分にサンプリングされず円でなく四角形のように見えてしまいました。
逆に小さくしすぎると有効なピクセル数が減って特徴がうまく現れてきません。

しれっと4で割っていますが、試行錯誤があります。

次回

前前回に計算した瞳面の強度分布に今回の瞳関数をかけてFFTかければ散乱光の結像計算完了!!と言いたかったのですが、一番重い計算が残っています。

検出NAが大きい条件での偏光の取り扱いです。
レンズによって光は各点でそれぞれに曲げられ、電場の振動方向が変わります。
これを直交する方向に分解して足しあげる必要があるのです。

• 次回記事(?)