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

以下の左図はベクトル場(y,x)のプロットです. 右図はこの図を45°回転させた様子です(表示範囲の違いのせいで見づらいですが).

この記事ではこんな感じでベクトル場を回転させて表示させる方法を紹介します. 記事にするほどのものではないでしょ, と思うかもしれませんが意外と考える時間が必要だったので, 同じような境遇にいる方のためにこのページにメモします.

なぜてこずったか

もちろん単に画像を回転させる方法は広く知られてます. 要は

$\begin{pmatrix} x'\\ y'\\ 1 \end{pmatrix} = \begin{pmatrix} \cos{\phi}&-\sin{\phi}&0\\ \sin{\phi}&\cos{\phi}&0\\ 0&0&1 \end{pmatrix} \begin{pmatrix} x\\ y\\ 1 \end{pmatrix} \tag{1} \label{1}$

の式で元画像(x,y)から $\phi$ だけ(x',y')に回転させることができます.

でも今やりたいのはベクトル場の画像を回転させる方法で本質的には同じですが少し工夫が必要だということに気づきました. しかもベクトルの回転とググっても大体出てくるのはベクトル場のrotの演算, ベクトル解析や電磁気学ででてきたあれだったりして, ただ無駄にネットの世界を漂流していました. もっとも結局最終的にたどり着いた資料はこの電磁気学の教科書だったわけですが....

本題

やることはべクトル画像を $\phi$ だけ回転させる, または元の座標xyから $\phi$ 回転した別の座標x'y'で表示させることです.

最初の図を例として, つまり以下のベクトルで表されるベクトル画像を考えていきます.

$\vec{L} = [ y, x ] \tag{2} \label{2}$

まず式(2)を以下のように単位基底ベクトルを使うとわかりやすくなります.

$\vec{L} = y\vec{x} + x\vec{y} \tag{3} \label{3}$

ここで $\vec{x}=[1,0]$ , $\vec{y}=[0,1]$ です.

さてこれを $\phi$ 回転させたい場合基本的には先の式(1)の変換式を使います. ただしここで注意が必要なのはスカラー部分だけでなく単位ベクトルの部分もこの式を使ってxyからx'y'に変換する必要があります.

つまりスカラー部分の変換
$x=\cos{\phi}x'-\sin{\phi}y', y=\sin{\phi}x' + \cos{\phi}y'$
に加えてベクトル部分の変換
$\vec{x} =\cos{\phi}\vec{x'} -\sin{\phi}\vec{y'},\vec{y}=\sin{\phi}\vec{x'} + \cos{\phi}\vec{y'}$
も行います.

これを式(3)に代入し, 計算を進めます.

$\vec{L} = (-\sin{\phi}x'+\cos{\phi}y')(\cos{\phi}\vec{x'} +\sin{\phi}\vec{y'}) +(\cos{\phi}x'+\sin{\phi}y')(-\sin{\phi}\vec{x'} +\cos{\phi}\vec{y'})$
$= (-\sin{2\phi}x'+\cos{2\phi}y')\vec{x'} +(\cos{2\phi}x'+\sin{2\phi}y')\vec{y'} \tag{4} \label{4}$

計算自体はこれで終わりです.

例としてこの式に $\phi=45$ °とし改めてxyで表すと

$\vec{L} = -x\vec{x} +y\vec{y} \tag{5} \label{5}$

となります. これを図にすると実際最初の右図のようになります.

内容は以上.

ベクトル画像の表示方法

これまで図で示してきたベクトル画像を作るPythonコードを以下示します(式(5)の場合). matplotlibのquiverという関数がキーです.

import numpy as np
import matplotlib.pyplot as plt

datanum=9
x = np.array([i-(datanum-1)/2 for i in range(datanum)])
y = np.array([i-(datanum-1)/2 for i in range(datanum)])

X, Y = np.meshgrid(x, y)

fig, ax = plt.subplots(figsize = (12, 7))

#最初の2つの引数が座標, 次の2つがベクトル.
ax.quiver(X, Y, -X, Y, linewidths=0.5, headlength=3, headaxislength=3)
 
ax.axis([-6, 6, -6, 6])

