数理解析研究所講究録
第 1829 巻 2013 年 54-63

54

高速な 3 次元再構成のための最適化アプローチ
筑波大学大学院システム情報工学研究科
正木俊行，久野誉人

Toshiyuki Masaki and Takahito Kuno

Graduate School of Systems and Information Engineering,
University of Tsukuba
概要

3 次元再構成問題を最適化手法を用いて解く FKarl らの研究
[4,5] をもとに，より高速な求解のための新たなモデルを提案する．
モデル化の際に残差の定義を見直すことで問題の非線形性を回避し，線形計画問題や最小二乗問題などの解きやすい問題へ帰着させることで，大幅な高速化に成功した．また，計算機実験により提案モデルの優位性を確かめるとともに，両モデルにより得られた解の間に成り立つ関係を示す．

Key words: Computer vision, multiple view geometry, linear
program- ming, second-order cone programming.

1はじめに
3 次元再構成 ( $3D$ -reconstruction) とは，2 次元の画像情報から 3 次元の座標情報を復元する問題のことである．この問題はコンピュータビジョンの分野において古くから扱われているスタンダードな問題である一方，近年の高速な計算機による処理を前提とした大規模な問題の提案もなされており，高速かつ高精度な解法が求められている．

本研究では，幾つかの 3 次元再構成問題を最適化問題へ帰着させるための手法を提案する．最適化問題の多くのクラスには専用の高速なソルバーが数多く存在しており，最適化問題への帰着にはこれらの高速なソルバーを利用できるという利点がある．関連する先行研究としては 2008

55
年の FKarl らによるもの [4, 5] があり，3 次元再構成問題を最適化問題としてモデル化している．この最適化問題は非線形構造を持つために求解が難しく，彼らは二次錐計画問題のソルバーと二分法を利用した反復解法を提案している．しかし，二分法の収束までにかかる反復処理がネックとなり，解が高速に求まらないという問題があった．
我々は問題構造の非線形性を避けるために新たなモデルを定義し，3 次元再構成問題を線形な最適化問題へと帰着させる．これにより問題構造を簡単化し，従来よりも大幅に高速な求解を目指す．さらに目的関数に用いるノルムを変更することで，問題を二次錐計画問題や線形計画問題，二次計画問題，最小二乗問題へと帰着させ，それぞれの問題に帰着した

ときの求解速度や解の精度の違いを数値実験により比較する．また，従来モデルと提案モデルそれぞれの解の間に成り立つ関係式を導出し，特定の条件下で両モデルの解が一致することを示す．

2

ピンホールカメラモデル

図 1: ピンホールカメラモデル

撮影される

$n$

個のオブジェクトの座標を，3 次元の同次座標系におけるベクトルゴ $=[x_{1}^{i}, x_{2}^{i}, x_{3}^{i}, 1]^{T}\in \mathbb{R}^{4}(i=1, \ldots, n)$ で表す．
$m$ 個のカメラに対応するカメラ行列を$P^{j}=\{\begin{array}{l}p_{1}^{j}p_{2}^{j}p_{3}^{j}\end{array}\}\in \mathbb{R}^{3\cross 4}, j=1, \ldots, m$とし，カメラ$P^{j}$によるオブジェクト$x^{i}$のイメージを $y^{ij}=[y_{1}^{ij}, y_{2}^{ij}, 1]^{T}$とする．ピンホールカメラモデルによる撮影の仕組みは図 2 のようになっており，三角形の相似の関係からすぐに$y_{k}^{ij}= \frac{p_{k^{X^{i}}}^{j}}{p_{3^{X^{i}}}^{j}}, k=1,2$が導かれる．

(1)

