【結像計算:第1回】結像サイズは何で決まるのか？ ― フーリエ光学をPythonで実装して理解する

Python
2026.05.30

はじめに

よく散乱光分布の計算をするのですが、ふとこれを結像したらどう見えるんだろう？と思ったので波動光学に基づいてPythonで計算してみました。

実際にやってみるとハマりどこがとても多い。
特にFFTをすると、実空間から周波数空間に移るため像大きさがそもそも何の単位なのか分からなくなります。

この記事を読めば、結像計算における“像のスケールの正体”を、フーリエ光学と数値実装の両面から整理することができます。

問題設定

散乱光分布の計算では、各(θ,φ)方向への散乱光強度が求まります。
これを波動光学的にフーリエ変換して結像してみたいと思います。(フラウンフォーファー回折)
方針としては
(θ,φ)方向情報を瞳面での分布に変換、フーリエ変換すれば良さそうです。

計算フロー

詳細はこうなります。

  1. 散乱光のダミーデータを準備(方向空間)　
  2. 瞳面座標に変換(k空間)
  3. 等間隔グリッドに再配置(FFTのため)
  4. 強度補正
  5. FFTして結像
  6. スケール計算

1. 散乱光のダミーデータを準備(方向空間)

散乱光電場の計算結果のダミーデータを作ります。
散乱光計算では、(θ,φ)に対して等間隔で散乱光強度を計算します。

import numpy as np
from scipy.interpolate import griddata #補間用
import matplotlib.pyplot as plt
plt.style.use("ggplot")

#度から弧度法への変換係数
deg = np.pi / 180

#dummyデータ生成
theta = np.linspace(0,90*deg,90)
phi = np.linspace(0,360*deg,360)
theta,phi = np.meshgrid(theta,phi)
E = np.cos(theta)

#コンター図作成
theta_cont = 90 - theta / deg #散乱角から仰角に。コンター図用に度に直す
phi_cont = phi / deg#コンター図用に度に直す
plt.figure(figsize=(6,5))
plt.contourf(phi_cont,theta_cont,E,levels=50)
plt.xlabel("azimuth angle(deg)")
plt.ylabel("elevation angle(deg)")
plt.show()

上方への散乱光強度を強くなるようにしました。
thetaは散乱角なので直上方向が基準ですが、コンター図にする際は90度から引いて
仰角に変換すればメルカトル図法っぽくなり直感的に理解しやすいです。
theta_phi_contour.jpg

2. 瞳面座標に変換(k空間)

theta,phiからk空間に座標を変換します。

#theta,phiからk空間に変換
kx = np.sin(theta)*np.cos(phi)
ky = np.sin(theta)*np.sin(phi)

plt.figure(figsize=(6,5))
plt.contourf(kx,ky,E,levels=50)
plt.xlabel("kx")
plt.ylabel("ky")
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.gca().set_aspect('equal')
plt.show()

電場は何も変えていないのですが、
座標をk空間に変換したことで輝点を真上から見たようなコンター図になります。
pupil_dist.jpg

ここで重要なのは、
k空間の座標がそのまま「瞳面の座標」に対応するという点です。

3. 等間隔グリッドに再配置(FFTのため)

(θ,φ)で等間隔にサンプリングしたデータをそのままk空間に変換すると、
瞳面上では点の密度が不均一になります。

しかしFFTは等間隔サンプリングを前提としているため、

一度、等間隔グリッドに再配置する必要があります。

# FFTするデータ数、2の累乗にする。
N=512 
#kを等間隔で作成
kx_im = np.linspace(-1,1,N)
ky_im = np.linspace(-1,1,N)
KX_im,KY_im =np.meshgrid(kx_im,ky_im)

#電場強度を補間
points = np.column_stack([kx.ravel(),ky.ravel()])
values = E.ravel()

#後でtheta使うので補間
values_theta=theta.ravel()

E_interpolate = griddata(points,values,(KX_im,KY_im),method='linear')
theta_interpolate = griddata(points,values_theta,(KX_im,KY_im),method='linear')

#nanが含まれるので0にする
E_interpolate[np.isnan(E_interpolate)]=0
theta_interpolate[np.isnan(theta_interpolate)]=0

E_interpolateが512×512の等間隔グリッドに補間したものです。
thetaも次で使うので同じように等間隔グリッドに補間しています。

4. 強度補正

以下の2点を考慮する必要があります。

  • レンズで拾える範囲に電場を制限する必要があります。
    ここでは仮にNA = 0.9以下が検出可能、それ以上検出できないので0とします。
  • レンズが適切に収差補正されている(アプラナート条件)場合は、$\sqrt{\cos\theta}$を掛ける必要があります。検出NAが小さいときは$\cos\theta$は1に近似できるので省略してもあまり影響がありません。
NA=0.9
mask = (np.sqrt(KX_im**2+KY_im**2) < NA).astype(int)
weight = np.cos(theta_interpolate)**0.5
E_interpolate *= mask*weight