ax.set_xticks(x)
ax.set_yticks(y)

ax.grid(which = "major", axis = "x", color = "black", alpha = 1, linestyle = "-", linewidth = 1)
ax.grid(which = "major", axis = "y", color = "black", alpha = 1, linestyle = "-", linewidth = 1)

plt.show()

マシンビジョンのようにカメラを使って距離的な測量をするのはイメージングの重要なアプリケーションの1つです.
精度よく行うためにはカメラ自体のイメージング系がどういった幾何的なパラメーターを持っているかを調べる必要がありますが, その目的のカメラキャリブレーションという確立した手法が画像処理の世界にあります. ここで使われているカメラレンズや収差の考え方が光学設計の考えとほほぼ同等ですが, 微妙な考え方の違いなどがあったりとおもしろいので, 特にディストーションについていろいろと考えてみました.

キャリブレーションの種類

最初にちょこっとだけ説明するとイメージングのキャリブレーションというと,
幾何的な補正
光学的な補正
の2つが正確にはあります. このうち光学的キャリブレーションは像面照度の均一性などの話で純粋な光学寄りの内容となっており, いわゆるカメラキャリブレーションの範囲ではないと思いますので割愛.

カメラモデルの概要

そのディストーションのモデルとして以下の式がよく用いられています(参考文献1). 式(1-1),式(1-2)が半径方向, 後の2つディストーションとが円周方向のディストーションと呼ばれます.

$\displaystyle x_r'=x(1+k_1 r^2+k_2 r^4+ k_3 r^6) \rightarrow \Delta x_r'= x_r' - x=x(k_1 r^2+k_2 r^4+ k_3 r^6) \tag{1-1} \label{1-1}$
$\displaystyle y_r'=y(1+k_1 r^2+k_2 r^4+ k_3 r^6) \rightarrow \Delta y_r'= y_r' - y=y(k_1 r^2+k_2 r^4+ k_3 r^6) \tag{1-2} \label{1-2}$

$\displaystyle x_t'=x+[2p_1 xy+p_2(r^2+2x^2)] \rightarrow \Delta x_t'= x_t' - x= 2p_1 xy+p_2(r^2+2x^2) \tag{1-3} \label{1-3}$
$\displaystyle y_t'=y+[p_1(r^2+2y^2)+2p_2 xy] \rightarrow \Delta y_t'= y_t' - y=p_1(r^2+2y^2)+2p_2 xy \tag{1-4} \label{1-4}$

半径方向のディストーションはもっとシンプルに

$\displaystyle r_r'=r(1+k_1 r^2+k_2 r^4+ k_3 r^6) \rightarrow \Delta r_r'= r_r' - r=k_1 r^3+k_2 r^5+ k_3 r^7 \tag{1-5} \label{1-5}$

とも書けます.

またカメラキャリブレーションではセンサ側の座標(X,Y)と物体側の座標(x,y)は以下のカメラ内部パラメーター行列を介して変換します. これは縦横方向の焦点距離( $f_x, f_y$ )の違い(アナモルフィック)と画像シフト( $c_x, c_y$ )の補正の効果がありますが, 軸対称な光学系でも実際の製造誤差によってはこれらのパラメータが誤差モデルとして捉えます.

$\begin{pmatrix} x\\ y\\ 1 \end{pmatrix} = \begin{pmatrix} f_x&0&c_x\\ 0&f_y&c_y\\ 0&0&1 \end{pmatrix} \begin{pmatrix} X\\ Y\\ 1 \end{pmatrix} \tag{1-6} \label{1-6}$

カメラキャリブレーションという作業は上記のディストーション係数とカメラ内部パラメーターの値を実際の測定から得る手法です. 以下これらモデルと収差理論の関係を考え, このカメラモデルの妥当性を理論的に得ようと思います.

半径方向のディストーション

半径方向の歪曲は比較的理解しやすいでしょう. 式(1-5)の通り, 画像中心からの距離の3乗, 5乗, 7乗に比例してズレが大きくなる歪み方です. 多くの文献で取り上げられていますのでこの記事ではこれ以上は扱いません.

円周方向のディストーション

