【このテーマの記事は、UAV写真測量について、日々の文献調査や研究で得た、PhotoScanに限らない一般的な情報を掲載していきます。用語の説明は「PhotoScanを極める」に譲ります。】

＜プロローグ＞

現在の標準的なUAV写真測量では、SfMの過程で内部パラメータを推定します。この「セルフキャリブレーション」におけるfの推定誤差が、大きな鉛直誤差を生むことがあります。対地高度100 mの撮影でfの推定が1‰ずれると、地表の高さを大雑把には10 cmほど見誤ることになります。

このfについて説明するとき、私は

投影中心と撮像面の距離を画素単位で表したもの
奥行き方向のスケールを表すパラメータ

のような表現を使っています。表現1はfの別称である「焦点距離」「画面距離」とよく合います。表現2はSfMにおける役割を強調した表現です。

しかし、fは直接的にはカメラモデルにおいて定義される1パラメータであり、上記はそれを「意訳」というか、わかりやすく表現したものです。ピンホールカメラは別として、レンズ歪みのあるカメラを使ったSfMにおいて、これらの表現は適切なのでしょうか。

＜歪みのあるカメラで、表現1は適切か＞

投影中心（光学中心; optical center）を、画像を構成する全ての光線が交わる点と定義するならば、

それはそもそも、理想的なピンホールカメラか、学校で習った理想的なレンズ（薄いレンズ）1枚のカメラにおいてしか存在しません。

SfMで歪みのあるカメラを扱う場合においても投影中心の語は用いられますが、これは現実の投影を

理想的なピンホールカメラによる投影と、レンズ歪みの影響に分離して捉えたとき（モデル化したとき）に生じる概念です。言い換えますと、歪みを考慮した「カメラモデル」上の概念になります。

例えば多くのSfM-MVSソフトでは、Brown (1971)などに由来する次のようなカメラモデルを用います

（ただし使われ方は逆方向です。織田 (2016)の2.4を参照）。

ここで、X, Y, Zは3次元空間のある点のカメラ座標系における座標、K1, K2, K3, P1, P2は歪みパラメータ、u,vは実際に観測される画素座標（画像中心基準）です。上式は、Agisoft Metashapeのカメラモデルにおいて、簡単のためcx = cy = 0, P3 = P4 = 0, B1 = B2 = 0とした場合の式です。ソフトによって、考慮するパラメータやその表現に多少の差があります（中野, 2016）。

式(1)のx, yは、歪みがなくf=1とした場合の画素座標です。一方、式(3)(4)で求められるu, vは、歪みの影響を受けた現実の画素座標です。つまり上記のカメラモデルは、

Step1: 式(1)によって理想的なピンホールカメラで撮影した場合の画素座標を求める
Step2: 式(3)(4)によって歪みの効果を加味している

という2段構えになっています。このStep1において、想像上のピンホールカメラの投影中心が使われるわけです。カメラ座標系自体が原点として投影中心を必要とします。あるカメラモデルを仮定し、内部・外部パラメータを適切に設定したときに、画像がもっともよく説明されるようなカメラ座標系の原点が、投影中心であると言うことも出来るでしょう。実際SfMではそのように投影中心を求めているわけです。

以上のように、歪みを考慮したSfMにおける投影中心はカメラモデル上の抽象的な概念であり、しかも、歪みを考慮する「前」の投影中心であって、歪みを1次近似したときの投影中心ですらありません。

従って表現1は、歪みを度外視する状況ではよいけれども、正確な表現とは言えないでしょう。

＜歪みのあるカメラで、表現2は適切か＞

「奥行き方向のスケール」という表現自体の意味は、こちらの記事を参照ください。

結論から書きますと、内部パラメータfが「奥行き方向のスケール」を表すという表現2は、表現1と同じく、理想的なピンホールカメラでは正しいけれども、現実のカメラでは正確ではありません。

上記のカメラモデルを使って数値例を作りましょう。

いま、カメラ座標系において点P (20, 0, 100)と点Q (20, 0, 20)を考えます。

カメラとの奥行き方向（Z方向）の距離は、点Pが点Qの5倍となっています。

はじめに腕ならしとして、理想的なピンホールカメラの場合を考えると、K1 = K2 = K3 = 0, P1 = P2 = 0ですので、式(1)(2)(3)(4)より点P, Qの画像への投影点p, qは、画素座標(u, v)にして

p (0.2 f, 0) , q (f, 0)

となります。もし、

点pの画素座標が (1000, 0)ならば、f = 5000

点qの画素座標が (1000, 0)ならば、f = 1000

が得られます。

このように、奥行き方向の距離の倍率と、fの倍率が等しくなるので、ｆは「奥行き方向のスケール」を表すと言えます。この段落の説明で「何ごっこ」をしているのかわかりにくい場合は、こちらの記事の視覚的な説明をご参照ください。

次に、歪みのある現実のカメラ・・・では複雑なので、K1のみ0でない値をもち、K2 = K3 = 0, P1 = P2 = 0の「単純な放射方向歪みをもつカメラ」を考えます。K1 = 0.1とします。このとき、式(1)(2)(3)(4)から求まる2投影点p, qの画素座標は次のようになります。

p: (u, v) = (0.2 * 1.004 * f, 0) = (0.2008 f, 0)

q: (u, v) = (1 * 1.1 * f, 0) = (1.1 f, 0)

従って、

点pの画素座標が (1000, 0)ならば、f = 4980.08

点qの画素座標が (1000, 0)ならば、f = 909.09

と計算されますが、この両者の比 (= 1.1 / 0.2008)は5.48であり、奥行き方向の距離の倍率5と一致しません。

さらに、K1 = 0.1のままでP, Qの位置を再設定してP (10, 0, 100)とQ (10, 0, 20)として同じ計算を行うと、

p: (u, v) = (0.1 * 1.001 * f, 0) = (0.1001 f, 0)

q: (u, v) = (0.5 * 1.025 * f, 0) = (0.5125 f, 0)

となり、u座標が同一の値だとしたときに求まるfの値の比は、5.12と、上記の5.48とはまた異なる値になります。

このように歪みのあるカメラの場合、fの比は奥行き方向の距離の倍率と一致しませんし、画像上の位置にも依存してしまいます。これが表現2は正しくない理由です。

＜グラフで見てみると・・・＞

上記の数値例を一般化して、P (Xp, 0, Zp)とおくと、pのu座標は

u = (Xp/Zp) * (1 + K1 * (Xp/Zp)^2) * f

となり、これより

f = u / ((Xp/Zp) * (1 + K1 * (Xp/Zp)^2))

を得ます。例えばXp = Zp, K1 = -0.25, u = 2736として描いたグラフが下図です。

奥行き方向の距離Zpとfの関係が非線形であり、比例していないことがわかります。

歪みがあるときの内部パラメータfとは何か