光学設計とその周辺、そしてたまに全く関係ないやつ

お久しぶりの記事ですが, 今回ご紹介したいのはおなじみMatlabに光学設計用のアドオンであるOptical System Design and Analysisがリリースされたというニュースです.
jp.mathworks.com

Matlabから光学設計ツールを制御する機能は従来からありましたが, 今回はMatlab自体に機能が追加されたということでインパクトのあるニュースだと思います.

操作方法としてまずコーディングベースの操作とこのアドオン用のGUIアプリであるOptical System Designer Appによる操作の2つの方法があります. それぞれで何ができるかは以下のページにまとめられていますが, 基本はコーディングで全ての操作ができ, 一部の機能はアプリの方ではできないというものです.
Get Started with Optical Design and Simulation - MATLAB & Simulink

GUIアプリの方ではどんな雰囲気になっているかは以下のような感じ. ちょっと見づらいかもしれません.

一番最初に紹介したリンクから各設定の操作方法がまとまっています. 機能としては2Dのレンズ図や3Dの可視化は可能, 収差曲線やスポットダイアグラムは使える一方, MTFとかは未実装のよう. 最適化は一応あるようですが, 公差解析などはなさそうです. また, レンズ(ミラー)設計自体に加えコーティングの設定とその設計自体も機能としてあるようです. トリプレットの設計例やズームレンズの設計事例も紹介されています.

厳密にいうとMatlabの既存のアドオンである画像処理用のImage Processing Toolboxにこの光学設計用のライブラリが追加されたというもので, 使用にするにはMatlab本体とこのImage Processing Toolbox, そしてもし最適化も利用したい場合はさらにOptimization Toolboxも必要になるようです. このImage Processing Toolbox内に組み込まれているように, 従来このToolboxには光学設計とかなり近い機能(PSFからOTFを求めるpsf2otfなど)が既にありました. 恐らくそういうところからユーザーからの要望が挙がるなどして実現した機能だと推測します.

コスパ的なことを言うと, 永久ライセンスでMatlab本体が35万ぐらい, Image Processing Toolboxが20万くらいとすると既存のCodeVやOpticstudioよりははるかに手ごろな価格になります. 価格に加え, Matlabの従来の機能である拡張的な関数や機能を利用し複雑な解析をするとか, GPUとかクラウドコンピューティングなどのMatlabプラットフォームも利用できる可能性もあるなど他のCADソフトにないメリットがあるのかも.

機能としてはまだまだ足りないところも多いですが, もし標準的な機能が追加されていけば実際の現場でも使われてくる可能性もあるかもしれません.

前回の記事では自由空間の回折による光波の伝搬を扱いました. 2回目の今回は実際にレンズがあった場合にどのように像面の強度パターンが回折の影響を受けるか扱います.
波動光学的にPSF, MTFをシミュレーションしてみる1(1/4) -自由空間の回折伝搬- - 光学設計とその周辺、そしてたまに全く関係ないやつ

光学系のモデル化

カメラレンズのように何枚も光学素子があるような光学系の中の回折伝搬を考える際, もちろん一つ一つの素子で個別に伝搬を計算するのが発想としてはシンプルです. 一方で, モデル化をして, 光学系自体を単純化することができれば手間という意味ではシンプルで, 計算量も少なく済むかもしれません.
そのモデルとして用いられるのが, 以下のような入射・射出瞳です. そもそも瞳とは何かは他文献を見てもらうとして, 射出瞳を用いることで, 複雑な光学系も単一のレンズ系として扱うことができ, 最初に射出瞳位置やそのサイズを光線追跡で求める作業は必要ですが, 全ての回折伝搬計算は射出瞳位置での1回だけで行えばよくなります. また理論的にも見通しが良くなり, いわゆる瞳の相関がOTFだとかそういう重要な考察が得られるわけです. この後者の手法をCodeVの開発元は射出瞳回折モデルと呼んでいたような. 本記事では射出瞳モデルと呼びます.