次の円周方向の歪曲はいろいろとやっかいです. 以下の図が式(1-3), (1-4)で左から $p_1$ の値があるとき( $p_2=0$ ), $p_2$ の値があるとき( $p_1=0$ ), $p_1$ と $p_2$ もどちらも値があるときの様子を示しています.

いろんな文献を調べましたがこの円周方向のディストーションについて詳しく説明している文献はあまり見つからず, 多くの画像処理の文献ではそもそも「実際の光学系では円周方向の歪曲が強く発生することはないから無視することが多い」と扱われる現状. そこで原典とされている参考文献2を見ると, 実際に偏心した状態に対して光線解析をすることで理論的に式を導いていますが式(1-3), (1-4)とは少し表し方が違っています. またこの文献では円周方向という呼び方はされていません.

なぜこれを円周方向と呼ぶかについては, 少しネットの世界をさまようと以下のQ&Aサイトがヒット.
optics - What is the "tangential" distortion of OpenCV actually tangential to? - Physics Stack Exchange

このページによると参考文献1の内容はその後いろんな方に修正を加えられることで, OpenCVの式になったらしく, その中でこの円周方向(接線方向)という呼び方になったようです. それでもこの質問者の方も言ってますが, どうもこの円周方向という言い方は誤解を生みそうな...

呼び方にこだわっても得るものはなさそうなので, 背景はわかったところでこのディストーションを頑張ってモダンな光学設計の観点から解釈しようと思います. 先述したとおり光学系が偏心した際に発生するディストーションですので, いわゆる偏心収差論を使ってこのディストーションの理論を再構築してみようと思います*1.

参考文献3,4によると, 光学系が偏心した際(単一の要素でも, ブロックとなっている要素でも)発生する収差は以下のリストにまとめられます. $\rho$ は動径成分, $H$ は画角, $E$ は偏心量(チルト, シフトまとめて)です. $E$ に比例する項に対して $E^2$ に比例する項は高次の偏心収差になるため, 通常はEの項のほうが影響が大きいです.

全部考えるのは大変なのでピンホールカメラの考えに則り $\rho$ に依存しない項だけ考えたいと思いますが,そうするとプリズム作用と2つの偏心歪曲, プリズム分光が対象になります. といってもプリズム作用は像位置のシフト, プリズム分光はRGB毎の像位置のシフトとしてカメラ内部パラメーターのほうで補正できそうなため, 偏心歪曲のみ考えてみようと思います. いや, どう考えても名前的にこの収差が一番関係していそうでしょという感じですけど.

名称	内容	依存性
プリズム作用	一律の像ずれ	$E$
軸上コマ	一律のコマ収差の発生	$E \rho ^2$
偏心非点収差	非対称な非点収差と像面の傾き(いわゆる方ぼけ)	$E \rho H$
偏心歪曲1	偏心時の非対称な歪曲	$E H ^2$
軸上非点収差	一律の非点収差の発生	$E^2 \rho$
偏心歪曲2	偏心時の非対称な歪曲(アフィン変換型偏心歪曲)	$E^2 H$
プリズム分光(色収差)	波長に依存した一律の像ずれ	$E$

この偏心歪曲は以下の式で表されます.

$\displaystyle \Delta Y=-\frac{E}{2} Y'^2 [(2+\cos{2\theta})V_{E(1,1)}-V_{E(2,1)}] - \frac{E^2}{2} Y'\cos{\theta}(3V_{E(1,2)}-2V_{E(2,2)}) \tag{2-1} \label{2-1}$
$\displaystyle \Delta X=-\frac{E}{2} Y'^2\sin{2\theta}V_{E(1,1)} - \frac{E^2}{2} Y'\sin{\theta}V_{E(2,2)} \tag{2-2} \label{2-2}$

ここで座標XYZは以下の図を参照. Y'が理想像高さ, Vとあるのが各収差係数です. $\theta$ が物体面(像面)の偏心の方向Yと物点のなす角度です(下左図). 以下の右図は $\theta=0$ となっているような状況です. またY'は画角角度 $\omega$ に比例します.