図示してみるとこんな感じで、ほとんど変化ありませんね。
image_comp2.jpg

5. FFTして結像

後はFFTするだけです。
瞳面の電場分布と像面の電場はフーリエ変換の関係にあります。

E_img = np.fft.fft2(E_interpolate)
#FFT結果を並べなおす
E_img = np.fft.fftshift(E_img)

nx=ny= N
x = np.fft.fftfreq(nx, d = 1/N)
y = np.fft.fftfreq(ny, d = 1/N)
x = np.fft.fftshift(x)
y = np.fft.fftshift(y)
X, Y = np.meshgrid(y, x)

#図示用のコード
plt.figure(figsize=(6,5))
plt.contourf(X,Y,np.abs(E_img)**2,levels=50)
plt.xlabel("X")
plt.ylabel("Y")
plt.gca().set_aspect('equal')
plt.show()

散乱光はある一点から発生していて、それを結像したので小さな1点に結像されています。
しかし、これはピクセル単位なので実際の大きさが分かりません。
更に、FFTしたので軸の単位は実空間$[\mathrm{m}]$から周波数$[\mathrm{m^{-1}}]$に移っているはずです。
image.jpg

6. 像面スケール

像面での座標xと瞳面上での周波数k、波長:λ、レンズの焦点距離をfとすると、変換式は下記です。

$$
x = \lambda f k
$$

更に以下のように波長情報を消去できます。

$$x = \lambda\times f\times k$$

$$ 　　= \lambda\times f\times \frac{\sin\theta}{\lambda}$$

$$ = f\times\sin\theta$$

$$　 = f \times k_{\mathrm{norm}}$$

波数kは規格化されているため、瞳面上の空間座標はレンズ焦点距離でスケールされるのです。

そしてコードは、、、といきたいのですが次節で顕微鏡を例にもう少し考えてみます。

顕微鏡の話

FFT結像で見たスケールの話は、実際の光学系でもそのまま出てきます。
その分かりやすい例が顕微鏡です。
顕微鏡の構成はこうです。

サンプル – 対物レンズ -　チューブレンズ

実際はチューブレンズの後に接眼レンズがあり、人の瞳にあうように平行光にするのですが省略します。
チューブレンズは固定で焦点距離が決まっています。
日本の主なメーカでは200mm、Olympusだけ180mmとしているみたいです。

対物レンズの焦点距離を$f_{obj}$、チューブレンズの焦点距離を$f_{tube}$とすると
倍率Mは以下と書けます。
$$
M = \frac{f_{\mathrm{tube}}}{f_{\mathrm{obj}}}
$$

対物レンズとチューブレンズの焦点距離の比で光学系の倍率が決まり、
顕微鏡はリボルバーで対物レンズを変えることで倍率を変えているのです。

7. スケール計算

顕微鏡のように、FFT結像する際にも対物レンズとチューブレンズのように検出と結像でのレンズ焦点距離を分けて考えておくと、理解度があがります。
瞳面のNAは<0.9など制限がありますが、これは対物レンズで決まります。

つまり、瞳面上の座標は対物レンズの焦点距離単位なのです。

FFTにより、瞳面での電場分布から像面の分布を求めるのは結像レンズの働きなので座標変換の係数に用いるレンズ焦点距離にはチューブレンズ(結像レンズ)の値を用いればいいのです。

対物レンズと結像レンズの焦点距離が20mmと200mmで、倍率10として計算してみます。

#像面のスケーリング
f_obj = 20 #対物レンズ焦点距離
f_img = 200 #結像レンズ焦点距離
lam = 0.532 *10**-3#波長

#瞳面での1ピクセルの物理的大きさを求める
dx = f_obj / (N/2)

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)

#グラフ化
plt.figure(figsize=(6,5))
pix =5

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],np.abs(E_img[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2])**2,levels=50)
plt.xlabel("X(mm)")
plt.ylabel("Y(mm)")
plt.gca().set_aspect('equal')
plt.show()

縦軸を長さ単位で計算することが出来ました。
倍率10倍で10μmくらいになっているので、妥当な大きさです。

scaled_img.png

ただし、像がピラミッドみたいな形になってしまいました。

実はこれは瞳面のNAの大きさの設定に起因するものです。記事も長くなってしまったので、次回以降で原因の詳細と対処方法を説明したいと思います。

参考文献

フーリエ光学の理論面は下記を参照ください。和書だとここまで詳しく書かれているものはほとんどありません。

次回

　次回以降でPSFや波動光学での収差の扱いにも触れながらどうすれば自然に結像を計算できるかを解説したいと思います。

次回記事

【結像計算:第2回】なぜ結像結果はピラミッド形状になったのか？PSFとゼルニケ多項式で原因を探る
はじめに　散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので波動光学に基づいてPythonで計算してみたシリーズの第2回です。前回FFTを使って波動光学的に結像計算をしましたが、ピラミッドみたいな...
yohaku-labs.com
2026.05.31

コメント

タイトルとURLをコピーしました