さて2種類のアプローチを紹介しましたが, 世にある光学設計ソフトではどちらの機能も実装されていることが多いです. ただ, カメラレンズとか顕微鏡のような一般的なアプリケーションでPSFやMTFを計算する目的なら, 射出瞳モデルを使うことが多いはず. 設計ソフトにあるMTF解析機能も基本的にこれに沿っています. 逆に, 特殊な状況*1では各面で個別に回折伝搬を計算することもあり, 例えばCodeVのビーム伝搬解析とかZemax OSの物理光学伝搬などの機能がこちらに相当します. 本記事の目的はMTFを得ることなので, 射出瞳モデルで考えていきます.

射出瞳モデルを使った光の伝搬とレンズのフーリエ変換作用

射出瞳モデルを使うことで, 物体からの光が単一のレンズで集光されて像面に光が伝搬されますが, このフォーカシングによる光場への影響を見てみます. 改めて上の図を見ると, 物体からの発散した球面波が瞳に入射し, レンズが像点に新たな収束する球面波を作ります. 以降レンズのある座標を $(\xi,\eta)$ で表します.
この物体-瞳の球面波は物体から瞳の距離をaとすると瞳直前での波面は以下のように表されます.

$\displaystyle U_1(\xi,\eta) = A \exp{ \left[ i \frac{\pi}{\lambda a} (\xi^2+\eta^2) \right] } \tag{1-1} \label{1-1}$

また瞳-像を伝搬する光はこの距離をbとすると瞳直後の波面は以下のように同じ形の式となります*2.

$\displaystyle U_2(\xi,\eta)= A' \exp{ \left[ -i \frac{\pi}{\lambda a} (\xi^2+\eta^2) \right] } \tag{1-2} \label{1-2}$

A->A'と振幅が変化したのはレンズによる透過率などを考慮したことによります.
この時, レンズによる波面の作用を $t(x,y)$ とすると, $U_2(\xi,\eta)=t(\xi,\eta) U_1(\xi,\eta)$ となるため, t(x,y)は係数比を $T_0$ とまとめると