まずEに比例する項を先に考えていきたいのですが, 式(2-1)で $V_{E(2,1)}$ の項は一方向の特徴のない項のため無視するとして, 極座標 $Y=r\cos{\theta}, X=r\sin{\theta}$ となりますので*2, これを式(2-1)に代入すると, 係数はまとめて,
$\displaystyle \Delta Y = A Y'^2 (2+\cos{2\theta}) =A Y'^2 (2+\frac{Y^2}{r^2}-\frac{X^2}{r^2}) =A Y'^2 \frac{3X^2+Y^2}{r^2} \tag{2-3} \label{2-3}$
$\displaystyle \Delta X=A Y'^2\sin{2\theta} =A Y'^2 2\frac{XY}{r^2} \tag{2-4} \label{2-4}$
となります.
理想像高Y'と物体側rの値は比例するはずなので, 係数も改めてまとめると最終的に以下の表式になります.
$\displaystyle \Delta Y = A' (3X^2+Y^2) \tag{2-5} \label{2-5}$
$\displaystyle \Delta X=A' 2XY \tag{2-6} \label{2-6}$
これは見事に最初の式(1-3), (1-4)で $p_1=0$ としたときに相当します. 実際の式(1-3), (1-4)はこれを偏心の方向を任意の方向にしたものに相当します.

次に $E^2$ に比例する項を考えますが,
$\displaystyle \Delta Y=C Y'\cos{\theta} =CY' \frac{Y}{r}=C' Y \tag{2-7} \label{2-7}$
$\displaystyle \Delta X=D Y'\sin{\theta} =DY' \frac{X}{r}=D' X \tag{2-8} \label{2-8}$
となり, 要はY方向のずれはそのYの値に比例するということで, これ結局どういうことかというと, 焦点距離がYとX方向別に変化する, つまりアナモルフィック系のようなディストーションの出方になります. よってこの項はディストーションモデルというよりは内部パラメーター行列で補正できるような現象となります. ただカメラキャリブレーションのほうは座標系の水平方向, 垂直方向の固定された方向のパラメーターになっている一方, 偏心収差のほうは偏心の発生する360全方向に向きがあるため, 標準的なopencvの機能では補正することは恐らくできず.

こんな感じで, 円周方向ディストーションを偏心収差論を使うことである程度説明することができました. さて光学設計の観点からカメラモデルを考察すると, 例えば式(2-1)の $V_{E(1,1)}$ は偏心が無い場合の軸対称光学系の歪曲収差係数と非点収差係数に依存します. つまり, カメラモデルでは半径方向のディストーション係数kと円周方向の係数pは独立したパラメーターですが, 光学設計の理論から厳密には独立したパラメーターとはならず, 極端なケースでは数値の安定性などの影響があるかもしれません. こんな感じで収差論に基づいたカメラキャリブレーションの最適化アルゴリズムの設計を考えるような研究があると面白いと思うのですが, 実際あったりするんでしょうか?

実際デジカメを使うような場合は円周方向のディストーションは無視するのは妥当だと思いますが, 例えばコンバーターレンズなどのアタッチメントレンズを付けるときや設計的に偏心をしているようなイメージングを使う場合は考慮する必要があるはずです.

また一部文献には円周ディストーションはセンサー平面とレンズが平行ではないとき, つまりチルトしているとき発生すると説明されてますが, レンズ系の一部が平行偏心しているときも発生します. 要はエレメントの1つでも光軸に対称ではないような状況で発生するよ, ということでございまする.

参考文献
1)OpenCV document, https://docs.opencv.org/4.x/dc/dbb/tutorial_py_calibration.html
2)Conrady, A. E. "Lens-systems, decentered." Monthly Notices of the Royal Astronomical Society, Vol. 79, p. 384-390 79 (1919): 384-390.
3)松居吉哉, 偏心の存在する光学系の3次の収差論, 日本オプトメカトロニクス協会
4)木村研一,秋山健志, 浜野博之:“防振光学系への収差論の応用", 第 19 回光学シンポジウム予稿集 (1994) pp.47-50. 10

*1:偏心収差論もだいぶクラシックですけど

*2:参考文献の表記と合わせるためにcosとsinがよくある定義とは逆になっています

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

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

ベクトル画像の回転方法