写真レンズのレイトレーシング

夏休み中に写真レンズをレイトレーシングすることを行っていたので、その方法論と知見についてまとめてみました。この記事では写真レンズに対するレイトレーシングについて説明した後、Pythonによる実装を行い、球面収差や光線の経路を表示することを行っています。写真レンズの性能評価の方法についてはまた別記事で紹介したいです。

写真レンズに対するレイトレーシング

元々写真レンズの設計では 光線追跡法(レイトレーシング) がレンズの性能評価のために古くから使われていました。

レンズ設計者はレイトレーシングによって物体から出た光線(レイ)の像面への到達位置を計算し、収差(Aberration) や PSF(点像分布関数), MTF(変調伝達関数) といった量を計算することで性能評価を行っています。

以下はPSFの例です。PSFは点光源から出た光が像面でどのように結像するのかを表しています。完全な点になればボケのない像が得られるわけですが、実際には収差の影響で拡がりのある形になります。レンズ設計者の目標はこの拡がりを出来るだけ抑えることになります。

3DCGの文脈ではレイトレーシングというと写実的なCGを計算するための方法論を意味しますが、レンズ設計の文脈では光線の屈折・反射経路の計算方法を差しているので、意味合いが異なる点に注意します。この記事で紹介するレイトレーシングは後者に該当します。

3DCGへの応用

ここで紹介するレイトレーシングは3DCGにも応用することができます。レイトレでは最初にレイを飛ばすところ(DirectX RayTracingでいうRay Generation)でカメラモデルが必要になりますが、ここに写真レンズのレンズトレーシングを組み込むことで、まるでその写真レンズで撮影したかのような画像を出力することができます。

以下はその例です。この例ではSpectral Renderingを行っているため波長によって屈折率が変化し、色収差も表現できています。また、レンズ面での反射も考慮することでレンズフレアも再現できています。

これを行うコードはGithubで公開しているので気になる方は試してみて下さい。

yumcyaWiz/LensSim

これは余談ですが、光学設計の文脈で言う照明解析はモンテカルロ法を使っているそうなので、実質パストレなんじゃないかって思ってます。

レンズデータの表現

前置きが長くなりましたが実装の説明に入っていきます。まず写真レンズをプログラム中でどのように表現するかについて説明します。

座標系

レンズを配置していく基準となる座標系について定義しておきます。

Z軸をレンズ中心を貫く 光軸(Optical Axis) に置きます。+Z方向を像面側、-Z方向を物体側とします。

Y方向を上下、X方向をYZ軸に直交する方向(右手系)として定義します。

座標原点は像面の位置に取ることにします。

レンズ面の表現

レンズ面は以下の2種類で表現されることが多いです。

球面レンズ(Spherical Lens)
非球面レンズ(Aspheric Lens)

球面レンズとは文字通りレンズ面が球面になっているレンズです。写真レンズの大半は球面レンズで構成されています。

非球面レンズはレンズ面が球面ではないレンズです。以前は加工の難易度が高く製造が難しかったそうですが、近年は写真レンズ、望遠鏡、スマホなど至るところで使われています。

今回は球面レンズに絞って説明を行いたいと思います。非球面レンズの取り扱いについてはまた別記事で紹介したいです。

球面レンズの表現

球面レンズは球の半径 $r$ と、開口半径(Aperture Radius) $h$ で表現されます。

$r$ が正のときは図のように物体側に凸であるとします。 $r$ が負のときは物体側に凹であるとします。

レンズ設計の文脈では $r$ を 曲率半径(Curvature Radius) と呼ぶことが多いです。この記事でも以後そのように呼びます。

レンズ系の表現

次に複数の球面レンズ面で構成されるレンズ系について見ていきます。

単レンズの場合

まずは単レンズの場合について考えます。

$d$ はレンズの 厚み(Thickness) を表します。図のようにして1つのレンズの形状は第1面の形状を表す $r_1$ , $h_1$ と、第1面と第2面の間隔 $d$ 、第2面の形状を表す $r_2$ , $h_2$ の5つのパラメーターで表現することができます。

さらにレンズはガラスの種類をパラメーターとして持ちます。レイトレーシングの計算上は 屈折率(IOR) が求まっていれば良いので、ここでは簡単に特定の波長(例えばd線の587.56nm)に対する屈折率 $n$ を与えることにします。

以上をテーブルにまとめると次のように表現できます。