$\displaystyle t(\xi,\eta)= \frac{U_2(\xi,\eta)}{U_1(\xi,\eta)}=T_0 \exp{\left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{1-3} \label{1-3}$

ここで結像の式 $1/a+1/b=1/f$ を使いました. ここからわかる通り, レンズによる位相への影響は物体像配置にかかわらず, そのレンズの固有パラメーターである焦点距離のみで表すことが出来ます.

レンズの有限の開口を表す関数を瞳関数と呼びますがこれを $p(\xi,\eta)$ とすると, この関数も式(1-3)にそのまま加えることが出来ます. また収差の影響もこの $t(\xi,\eta)$ に含めることができ, 以下の図のように収差というのは射出瞳位置における理想球面の波面と実際の収差がある波面の光路差となるため, これを $W(\xi,\eta)$ とすると, 先の瞳関数も合わせて, 以下式(1-4)のようにまとめられます*3.

$\displaystyle t(\xi,\eta)= \frac{U_2(\xi,\eta)}{U_1(\xi,\eta)}=T_0 p(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \exp{ \left[ -i W(\xi,\eta)\right] }$
$\displaystyle = T_0 P(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{1-4} \label{1-4}$

となります.

点像分布関数 PSF(Point Spread Function)と周波数応答関数

以上の内容をベースに点像分布関数PSF, つまり点光源を結像させた際の像面での広がり具合, を見ていきます.
まず結像にはコヒーレント結像とインコヒーレント結像の2種類があります. 通常PSFとかMTFとかいうときはインコヒーレント結像を想定することが多いのですが, 進行の都合上コヒーレント結像をまず考えます. コヒーレント結像とは物体からの光波の振幅と位相がそのまま光学系を通って像面で干渉・回折し, 振幅の重ね合わせで像を形成する状況です. つまり像面のある点における振幅, 強度は物体面全体からの影響を受けるということです.
以下の図のような状況を考えていきますが, 物体面(座標 $xy_1$ )の位置( $x_o,y_o$ )にある点 $P_1$ から距離 $z_1$ の位置に焦点距離fの単レンズがあり, この点 $P_1$ から伝搬する光波の振幅と位相がそのまま光学系を通ってレンズから距離 $z_2$ を通って像面(座標 $xy_2$ )に干渉・回折し位置( $x_i,y_i$ )にある点 $P_2$ に像が形成されています.

重ね合わせ積分を使うと, 像面における振幅 $U_2(x_i,y_i)$ は, 物体面の振幅 $U_1(x_o,y_o)$ を使って

$\displaystyle U_2(x_i,y_i)= \iint_{\infty} U_1(x_o,y_o) h(x_i,y_i;x_o,y_o) dx_o dy_o \tag{2-1} \label{2-1}$

となります. ここでhはインパルス応答であり, 点像つまりPSFそのものです. ここからの目的はこのPSFであるhをこれまでの回折計算結果を使って理論的に求めていきます.

物体からの1点が作る振幅分布は $U_1=\delta(x_o, y_o)$ とデルタ関数を使うと, 像面での振幅分布 $U_2$ 自体が点像分布関数PSFとなります.

前回の記事含めこれまでの内容を振り返りP1からP2の伝搬をステップ毎にまとめると,
1)物体面のデルタ関数が距離z1の分フレネル伝搬し, レンズまで到達(U1->Ui)

>デルタ関数と前回記事の式(2-3)のインパルス応答を使い式(1-3)の畳み込み積分で計算

2)レンズ位置で式(1-4)で表される位相関数 $t(\xi,\eta)$ が付加(Ui->Ui')

>1)の振幅と式(1-4)のレンズ作用tの掛け算を計算

3)その光が再度距離z2の分像面まで到達(Ui'->U2)

>2)の振幅を入力として再度1)の計算

をすることになります.

まず1)の計算をすると, 物体面上の点P1が座標 $(x_o,y_o)$ にあるとすると, 到達面のレンズ座標 $(\xi,\eta)$ 上の振幅位相 $U_i(\xi,\eta)$ は, 物体強度分布を表す入射光波はデルタ関数を使って $U_1(x_1,y_1)=A \delta(x_1-x_o, y_1-y_o)$ となるため,

$\displaystyle U_i(\xi,\eta)= A \delta(x_1-x_o, y1-y_o) \ast \frac{\exp{(ikz_1)}}{i \lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi^2+\eta^2) \right] }$
$\displaystyle = \frac{A \exp{(ikz_1)}}{i\lambda z_1} \iint_{\Sigma} \delta(x_1-x_o, y_1-y_o) \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_1)^2+(\eta-y_1)^2 \right] } d\xi d\eta$
$\displaystyle = \frac{A \exp{(ikz_1)}}{i\lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right] } \tag{2-2} \label{2-2}$

と計算できます.

次に2)の計算ですが, これは単に式(1-4)を掛け算をするだけです. ということでレンズ直後の振幅位相分布 $U_i'(\xi,\eta)$ は

$\displaystyle U_i'(\xi,\eta)=t(\xi,\eta) U_i(\xi,\eta)$
$\displaystyle =T_0 \frac{A \exp{(ikz_1)}}{i\lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right]} P(\xi,\eta)\exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{2-3} \label{2-3}$

です.

最後に3)の $U_2(x_i,y_i)$ 計算ですが, 1)の計算と同様に,