56
表 1:3 次元再構成問題の未知変数と既知変数$\frac{problem\# ofcameras\neq ofobjectsknownunknown\neq ofvariables}{triangulationm\geq 21P^{j},y^{1j}x^{1}3}$camera-resectioning

1$n\geq 6$$x^{t},y^{l1}$$P^{1}$

11known-rotation-problem

$m$

$n$

$R^{j},y^{ij}$$x^{i},$ $t^{j}$$3(n-1)+3m-1$

3

本研究であつかう問題
本研究であつかう 3 つの 3 次元再構成問題の未知変数と既知変数をまとめたものを表 1 に示す．ただし，回転行列$t^{j}\in \mathbb{R}^{3}$はそれぞれカメラ行列$P^{j}$$R^{j}\in \mathbb{R}^{3\cross 3}$と並進ベクトルの向きと座標に対応するパラメタで，$P^{j}=[R^{j}|t^{j}]$と書けるものとする．

triangulation は，カメラの位置$y^{1j}$$R^{j}$, 方向$t^{j}$とオブジエクトのイメージが与えられた状態で，オブジェクトの 3 次元座標である．camera-resectioning は，オブジエクトの座標$x^{1}$$x^{i}$を推定する問題とそのイメージが与えられた状態で，カメラの位置 , 方向 を推定する．knownとイメージが与えられた状態
rotation-problem は，カメラの方向で，カメラの位置 およびオブジェクトの座標 を推定する．それぞれ$R^{1}$$y^{i1}$$t^{1}$$R^{j}$$y^{ij}$$x^{i}$$t^{j}$の変数の数に着目すると，triangulation と camera-resectioning は $m,$$n$にかかわらず変数の数が決まっているのに対し，known-rotation-problem
は $m,$4$n$の増加にともなって変数の数も増加することがわかる．

従来モデル
いま，$x^{i},$$R^{j},t^{j}$のいずれかに未知数が含まれるとし，等式 (1) よりノルムを用いて残差$\gamma^{ij}$$L_{p}$を定義する．

$\gamma^{ij} =\Vert r_{1}^{ij}, r_{2}^{ij}\Vert_{p},$$r_{k}^{ij} =y_{k}^{ij}-\underline{p_{k^{X^{i}}}^{j}} k=1,2$(2)’$p_{3^{X^{i}}}^{j}$

この残差 (2) を用いて，3 節で述べた 3 次元再構成問題を次のような最57
適化問題として定式化する．
$\min.

\Vert\gamma^{11}, \ldots, \gamma^{nm}\Vert_{q}$s.t.

$\gamma^{ij}=\Vert r_{1}^{ij}, r_{2}^{ij}\Vert_{p}$$i=1,$$\ldots,$(3)$n$$r_{k}^{ij}=y_{k}^{ij}-\underline{p_{k^{X^{i}}}^{j}} k=1,2 j=1, \ldots, m$’$p_{3^{X^{i}}}^{j}$

たとえばカメラの位置$R^{j}$

数として与え，未知変数, 方向$x^{1}$$t^{j}$

およびオブジェクトのイメージ$y^{1j}$を定について問題 (3) を解くことで，triangulation問題におけるオブジェクトの推定座標$x^{1}$を得るのである．

先行研究では $p=2,$ $q=\infty$ のときの問題 (3) の解法を与えている．
$p=2,$ $q=\infty$ のとき，問題 (3) は次の等価な問題に書き換えられる．
$\min.$

s.t.

$\gamma$$\Vert y_{1}^{ij}p_{3}^{j}x^{i}-p_{1}^{j}x^{i},$ $y_{2}^{ij}p_{3}^{j}x^{i}-p_{2}^{j}x^{i}\Vert_{2}\leq\gamma p_{3}^{j}x^{i},$$\{\begin{array}{l}i=1, \ldots, nj=1, \ldots, m\end{array}$(4)(ただし，全てのオブジェクが全てのカメラの前面に位置している，すな
わち

残差が$P_{3^{X^{i}}}^{j}>0$$\gamma$が未知変数であることから，この問題の制約式の右辺は未知変数$R^{j},t^{j}$$x^{i},$であると仮定する)のいずれであっても非線形である．制約条件が非線形の問題を解くことは困難であるため，残差 を固定し，問題 (4) の制約条件の実$\gamma$行可能性のみに着目する．Algorithm 1 に示すように，制約条件の実行可を二分法を用いて探索することで，問題 (4) の解が得られる．このとき，残差 が定数であれば問題 (4) は二次錐計画問題
能性を満たす最小の$\gamma$$\gamma$(Second Order Cone Programming, SOCP) と呼ばれる最適化問題になる
ことから，Algorithm 1 の実装には既存の二次錐計画ソルバーを利用することができる．二次錐計画ソルバーとしては，SeDuMi[8] などがある．

5提案モデル
従来の解法では，3 次元再構成問題を二分法と最適化ソルバーの組み合わせにより解いていた．しかし，二分法が収束するまで二次錐ソルバーを繰り返し呼び出す必要があることから，問題サイズが大きくなるにっれて求解に時間がかかってしまうという問題があった．問題を解くために二分法が必要なのは，問題 3 の制約条件が非線形であることに原因がある．

58$\overline{\frac{A1gorithm1 二分法による問題 (3) の解法 (p=2,q=\infty)}{Require:aninterva1[\gamma_{\ell},\gamma_{u}]knowntocontaintheoptima1valuer^{*}and}}$a tolerance$\epsilon>0.$

repeat$\gammaarrow(\gamma_{\ell}+\gamma_{u})/2$;

check if (4) is feasible or not, by solving an associated SOCP problem;
if (4) is feasible then$\gamma_{u}arrow\gamma$

else$\gamma_{\ell}arrow\gamma$

end ifuntil$\gamma_{u}-\gamma_{\ell}\leq\epsilon$;

そこで我々は，先行研究で用いられていた残差 (2) の定義を見直し，3次元再構成問題を線形制約を持つ最適化問題へと帰着させることを考える．式 (1) を変形して得られる次の式$y_{k}^{ij}p_{3}^{j}x^{i}=p_{k^{X^{i}}}^{j},$から，残差$\delta^{ij}$(5)$k=1,2$を次のように定義する．
$\delta^{ij} =\Vert d_{1}^{ij}, d_{2}^{ij}\Vert_{p},$(6)$d_{k}^{ij} =y_{k}^{ij}p_{3}^{j}x^{i}-p_{k}^{j}x^{i}, k=1,2$

3 次元再構成問題においてはイメージ$y^{ij}$は既知であるので，$P^{j_{X^{j}}}$の各列が線形である限り曜は線形であることに注意されたい．この新たな残差(6) を用いれば，3 次元再構成問題を次の最適化問題に帰着させることができる．
$\min.

\Vert\delta^{11}, \ldots, \delta^{nm}\Vert_{q}$$s.t.$$\delta^{ij}=\Vert\dot{d}_{1}^{j},$$d_{2}^{ij}\Vert_{p}$$i=1,$$\ldots,$(7)$n$$d_{k}^{ij}=y_{k}^{ij}p_{3}^{j}x^{i}-p_{k}^{j}x^{i}, k=1,2 j=1, \ldots, m$

問題 (3) と問題 (7) の違いは，2 本目の制約式の線形性にある．項$p_{1}^{j}x^{i},$$p_{2}^{j}x^{i},$は，triangulation, camera-resectioning, known-rotation-problem のすべての問題で，未知変数について線形である．known-rotation-problem では$P^{j}$の並進ベクトル$t^{j}$が未知変数なので一見線形でなく見えるが，式を展開すれば$P^{j}x^{i}=[R^{j}|t^{j}]x^{i}=R^{j}\{\begin{array}{l}x_{1}^{j}x_{2}^{j}x_{3}^{j}\end{array}\}+t^{j}$$p_{3^{X^{i}}}^{j}$59
表 2: 問題 (3) の解法

表 3: 問題 (7) の解法$\frac{p\backslash q\infty}{\infty 二分法\ 線形計画}$$\frac{p\backslash q\infty 12}{\infty 線形計画線形計画二次計}$

1

二分法 &線形計画

1

線形計画線形計画二次計画

2

二分法 &二次錐計画

2

二次錐計画二次錐計画最小二乗となることから線形であることがわかる．よって (7) の制約式
$d_{k}^{ij}=y_{k}^{ij}p_{3}^{j}x^{i}-p_{k}^{j}x^{i}, k=1,2$は線形である．ところが，(3) の制約式
$r_{k}^{ij}=y_{k}^{ij}-\underline{p_{k}^{j}x^{i}} k=1,2$’$p_{3^{X^{i}}}^{j}$は，未知変数を含む式$p_{k}^{j}x^{i}$と$p_{3^{X^{i}}}^{j}$の分数形式であることから線形ではない．

この違いは問題の複雑さに大きく関係しており，たとえば問題 (7) において $p=q=1$とした場合には次のような線形計画問題 (Linear Program-ming, $LP$ ) に帰着できる．
$\min.

\sum_{i,j,k}|\delta_{k}^{ij}|$st.

$d_{k}^{ij}=y_{k}^{ij}p_{3}^{j}x^{i}-p_{k^{X^{i}}}^{j},$$k=1,2\{$$i=1,$$\ldots,$$n$(8)$j=1, \ldots, m$

問題 (8) は二分法を用いることなく，線形計画ソルバーのみで解くことができる．他にも $p=1,$ $q=2$ とすれば二次計画問題 (Quadratic Programming,$QP$ )$p=q=2$ とすれば最小二乗問題 (Least Square Problem,として，LSP) として解くことができる．
問題 (3) と問題 (7) の解法をまとめたものが表 2,$q=\infty$

表 3 である．表 2 がの場合のみ示してあるのは，二分法により変数を固定する関係上必ず $q=\infty$ である必要があり，$q=1,2$ のときには解けないためである．

6

ふたつの問題の関係問題 (3) と (7) により得られた解の間には次の定理が成り立っ．

60定理

問題 (3), (7) の一方の問題の最適解は，もう一方の問題の実行可能解である．また，それぞれの最適解をれたときのそれぞれの問題の目的関数値をとし，解$\overline{\Omega},\tilde{\Omega}$$\Gamma(\Omega),$$\triangle(\Omega)$$\triangle(\tilde{\Omega}) \leq \triangle(\overline{\Omega}) \leq L\cdot\triangle(\tilde{\Omega})$$\Omega$が与えらとするとき，,$\Gamma(\overline{\Omega}) \leq \Gamma(\tilde{\Omega}) \leq L\cdot\Gamma(\overline{\Omega})$が成り立つ．ただし， は$L$$L= \frac{\max_{i,j}p_{3}^{j}x^{i}(\overline{\Omega})}{\min_{i,j}p_{3}^{j}x^{i}(\tilde{\Omega})}$であたえられる．

また，以下のいずれかを満たすとき，二つの問題の最適解は一致する．

1.

$\Gamma(\overline{\Omega})=0$2.

$\triangle(\tilde{\Omega})=0$

3. $L=1$

この定理から，一方の問題で得られた最適解をもう一方の問題の近似解と見なせ，その精度の上界がた，どちらかの問題に残差$O$$L$で与えられると考えることができる．まの最適解が存在する場合には，二つの問題の最適解が必ず一致することに注目されたい．

7

数値実験
known-rotation-problem をベンチマークとして数値実験を行った．実

装には MATLAB version $7.9, R2009b)$ $[7]$ 用いている．SOCP ソルバー$($$LP$ ソルバーとして GLPK[ ] およびその MATLABとして SeDuMi[8] を，$I$$QP$ ソルバーとして MATLAB optiインターフェース GLPK MEX[2] を，mization toolbox から quadprog 関数を，LSM ソルバーとして MATLAB
標準の mldivide 関数を利用している．

数値実験の手順を以下に示す．

1. 乱数を用いてオブジェクト座標のとき，$p_{3^{X^{i}}}^{j}>0$$x^{i}$, カメラ行列$P^{j}$を生成する．こを満たすため，カメラが全オブジェクトをぐるりと取り囲むような構造を与えた．

61
表 4: 実験結果 (試行回数:10 回)
残差$\gamma$

pq

221

計算時間 $0$)

問題

誤差$\infty$—分法 &二次錐計画

0.000706

11.551$\infty$

二次錐計画$\infty$

線形計画

0.000761
0.001568
0.000813
0.000684

0.479
0.196
0. $399$
0.011$\delta$

11

線形計画

22

最小二乗

2. イメージ$y^{ij}$

速度$((_{\mathbb{E}}^{g_{i}})$

lx
$24x$

$59x$
$29x$

$1050x$を次式より求める．

$y^{ij}=[\frac{p_{1}^{j}x^{i}}{p_{3^{X^{\dot{x}}}}^{j}}, \frac{p_{2}^{j}x^{i}}{p_{3^{X^{i}}}^{j}}]^{T}$

3. カメラの回転行列$R^{j}$とイメージ$y^{ij}$に適当な誤差を与えたものをパラメータとして，known-rotation-problem を解く．

4. 得られた最適解$\overline{x}^{i},\overline{t}^{j}$と真値との ($L_{2}$

ノルムによる) 平均誤差$\frac{1}{n}\sum_{i=1}^{n}\Vert\overline{x}^{i}-x^{i}\Vert_{2}+\frac{1}{m}\sum_{j=1}^{m}\Vert\overline{t}^{j}-t^{j}\Vert_{2}$を評価する．
オブジェクト数 $n=30$ , カメラ数 $m=10$ としたときの結果を表 4 に示す．この結果から，提案したモデル (残差 ) が既存の手法にくらべ精度．
$\delta$

速度ともに優れたパフオーマンスを示していることがわかる．特に最小二乗問題として解いた $p=q=2$ の提案モデルでは，従来手法の 1000 倍という速度が出ている．これだけの高速化が実現できた理由は，一つには二分法のための反復の有無によるところが大きい．もうーつの理由としては，問題そのものがかなり解きやすい形になっていることが挙げられる．特に $p=q=2$ の提案モデルでは制約条件が存在しないため，提案モデルを用いた他のモデルに比べても高速である．

また図 7 は，カメラ数 $m=40$ , オブジェクト数を $n=5,10,50,100,$$\ldots,$300 と変化させたときの，CPU 時間の変化をプロットしたものである．提案手法とは，表 4 で最も高速であった $p=q=2$ のモデルを指す．オブジエクト数が増加しても，提案手法は従来手法よりも安定して高速であることがわかる．最も大きなサイズである $n=300,$$m=40$ の問題でも，

62
700
600
$\prime$

血$\underline{}$$\underline{}$

提案手法$\prime$$500-$,-巳従来手法$\underline{\overline{\prime r}}$

$400$

官誕

$\prime\blacksquare’$$\Delta\cup\supset 300--\bullet^{=}\neq-$$-,$
200$100 -\sim’-\bullet.-’$$0$$0 50 100 150 200 250 300$$n$

図 2: オブジェクト数$n\ovalbox{\tt\small REJECT}$

こ対する計算時間の変化 (試行回数: 10 回)

従来手法の 603 秒に対して提案手法は 37 秒となっている．

参考文献$http.\cdot//www$Programming

Linear[1] GLPK(GNUKit). gnu. $org/soflware/glpk/.$[2] GLPKMEX-A Matlab MEX Interface For The GLPK Library

glpkmex. sourceforge. $net/.$$http.\cdot//$[3] Hartley, R. . and Sturm, P. , Triangulation In COMPUTER ANAL$I$

YSIS OF IMAGES AND PATTERNS, pages 190-197, 1995.
[4] Kahl, F. , Multiple view geometry and the$L_{\infty}$-norm. In Int. Conf.

Computer Vision,pages 1002-1009, Beijing, China, 2005.
[5] Kahl, F. and Hartley, R. . ,$I$

Multiple View Geometry Under the-norm. In IEEE Transactions on Pattem Analysis and MachineIntelligence, pages 1603-1617, September 2008.

$L_{\infty}$[6] Olsson, C. and Kahl, F. , Generalized Convexity in Multiple View

Geometry. In JOURNAL OF MA THEMA TICAL $IMA$ GING AND
VISION, pages 35-51, September 2010.

63[7] MATLAB $www.$ mathworks. $com/products/matlab/.$[8] SeDuMi $http.\cdot//$sedumi. $ie.$ lehigh. $edu/$[9] Seo, Y. and Hartley, R. . , Sequential$I$

Triangulation. In Computer Vision-$L_{\infty}$

Norm Minimization for

ACCV 2007, vol. 4844, pages

322-331, 2007.
[10] SPDA In $http.\cdot//sdpa$ . indsys. chuo-u. $ac.jp/sdpa/$