Curvature Radius	Aperture Radius	Thickness	IOR
$r_1$	$h_1$	$d$	$n$
$r_2$	$h_2$	0	1

第2面のThicknessは次の面が存在しないので0になっています。また屈折率は次の面が空気なので1になります。

レンズ系

複数の単レンズが集まったレンズ系の表現について考えます。

これは先ほどの単レンズの表現を単に並べていくだけです。例えば以下のような ダブルガウス型レンズ(Double Gauss Lens) は次のようなテーブルで表現することができます。

Curvature Radius	Aperture Radius	Thickness	IOR
29.475	12.6	3.76	1.67
84.83	12.6	0.12	1
19.275	11.5	4.025	1.67
40.77	11.5	3.275	1.699
12.75	5.705	18	1
...	...	...	...

ここで図中にある青線の部分は 絞り(Stop) に対応します。絞りは一般的に平面ですが、 $r$ が無限大の球面と考えればレンズ系全体を球面だけで考えることができます。ただし、数値計算上は $r$ が非常に大きいと計算誤差が大きくなることがあるため、平面フラグをつけて球面と平面を区別する方が良いです。絞りだけではなく、平面となっているレンズ面も同様に扱われます。

写真レンズの設計データは一般的にこのようなテーブル形式で表現されることが多いです。

ここで紹介したダブルガウス型レンズの設計データはScenes for pbrt-v3から入手することができます。

レイトレーシング

まずはアルゴリズムの全体的な流れについて見ます。

写真レンズに対するレイトレーシングは次のような流れで行われます。

球面(or 平面)とレイの衝突位置を計算。同時に法線も算出
スネルの法則に従い屈折方向を計算
レイを屈折方向に飛ばし、1に戻る

これを繰り返せば最終的にレイは最終面から出発し、像面との交差位置を計算することができます。

なお屈折時にはフレネル反射や全反射が起こる可能性があります。これを考慮すればレンズフレアを再現することができますが、実装が少し複雑になるので、ここでは考慮しないことにします。

各手順について説明していきます。

レンズ系の位置

まずレンズ系の位置を適当な位置に配置し、それぞれのレンズの中心位置を決定することを考えます。ここで言うレンズの中心位置とはレンズ面が光軸(z軸)と交わる位置を言います。

いま最終面から像面までの距離を $s'$ とします。この値を最終面のThicknessとすることにします。