$\displaystyle U_2(x_i,y_i)=U_i'(\xi,\eta) \ast \frac{\exp{(ikz_2)}}{i \lambda z_2} \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i^2+y_i^2) \right] }$
$\displaystyle = T_0 \frac{A \exp{ \{ik (z_1+z_2) \} } }{(i\lambda)^2 z_1 z_2} \iint_{\infty} P(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] }$
$\displaystyle \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right] } \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i-\xi)^2+(y_i-\eta)^2 \right] } d\xi d\eta$
$\displaystyle = T_0 \frac{A \exp{ \{ik (z_1+z_2) \} } }{(i\lambda)^2 z_1 z_2} \exp{ \left[ i \frac{\pi}{\lambda z_1} (x_o^2+y_o^2) \right] } \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i^2+y_i^2) \right] }$
$\displaystyle \iint_{\Sigma} P(\xi,\eta) \exp{ \left[ i \frac{\pi}{\lambda} \left( \frac{1}{z_1}+\frac{1}{z_2}-\frac{1}{f} \right) (\xi^2+\eta^2) \right] } \exp{ \left[-i \frac{2\pi}{\lambda} \left( \left( \frac{x_o}{z_1}+\frac{x_i}{z_2} \right) \xi+\left( \frac{y_o}{z_1}+\frac{y_i}{z_2} \right) \eta \right) \right] } d\xi d\eta \tag{2-4} \label{2-4}$

ここで結像の式 $1/z_1+1/z_2=1/f$ を使えば, 瞳座標の2乗の項はキャンセルされるため, 式(2-4)はフーリエ変換の形になる. これがレンズのフーリエ変換と呼ばれる所以です.
また通常最終的には振幅ではなくその絶対値2乗である強度が重要なため, 積分の前にある $\exp{ \{ik (z_1+z_2) \} }$ や $\exp{ \left[ (x_i^2+y_i^2) \right] }$ は無視できます. 注意が必要なのは物体座標に関係する項 $\exp{ \left[ (x_o^2+y_o^2) \right] }$ です. というのは式(2-1)の畳み込み積分の積分変数と同じ変数となっているため, 単純に無視することはできません. 良く説明されるのがMTFのよい光学系なら物体点と像点がほぼ1対1となるため, 事実上無視できるということです. そんなに単純な話では実際はないですが, インパルス応答の絶対値2乗が伝達関数となるインコヒーレント結像ではどうせキャンセルされるため, これ以上は深堀せず以降は無視して考えていきます*4.
以上の議論をまとめ, この光学系の倍率mが $m=-z_2/z_1$ となることも利用すると式(2-4)を記述しなおすと,