座標系の取り方から、像面は原点に置かれているのでした。光軸は $z$ 軸に一致しているので、最終面の中心位置は $(0, 0, -s')$ となります。

他の面についてはThicknessを最終面から最前面に向かって足し合わせていくことで $(0, 0, -z)$ のように与えられます。

球面とレイの衝突位置・法線の計算

レイトレでおなじみの計算です。ベクトル方程式を作って2次方程式に帰着して解析解を求め、それを数値計算します。

考えている球面の曲率半径を $r$ , 開口半径を $h$ とします。いまレンズ系を適当な場所に配置し、各レンズ面の中心位置 $\vec{l}$ が求まっているとします。

このとき球面の中心位置 $\vec{c}$ は

\vec{c} = \vec{l} + (0, 0, r)

となります。

球面上の点を $\vec{p}$ とすると、球面を表す方程式は以下で与えられます。

\|\vec{p} - \vec{c}\| = r

レイの始点を $\vec{o}$ 、レイの方向を $\vec{d}$ とします。レイが球面と交差すると仮定すると $\vec{p} = \vec{o} + t\vec{d}$ となります( $t$ は実数)。これを上の式に代入すると

\|\vec{o} + t\vec{d} - \vec{c}\| = r

両辺を二乗すると次のような $t$ に関する二次方程式が得られます。

\|\vec{d}\|^2 t^2 + 2\vec{d}\cdot(\vec{o} - \vec{c})t + \|\vec{o} - \vec{c}\|^2 - r^2 = 0

解の公式より $t$ を次のように計算できます。( $\vec{d}$ は正規化されていると仮定します。つまり $\|\vec{d}\| = 1$ )

\begin{aligned} b &= \vec{d}\cdot(\vec{o} - \vec{c}) \\ c &= \|\vec{o} - \vec{c}\|^2 - r^2 \\ D &= b^2 - c \\ t_1 &= -b - \sqrt{D} \\ t_2 &= -b + \sqrt{D} \end{aligned}

ここで $D$ は判別式です。 $D < 0$ は解なしで、レイと球面が交差しない状態に対応するのでレイトレーシングをここで終了します。

$D > 0$ の場合はレイと球面が2箇所で交差するので、 $t_1$ と $t_2$ の２つの解が存在します。どちらの解を選ぶかですが、これはレイの進行方向 $\vec{d}$ の $z$ 成分 $d_z$ と曲率半径 $r$ の符号によって次のように選びます。

$d_z$ の符号	$r$ の符号	解
$+$	$+$	$t_1$
$+$	$-$	$t_2$
$-$	$-$	$t_1$
$-$	$+$	$t_2$

以上より $t$ の値が求まるので、レイと球面の交差位置 $\vec{p}$ を

\vec{p} = \vec{o} + t\vec{d}

によって計算できます。

次に衝突位置が開口半径以下であるかをチェックします。以下の式が成立していればレイは球面と交差しています。

p_x^2 + p_y^2 \le h^2

この式が成り立っていない時はレイトレーシングを終了します。

最後に法線 $\vec{n}$ は球面の中心位置 $\vec{c}$ を使って次のように計算できます。

\vec{n} = \mathrm{Normalize}(\vec{p} - \vec{c})

$\mathrm{Normalize}$ は正規化を意味しています。ここで法線の向きに注意する必要があります。ここでは法線の向きをレンズ面から入射側への方向に合わせることにします。以下の図で言うと点線になっているものが正しい法線に対応し、実線は上記の計算式で計算した法線に対応します。

図から分かる通り、 $r$ の正負で法線の向きは逆になってしまいます。法線の向きが間違っていると屈折方向の計算がおかしくなってしまうので、必ず正しい向きを選ぶ必要があります。

解決策として、レイの方向ベクトルと法線の内積を取り、その符号の正負で法線の向きを逆にする方法があります。例えば上の図では計算した法線の向きが逆になっていますが、このとき内積の値は正になっています。そこで内積の値が正のときは法線を-1倍したものを新たに法線として返すようにします。

平面との交差位置・法線の計算

続いて平面とレイの交差位置の計算です。

平面上の点を $\vec{p}$ 、平面の中心を $\vec{l}$ 、法線を $\vec{n}$ とすると、平面は以下のようなベクトル方程式で表現できます。

(\vec{p} - \vec{l})\cdot\vec{n} = 0

始点 $\vec{o}$ , 方向 $\vec{d}$ を持つレイが平面と交差すると仮定し、 $\vec{p} = \vec{o} + t\vec{d}$ ( $t$ は実数)とすると

(\vec{o} + t\vec{d} - \vec{l})\cdot\vec{n} = 0

これを $t$ に関して解くと

t = -\frac{(\vec{o} - \vec{l})\cdot\vec{n}}{\vec{d}\cdot\vec{n}}

特に平面レンズの場合を考えて、 $\vec{l} = (0, 0, l_z)$ , $\vec{n} = (0, 0, 1)$ のとき

t = -\frac{o_z - l_z}{d_z}

となります。

あとは球面の時と同様に $t$ の値から交差位置 $\vec{p} = \vec{o} + t\vec{d}$ を計算し、その点が開口半径以下であるかをチェックすれば良いです。

法線は $(0, 0, 1)$ を用いれば良いです。ただし球面の時と同様にレイの方向ベクトルに応じて逆転させたりしなかったりする必要があります。

屈折方向の計算

スネルの法則(Snell's Law) に従って計算します。 $n_1$ を入射側媒質の屈折率, $n_2$ を屈折後の媒質の屈折率、 $\theta_1$ を入射角、 $\theta_2$ を屈折角とすると、スネルの法則は次のような式でした。

n_1\sin{\theta_1} = n_2\sin{\theta_2}

この式を使って屈折方向を計算することを考えます。入射方向を $\vec{v}$ 、媒質の界面の法線を $\vec{n}$ 、屈折方向を $\vec{t}$ とします。また、方向ベクトルの水平成分と垂直成分をそれぞれ下付きの $_h$ , $_p$ と表すことにします。方向ベクトルは正規化されていると仮定します。

ここで入射方向 $\vec{v}$ はレイの方向 $\vec{d}$ を-1倍したものになっていることに注意してください。こっちの方が導出が見やすいのでそうしました。

図から入射ベクトルの水平成分 $\vec{v_h}$ の大きさは $\sin{\theta_1}$ に、垂直成分 $\vec{v_p}$ の大きさは $\cos{\theta_1}$ に等しいことが分かります。 $\vec{t}$ についても同様です。この関係にスネルの法則を適用します。

まず内積から

\cos{\theta_1} = \vec{v}\cdot\vec{n}

と計算できます。これを用いて $\vec{v}_h$ は

\begin{aligned} \vec{v}_h &= \vec{v} - (\cos{\theta_1})\vec{n} \\ &= \vec{v} - (\vec{v}\cdot\vec{n})\vec{n} \end{aligned}

と計算できます。 $\|\vec{v_h}\| = \sin{\theta_1}$ です。

スネルの法則より

\|\vec{t_{h}}\| = \sin{\theta_2} = \frac{n_1}{n_2}\sin{\theta_1} = \frac{n_1}{n_2}\|\vec{v_h}\|

$\vec{t_h}$ の向きは $\vec{v_h}$ の逆なので

\begin{aligned} \vec{t_h} &= -\frac{n_1}{n_2}\vec{v_h} \\ &= -\frac{n_1}{n_2}(\vec{v} - (\vec{v}\cdot\vec{n})\vec{n}) \end{aligned}

と計算できます。

次に $\vec{t_p}$ を求めます。図より

\vec{t_p} = -(\cos{\theta_2})\vec{n}

$\cos{\theta_2}$ は三角関数の公式から $\sin{\theta_2}$ を用いて次のように計算できます。

\begin{aligned} \cos{\theta_2} &= \sqrt{1 - \sin^2{\theta_2}} \\ &= \sqrt{1 - \|t_h\|^2} \end{aligned}

ここで $\|\vec{t_h}\| > 1$ のときは全反射に対応し、その時はレイトレーシングを終了します。まとめると

\vec{t_p} = -(\sqrt{1 - \|t_h\|^2})\vec{n}

最後に $\vec{t}$ は次のように計算できます。

\begin{aligned} \vec{t} &= \vec{t_h} + \vec{t_p} \\ &= -\frac{n_1}{n_2}\vec{v} + (\frac{n_1}{n_2}(\vec{v}\cdot\vec{n}) - \sqrt{1 - \|t_h\|^2})\vec{n} \end{aligned}

屈折後のレイの生成

上記の手順より球面との衝突位置 $\vec{p}$ , 屈折後の方向 $\vec{t}$ が求まったので、次のレイは始点を $\vec{p}$ 、方向を $\vec{t}$ として飛ばすだけです。

ここで紹介したレイトレーシングの方法はレンズ設計の教科書には書かれていない方法で、あまり一般的ではないようです。レンズ設計の教科書ではレイを違う形式で表現し、成分ごとに計算式があって煩雑な印象があります。個人的にはここで紹介した方法の方がベクトルを使ってシンプルに書けて好きです。

Pythonによる実装

写真レンズのレイトレーシングをPythonを用いて実装していきます。ベクトルの表現にはnumpyを用いることにします。

ここで紹介するPython実装はyumcyaWiz/raytracing-for-lensで公開しています。実行例はJupyter Notebookとして見ることができます。気になる方は色々試してみて下さい。

レイの表現

レイを表現するRayクラスを用意します。

レイ上の始点から距離 $t$ にある点の位置を計算できると何かと便利なため、それ用のposition()も用意しておきます。

class Ray:
    # レイを表現する
    # origin: レイの始点
    # direction: レイの方向ベクトル
    def __init__(self, origin: np.ndarray, direction: np.ndarray):
        self.origin = origin
        self.direction = direction

    def __repr__(self):
        return "origin: {0}, direction: {1}".format(self.origin, self.direction)

    # レイ上の位置を計算する
    # t: 始点からの距離
    def position(self, t: float) -> np.ndarray:
        p = self.origin + t * self.direction
        return p

ユーティリティ

ベクトルを正規化するnormalize()と、レイの方向に応じて必要があれば法線をひっくり返してくれるflip_normal()を用意しておきます。

# ベクトルを正規化する
def normalize(v: np.ndarray) -> np.ndarray:
    return v / np.linalg.norm(v)

# レイの方向に応じて法線を逆転させる
# v: レイの方向ベクトル
# n: 法線
def flip_normal(v: np.ndarray, n: np.ndarray) -> np.ndarray:
    if np.dot(v, n) < 0:
        return n
    else:
        return -n

屈折方向の計算

受け取った入射方向と法線、屈折率から屈折方向を計算するrefractを用意しておきます。全反射が起こる時はNoneを返すようにします。先ほどの式をそのまま書いただけです。

# 屈折ベクトルを計算する
# 全反射の場合はNoneを返す
# v: 入射ベクトル
# n: 法線
# n1: 入射側媒質の屈折率
# n2: 出射側媒質の屈折率
def refract(v: np.ndarray, n: np.ndarray, n1: float, n2: float) -> Optional[np.ndarray]:
    # 屈折ベクトルの水平方向
    t_h = -n1 / n2 * (v - np.dot(v, n)*n)

    # 全反射
    if np.linalg.norm(t_h) > 1:
        return None

    # 屈折ベクトルの垂直方向
    t_p = -np.sqrt(1 - np.linalg.norm(t_h)**2) * n

    # 屈折ベクトル
    t = t_h + t_p

    return

レンズデータの表現

レンズデータはテーブル形式で表現できることを先ほど見ました。ここでは以下のようなcsvファイルでレンズを表現することにします。

r	h	d	ior
29.475	12.6	3.76	1.67
84.83	12.6	0.12	1
19.275	11.5	4.025	1.67
40.77	11.5	3.275	1.699
12.75	5.705	18	1
...	...	...	...

rは曲率半径, hは開口半径, dは厚み, iorは屈折率に対応します。データは物体側から像面側に向かって並べられています。

絞りのようにレンズ面が平面であるものは $r = 0$ で表現することにします。

csvファイルの読み込みにはpandasを利用します。

レンズ面の実装

レンズ面はLensSurfaceクラスで表現することにします。メンバ変数として曲率半径と開口半径、屈折率、それにレンズ面の中心位置のz座標を表すzを持たせておきます。zの値は他の面の情報もないと決定できないので、後ほどレンズ系を読み込む段階で初期化するようにします。

class LensSurface:
    # レンズ面を表現する
    # r=0で平面を表現する
    # r: 曲率半径
    # h: 開口半径
    # d: 次の面までの距離
    # ior: 屈折率
    def __init__(self, r: float, h: float, d: float, ior: float):
        self.z = 0
        self.r = r
        self.h = h
        self.d = d
        self.ior = ior

レイとレンズ面の交差位置・法線を計算するintersect()を実装します。 $r$ の値が0のときは平面との交差判定、それ以外のときは球との交差判定を行います。交差しないときはNoneを返すようにします。

    # レイとの交差位置, 法線を計算する
    # 交差しない場合はNoneを返す
    # ray: レイ
    def intersect(self, ray: Ray) -> Optional[Tuple[np.ndarray, np.ndarray]]:
        if self.r != 0:
            # 球面との交差

            # レンズの中心位置
            center = np.array([0, 0, self.z + self.r])

            # 判別式
            b = np.dot(ray.direction, ray.origin - center)
            c = np.linalg.norm(ray.origin - center)**2 - self.r**2
            D = b**2 - c

            # D < 0の場合は交差しない
            if D < 0:
                return None, None

            # tの候補
            t_1 = -b - np.sqrt(D)
            t_2 = -b + np.sqrt(D)

            # 適切なtを選択
            t = None
            if ray.direction[2] > 0 and self.r > 0:
                t = t_1
            elif ray.direction[2] > 0 and self.r < 0:
                t = t_2
            elif ray.direction[2] < 0 and self.r < 0:
                t = t_1
            else:
                t = t_2

            # 交差位置
            p = ray.position(t)

            # 交差位置が開口半径以上なら交差しない
            if p[0] ** 2 + p[1] ** 2 > self.h ** 2:
                return None, None

            # 法線
            n = flip_normal(ray.direction, normalize(p - center))

            return p, n
        else:
            # 平面との交差

            # 交差位置
            t = -(ray.origin[2] - self.z) / ray.direction[2]
            p = ray.position(t)

            # 交差位置が開口半径以上なら交差しない
            if p[0] ** 2 + p[1] ** 2 > self.h ** 2:
                return None, None

            # 法線
            n = flip_normal(ray.direction, np.array([0, 0, -1]))

            return p, n

レンズ系の実装

レンズ系はLensSystemクラスで表現することにします。

レンズ系はメンバ変数としてLensSurfaceのリストを持つようにします。コンストラクタではレンズデータのcsvへのファイルパスを引数として取り、LensSurfaceのリストを初期化した後、各レンズ面の中心位置のz座標を初期化していきます。

class LensSystem:
    # レンズ系を表現する
    # filepath: csvファイルのファイルパス
    def __init__(self, filepath: str):
        # レンズデータの読み込み
        self.df = pd.read_csv(filepath)

        # レンズ面の生成
        self.lenses = []
        for i in range(len(self.df)):
            self.lenses.append(LensSurface(
                self.df.iloc[i]["r"],
                self.df.iloc[i]["h"],
                self.df.iloc[i]["d"],
                self.df.iloc[i]["ior"]
            ))

        # 各レンズ面の位置を計算
        z = 0
        for i in reversed(range(len(self.df))):
            z -= self.lenses[i].d
            self.lenses[i].z = z

    def __repr__(self):
        return str(self.df)

次は物体側から像面側に向かってレイトレーシングを行うraytracing_from_object()を実装します。この関数は途中のレイを全て保存してリストとして返すようにしています。こうすることで光が通る経路の様子が分かります。

実装はレイトレーシングの部分で説明した通りの流れになっています。注意点としてはrefract()に-ray.directionを入れることや、n1やrayの更新を忘れずにということぐらいです。

    # 物体側から像側に向かってレイトレーシングを行い、光線経路を返す
    # ray_in: 入射レイ
    def raytrace_from_object(self, ray_in: Ray) -> List[Ray]:
        n1 = 1
        ray = ray_in
        rays = [ray]
        for i in range(len(self.lenses)):
            lens = self.lenses[i]

            # レンズとの交差位置, 法線を計算
            p, n = lens.intersect(ray)
            if p is None or n is None:
                return None

            # 屈折方向を計算
            n2 = lens.ior
            t = refract(-ray.direction, n, n1, n2)
            if t is None:
                return None

            # レイを更新
            ray = Ray(p, t)
            rays.append(ray)

            # 屈折率を更新
            n1 = n2

        return rays

LensSystemを使ってみる

以上でLensSystemが実装できたので、早速使ってみることにします。

レンズデータとして先ほどのダブルガウス型レンズを用いることにします。このレンズの焦点距離は50.3581mmで、最終面から焦点までの距離(BFL)は36.1059mmと与えられています。

最終面から像面までの距離はBFLである36.1059mmとしておきます。こうすることで光軸に並行に入射した光は収差がゼロと仮定すれば、全て像面中心に集まることになります。これを確認してみましょう。

球面収差の確認

まずはLensSystemを作成します。

lsys = LensSystem('data/dgauss50mm.csv')

続いてレイトレしてみます。入射レイは始点を $(0, 1, -1000)$ , 方向を $(0, 0, 1)$ としておきます。

rays = lsys.raytrace_from_object(Ray(np.array([0, 1, -1000]), np.array([0, 0, 1])))

raysの中身はこんな感じです。originは光線の始点, directionは光線の方向を表しています。

[origin: [    0     1 -1000], direction: [0 0 1],
 origin: [  0.           1.         -68.12893159], direction: [-0.         -0.01361615  0.9999073 ],
 origin: [  0.           0.94895732 -64.38059204], direction: [-0.         -0.01524392  0.9998838 ],
 origin: [  0.           0.94685399 -64.2426296 ], direction: [-0.         -0.02884308  0.99958395],
 origin: [  0.           0.83113927 -60.2324273 ], direction: [-0.         -0.02869874  0.99958811],
 origin: [  0.           0.73674378 -56.94459626], direction: [-0.         -0.00833916  0.99996523],
 origin: [  0.          0.6893449 -51.2609   ], direction: [-0.         -0.00833916  0.99996523],
 origin: [  0.           0.65193973 -46.77556853], direction: [0.         0.01173173 0.99993118],
 origin: [  0.           0.66602006 -45.57545957], direction: [0.         0.01080045 0.99994167],
 origin: [  0.           0.7313281  -39.52902271], direction: [ 0.         -0.00571145  0.99998369],
 origin: [  0.           0.73016447 -39.32529009], direction: [ 0.         -0.00402404  0.9999919 ],
 origin: [  0.           0.71723546 -36.11237456], direction: [ 0.         -0.01985852  0.9998028 ]]

raysの最後尾にあるレイがレンズ系の最後尾から出射したレイに対応します。これを使って像面との交差位置を計算してみます。

# 像面との交点
t = -rays[-1].origin[2] / rays[-1].direction[2]
p = rays[-1].position(t)
print(p)

出力

[ 0.0000000e+00 -4.4162616e-05  0.0000000e+00]

y成分に注目すると、わずかに原点からズレていることが分かります。これはレンズが 球面収差(Spherical Aberration) を持つためです。

レイの入射高が変わると球面収差の量は変化します。これをグラフにしてみましょう。

一般的に球面収差はx軸に像面との交差位置のy座標、y軸にレイの入射高をとって表示することが多いです。LensSystemにこれをプロットするplot_spherical_aberration()を実装します。

    # 球面収差をプロットする
    def plot_spherical_aberration(self):
        graph_x = []
        graph_y = []

        for i in range(50):
            # 入射高
            u = (i / 50)
            height = 0.5 * self.lenses[0].h * u
            graph_y.append(height)

            # レイトレ
            rays = self.raytrace_from_object(Ray(
                np.array([0, height, -1000]),
                np.array([0, 0, 1])
            ))

            # 像面との交点
            t = -rays[-1].origin[2] / rays[-1].direction[2]
            p = rays[-1].position(t)
            graph_x.append(p[1])

        ax = plt.plot(graph_x, graph_y)
        plt.grid()
        plt.title('Spherical Aberration')
        plt.xlabel('$y\mathrm{[mm]}$')
        plt.ylabel('Height$\mathrm{[mm]}$')
        return ax

プロットしてみると次のようになります。

入射高が0mm付近と、10mm付近で球面収差の量が小さくなっていて、10mm以降では急速に球面収差の量が大きくなっていることが分かります。これは写真レンズでよく見られる形で、収差の補正を行った結果こういう形になります。

収差には球面収差以外にもコマ収差, 非点収差, 像面湾曲, 歪曲収差といったものがあります。これらの収差による影響もレイトレーシングによってグラフ化することができます。最初にも述べましたが、こういったレンズの性能評価についてはまた別記事で紹介したいと思います。

光路図の表示

球面収差のグラフだけを見てもレイトレしている実感が湧かないと思うので、光路図を表示してみることにします。

まずはレンズ系を表示するところから実装していきます。

レンズ面のプロット

絞りに関してはただのlineplotを利用すればよいです。開口半径以下の場所を青線で囲って表示することにします。

            # レンズ面のプロット
            # 絞りの場合
            if lens.r == 0:
                ax.plot([lens.z, lens.z], [lens.h, 1.2*lens.h], c='blue')
                ax.plot([lens.z, lens.z], [-lens.h, -1.2*lens.h], c='blue'

球面はmatplotlib.patchesにあるArcを利用します。Arcは楕円の中心位置、長径、短径、回転角、制限角などを必要とします。中心位置は $(0, 0, z + r)$ とすればよく、長径・短径は $2|r|$ , 回転角は $r$ の正負に応じて0または180, 制限角は $\sin^{-1}\frac{h}{r}$ で計算すれば良いです。

            # 球面レンズの場合
            else:
                theta = abs(np.degrees(np.arcsin(lens.h / lens.r)))
                angle = 180 if lens.r > 0 else 0
                arc = patches.Arc((lens.z + lens.r, 0), 2 * abs(lens.r),
                                  2 * abs(lens.r), angle=angle, theta1=-theta, theta2=theta)
                ax.add_patch(arc)

レンズ系のプロット

上のコードで個別の面を表示できるので、これをレンズ面全体に渡ってfor文を回してレンズ系を表示させます。しかし、これだけだと以下のような違和感のあるプロットが得られます。

違和感を感じるのはレンズ枠が表示されていないためです。これを表示することを考えます。

レンズ枠を表示する必要があるのはレンズが空気面を挟まずにくっついている部分です。なのでこれを判定してプロットするかしないかを分けます。

            if i > 0:
                lens_prev = self.lenses[i - 1]

                # 前面 or 現在の面が絞りならスキップ
                if lens.r == 0 or lens_prev.r == 0:
                    continue

                # 前面が空気ならスキップ
                if lens_prev.ior == 1:
                    continue

レンズ枠は前のレンズ面と次のレンズ面の間を平行線でつなぐことで表示します。このときレンズ面の開口半径が高い方のy座標の値を用いてプロットすることにします。compute_lでx座標の値を計算しています。

                def compute_l(z, r, h):
                    return r - np.sqrt(r**2 - h**2) if r > 0 else -(np.abs(r) - np.sqrt(r**2 - h**2))

                zp = lens_prev.z
                hp = lens_prev.h
                lp = compute_l(lens_prev.z, lens_prev.r, lens_prev.h)
                z = lens.z
                h = lens.h
                l = compute_l(lens.z, lens.r, lens.h)

                if lens.h > lens_prev.h:
                    ax.plot([zp + lp, z + l], [h, h], c='black')
                    ax.plot([zp + lp, z + l], [-h, -h], c='black')
                    ax.plot([zp + lp, zp + lp], [hp, h], c='black')
                    ax.plot([zp + lp, zp + lp], [-hp, -h], c='black')
                else:
                    ax.plot([zp + lp, z + l], [hp, hp], c='black')
                    ax.plot([zp + lp, z + l], [-hp, -hp], c='black')
                    ax.plot([z + l, z + l], [hp, h], c='black')
                    ax.plot([z + l, z + l], [-h, -hp], c='black'

レンズ枠を加えてプロットし直すと次のような図が得られます。いい感じです。

光路の表示

光路は何本かのレイをレイトレーシングして得られたraysの始点を線分で結んで表示すればよいです。このとき最終面から像面までの光路も表示すると見栄えが良いので、像面との交差位置を計算し、そこまで線分を延長します。

optical_path_diagram()に実装していきます。最初にレンズ系をプロットして、その後に光路をlineplotでプロットしています。

    # 光路図をプロットする
    # n_rays: レイの本数
    def optical_path_diagram(self, n_rays=10):
        # レンズ系の表示
        ax = self.plot()

        # 光路図のプロット
        for i in range(n_rays):
            u = 2 * (i + 0.5) / n_rays - 1
            h = self.lenses[0].h

            # 入射レイ
            ray_direction = np.array([0, 0, 1])
            ray_origin = np.array([0, u*h, self.lenses[0].z]) - ray_direction
            ray_in = Ray(ray_origin, ray_direction)

            # レイトレ
            rays = self.raytrace_from_object(ray_in)

            # zとyの抽出
            line_x = list(map(lambda x: x.origin[2], rays))
            line_y = list(map(lambda x: x.origin[1], rays))

            # 像面までの光路を追加
            if len(line_x) == len(self.lenses) + 1:
                if rays[-1].direction[2] != 0:
                    # 像面との交差位置を計算
                    t = -rays[-1].origin[2] / rays[-1].direction[2]
                    p = rays[-1].origin + t*rays[-1].direction

                    # zとyを追加
                    line_x.append(p[2])
                    line_y.append(p[1])

            # プロット
            ax.plot(line_x, line_y, c='lime')

        return ax

実行すると次のような図が得られます。見栄えが良いし、レイトレしている実感が湧きますね。

レンズ設計データの入手

最後にレンズの設計データを入手する方法について紹介しておきたいと思います。

大抵の写真レンズの設計データは 特許(Patent) として公開されています。従って特許データベースのようなものを検索してPDFなりを見れば良いです。

しかし、特許PDFは設計データだけが書いてあって、具体的に〇〇という名前で市販していますとかは書いていないので、一致するものを探すのに結構苦労します。私はSigmaのレンズが好きなので、Sigma 30mm F1.4 DC HSMというレンズを再現したくて対応する特許を探していたのですが、レンズ構成図やレンズ枚数, 画角などの仕様を特許文書の情報と比較して一致するものを探していました。

こんなことをしたくないという人は、以下の海外サイトがおすすめです。有名なオールドレンズや最近の写真レンズに対応する特許番号が表でまとめられていて、かなりの数があります。

Lens patents

今回紹介したダブルガウス型レンズは以下のpbrtシーンから引っ張ってきました。これも元はどっかの特許データから引っ張ってきているようです。ダブルガウス型以外にも広角レンズや魚眼レンズのデータもあるので、是非試してみて下さい。

Scenes for pbrt-v3

参考資料

この分野で参考となる資料について説明します。

以下は物理ベースレンダリングで有名なpbrtです。3DCGの分野でカメラモデルを組み込みたい時に参考になります。ここで紹介したレイトレーシングの方法は概ねこれに則っています。

Physically Based Rendering: Camera Models

次にレンズ設計の分野で参考になる本です。光学の理論的なところや写真レンズの性能評価、設計方法について知ることができます。日本語で読めるものだと以下がおすすめです。

レンズ設計法

洋書だと以下の本がよくまとまっています。この分野では標準的な教科書なのではないでしょうか。特にModern Optical Engineeringは網羅的なのに読みやすくて良かったです。