$\displaystyle U_2(x_i,y_i)=h(x_i,y_i;x_o,y_o)$
$\displaystyle = \frac{A'}{\lambda^2 z_1 z_2} \iint_{\infty} P(\xi,\eta) \exp{ \left[-i \frac{2\pi}{\lambda z_2} \left[ (x_i-mx_o) \xi+ (y_i-my_o) \eta \right] \right] } d\xi d\eta \tag{2-5} \label{2-5}$

つまり, インパルス応答は像座標 $(mx_o, my_o)$ を中心とした開口関数pのフラウンホーファー回折像となります.

PSFをもとめるだけならこの式(2-5)で終わりなのですが, OTFを求めるには式(2-1)が畳み込みの形で表す必要がありますので, もうちょい式(2-5)を変えていきましょう.
まず式(2-5)を理想的な結像状態で幾何光学的に表すとどうなるかを見てみます. $\xi' = \xi/\lambda z_2, \eta'=\eta/\lambda z_2$ とし式を変形すると

$\displaystyle h(x_i,y_i;x_o,y_o)= A' m \iint_{\infty} P(\lambda z_2 \xi', \lambda z_2 \eta') \exp{ \left[-i 2\pi \left[ (x_i-mx_o) \xi'+ (y_i-my_o) \eta' \right] \right] } d\xi' d\eta' \tag{2-6} \label{2-6}$

幾何光学かつ理想的な結像となると瞳関数は瞳座標すべてで1, $\lambda \to 0$ となるため, この時式(2-6)はデルタ関数の定義そのものとなるため,

$\displaystyle h(x_i,y_i;x_o,y_o) \to A' m \delta( x_i-mx_o, y_i-my_o) \tag{2-7} \label{2-7}$

となります.
これを式(2-1)に代入すると,

$\displaystyle U_{2g}(x_i,y_i) = \frac{1}{m}U_1(x_i/m,y_i/m) \tag{2-8} \label{2-8}$

です. 回折も考慮した振幅U_2とこの幾何光学的な振幅を区別するためこの式(2-8)はU_2gとしています.
$mx_o=x_o', my_o=y_o', h'=h/m$ とし式を変形すると式(2-8)も使うと,

$\displaystyle U_2(x_i,y_i)= A' \iint_{\infty} h'(x_i-x_o',y_i-y_o') \frac{1}{m} U_1\left( \frac{x_o'}{m},\frac{y_o'}{m} \right) dx_o' dy_o' = h'(x_o',y_o') \ast U_{2g}(x_o',y_o') \tag{2-9} \label{2-9}$

強度の場合は

$\displaystyle I_2(x_i,y_i) = | U_2(x_i,y_i)|^2= |h'(x_o',y_o') \ast U_{2g}(x_o',y_o') |^2 \tag{2-10} \label{2-10}$

です.
また改めてh'は

$\displaystyle h'(x_i-x_o',y_i-y_o')= A' \iint_{\infty} P(\lambda z_2 \xi', \lambda z_2 \eta') \exp{ \left[-i 2\pi \left[ (x_i-x_o') \xi'+ (y_i-y_o') \eta' \right] \right] } d\xi' d\eta' \tag{2-11} \label{2-11}$

です.

以上の内容をまとめると, 式(2-9)のように回折限界の光学系がつくる理想像は幾何光学に理想な像とインパルス応答の畳み込みとなる, そしてそのインパルス応答は式(2-6)のように回折の効果は瞳関数のフラウンホーファー回折像となる.

ここまではコヒーレント結像の場合でした. ではインコヒーレント結像ではどうなるでしょうか. インコヒーレント結像では物体各点からの振幅はインコヒーレントゆえに像面で干渉しません. その代わり, 各点からの強度が像面でそのまま足し算として最終的な強度を作ります.
つまり非常にざっくりと言えば式(2-9)の振幅の関係がそのまま強度の関係となります. 式(2-6)のインパルス応答自体は振幅に対する作用なので, 強度への作用はは式(2-6)を二乗だとすると,

$\displaystyle I_2(x_i,y_i)= | A'|^2 \iint_{\infty} |h'(x_i-x_o',y_i-y_o')|^2 I_{2g}(x_o',y_o' dx_o' dy_o' = |h'(x_o',y_o') |^2 \ast | I_{2g}(x_o',y_o')|^2 \tag{2-12} \label{2-12}$

となります.

光学伝達関数OTFとMTF

改めて前節の内容をまとめると, コヒーレント結像の場合は式(2-9)のように振幅の入出力が畳み込み積分で表される, つまり線形時不変システムだということ. インコヒーレントの場合は式(2-11)のように強度となります.
それぞれの式をフーリエ変換すると

$\displaystyle G_{U2}(f_X,f_Y)=H(f_X,f_Y)G_{Ug}(f_X,f_Y) \tag{3-1} \label{3-1}$
$\displaystyle G_{I2}(f_X,f_Y)=H'(f_X,f_Y)G_{Ig}(f_X,f_Y) \tag{3-2} \label{3-2}$

となります.

ここでHとH'はそれぞれの伝達関数となり, 光学ではそのまま光学伝達関数(OTF, Optical transfer function)と呼びます. コヒーレント結像の方のHを特に振幅伝達関数(ATF)と呼ぶこともあります.
インコヒーレント結像の場合さらにOTFの絶対値を取ったものを変調伝達関数(MTF, Modulation transfer function)と呼びます.
ここでさらにインコヒーレント結像の伝達関数は正規化することが普通なので, 式(3-2)のH'は

$\displaystyle H'(f_X,f_Y)=\frac{\mathcal{F} [ |h|^2 ] }{\iint_{\infty} |h'(x_i,y_I)|^2 dx_i dy_i} \tag{3-3} \label{3-3}$

となります. OTFを絶対値と位相部分で分解すると, 以下の式(3-4)が得られます. PTFが位相部分となりあまり光学で議論されることはないですが位相伝達関数と呼びます(Phase Transfer Function).

$\displaystyle H'(f_X,f_Y)=OTF(f_X,f_Y)=MTF(f_X,f_Y) \exp{ \left[ i PTF(f_X,f_Y) \right] } \tag{3-4} \label{3-4}$

数値計算例

せっかくなので数値計算例を示したいと思います. いったん波面収差はないとして(W=0), 焦点距離100mm, 瞳径5mm, F値20の軸上のPSFとMTFを計算するPythonコードが以下です.

import numpy as np
from matplotlib import pyplot as plt

def circ(t, width):
    #矩形開口を作成.
    return np.where(np.abs(t) <= width / 2, 1.0, 0.0)

M=1024
L=1e-3
dx=L/M

x=np.arange(-L/2,L/2-dx,dx)
y=x

wl=0.55e-6
k=2*np.pi/wl
Dxp=5e-3
wxp=Dxp/2
zxp=100e-3
lz=wl*zxp

fx=np.arange(-1/(2*dx),1/(2*dx),1/L)
Fx,Fy=np.meshgrid(fx,fx)

H=circ(np.sqrt(Fx*Fx+Fy*Fy),2*wxp/lz)

h2=np.abs(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(H))))**2

OTF=np.fft.fft2(np.fft.fftshift(h2))
MTF=np.abs(OTF)
MTF=np.fft.ifftshift(MTF)

PSF_cross_section= h2[int(M/2), :]
MTF_cross_section= MTF[int(M/2), :]

mask = fx >= 0
fx_filtered = fx[mask]
MTF_Normalized_temp = MTF_cross_section[mask]
MTF_Normalized = MTF_Normalized_temp/MTF_Normalized_temp[0]

# 2D PSF plot
plt.imshow(h2,extent=[x.min(), x.max(), y.min(), y.max()])
plt.colorbar(label="Value")
plt.title("2D PSF")
plt.show()

# PSF(cross section) plot
plt.scatter(x, PSF_cross_section)
plt.title("PSF")
plt.xlabel("cyc/m")
plt.ylabel("Intensity")
plt.show()

# MTF plot
plt.scatter(fx_filtered, MTF_Normalized)
plt.title("MTF")
plt.xlabel("cyc/m")
plt.ylabel("MTF")
plt.show()

まずPSFは以下の結果. 無収差なのでこのスケールだとグラフにする意味もあまりないですが, 瞳径を変えれば結果も変わりますので試してみください.

以下はMTFです.

今回はここまで. 次回は収差があるときを詳しく扱う予定です.

*1:例えば絞りが特殊な形状しているとか中間像があるときなど.

*2:この時プラスマイナスが違うことに注意.

*3:式(1-4)の内,瞳関数と波面収差の項だけを一般化瞳関数と呼びます.

*4:逆にいえばコヒーレント結像では要注意. ここらへんはGoodmanの本に説明があります.

光学設計とその周辺、そしてたまに全く関係ないやつ

学んだことを書き留めていきたいと思いますが、ありふれたことを書いても人類の進歩に貢献しないので、専門的な事柄をメインにしたいと思います。なお私の専門とは光学設計とか画像処理とかです。

Matlabに光学設計用のアドオンが登場 (Optical System Design and Analysis)

光学素子製造における補正(ターゲット)寸法の例

レンズ, フィルター厚み

曲率

熱膨張を考慮した金型設計

波動光学的にPSF, MTFをシミュレーションしてみる2(2/4) -光学系のモデル化とPSF, MTF-

光学系のモデル化

射出瞳モデルを使った光の伝搬とレンズのフーリエ変換作用

点像分布関数 PSF(Point Spread Function)と周波数応答関数

光学伝達関数OTFとMTF

数値計算例