数理解析研究所講究録
第 1773 巻 2012 年 68-76

68

Faster Algorithms for Computer Vision
Toshiyuki Masaki and Takahito Kuno
Graduate School of Systems and Information Engineering,
University of Tsukuba, Ibaraki 305-8573, Japan

Abstract

In this paper, we make a modification to Karl and Hartley‘s formulationof problems in computer vision [3, 4], and show the resulting norm minimization problem can be solved by solving a sequence of LP problems.
We also propose an approximation of the norm minimization problem,
which reduces to one LP or SOCP problem.
Key words: Computer vision, multiple view geometry, linear programming, second-order cone programming.

1.

Introduction

In recent years, problems dealt with in computer vision are larger in size andrequire much more computational time to solve than before. For example,
while the traditional triangulation problem has only three variables, the number of variables in the latest structure-and-motion problem amounts to several
hundreds in some cases. Thus, there is a great need for faster algorithms forlarge scale problems in the area of computer vision. In response to this, Karland Hartley [3, 4] developed a framework for solving geometric problems, suchas triangulation, camera-resectioning, homography-estimation and structureand-motion. In their framework, those problems are formulated into annorm minimization problem and solved by solving a sequence of second-ordercone programming (SOCP) problems.
In this paper, we show that, if Karl and Hartley‘s formulation is slightly
modified, the norm $1ni_{I}iiInizati_{oI}i$ problem can be solved by solving a sequenceof linear programming (LP) problems. In addition, we propose an approximation of the norm minimization problem, which reduces to one LP or SOCP
$L_{\infty}$

69

Figure 1: Geometry of a pinhole camera.
problem. Lastly, we report some numerical results, which demonstrate thesuperiority of our approach to Karl and Hartley‘s in both efficiency and accu-racy.

2.

Problems in computer vision

In this section, we first take triangulation as a typical example and illustrate
how it can be formulated into an optimization problem. We also show thatmany other problems in computer vision take the same form. Essential tothose formulations is the pinhole camem model.
Pinhole camera model: The pinhole camera model describes the relationship between the coordinates of a $3D$ point and its projection onto the image
plane of an ideal pinhole camera, where the camera aperture is a pinhole andno lenses are used to focus light. The geometry related to the mapping of apinhole camera is illustrated in Figure 1. Let us denote the object point by$V=(X, Y, Z)^{T}$ in the $3D$ coordinate system with its origin at the camera aperture O. Light emanating from V passes through and projects an inverted
image $v=(x, y)^{T}$ on the image plane, which is parallel to the X-Y plane andlocated at the focal length $f(>0)$ fromin the negative direction of theaxis. Let A $=(0,0, Z)^{T},$ $B=(0,0, -f)^{T}$ and $C=(x, y, -f)^{T}$ . Since thetriangle OAV is similar to the triangle OBC, we have $(x, y)^{T}=(f/Z)(X, Y)^{T}$ ,$O$

$O$

$Z$

70

Figure 2: Thriangulation using two cameras.

or equivalently$\{\begin{array}{l}xy1\end{array}\}=\frac{f}{Z}\{\begin{array}{l}XYZ/f\end{array}\}$in homogeneous coordinates. It should be also noted that the image
invariant under scaling of V. We denote this by$\{\begin{array}{l}xy1\end{array}\}\sim\{\begin{array}{l}XYZ/f\end{array}\}=\{\begin{array}{llll}1 0 0 00 1 0 00 0 1/f 0\end{array}\}\{\begin{array}{l}XYZ1\end{array}\}$,and say that $(x, y, 1)^{T}$ is equivalent, or proportional, to $(X, Y, Z/f)^{T}$ . The
matrix in (1) is called the camem matrix.

$v$is(1)$3\cross 4$

Triangulation:
Triangulation (or reconstruction) is the process of determining the $3D$ coordinates of the object V, given its projection onto two, ormore, images captured by pinhole cameras. In theory, the triangulation problem is quite trivial. Each image of V corresponds to a half-line in the $3D$
space such that all points on the line are projected to . Therefore, V mustlie on the intersection of those lines, and we must be able to calculate its coordinates analytically from a pair of different images. In practice, however,various types of noise, such as geometric noise from lens distortion or interest
point detection error, lead to inaccuracies in the measured image coordinates.
As a result, lines associated with different images of V do not always intersect
in the $3D$ space, as in Figure 2.
$v$

$v$

71
Suppose that $V=(X, Y, Z)^{T}$ is in an arbitrary $3D$ coordinate system, and$M$ .
that there are $M$ imagesof V captured by cameras $i=1,$
Let us denote the ith camera matrix by$v_{i}=(x_{i}, y_{i})^{T}$$\ldots,$$Q_{i}=\{\begin{array}{llll}1 0 0 00 1 0 00 0 1/f_{i} 0\end{array}\}$,

where $f_{i}(>0)$ is the focal length of camera . Note that V is denoted as$B_{n}V+t_{i}$ for some rotation matrix
and a translation vector in the $3D$
coordinate system with the origin at the focal pointof camera . Hence,from (1), we have$i$

$B_{\eta}$$t_{i}$$i$$O_{i}$$\{\begin{array}{l}v_{i}1\end{array}\}\sim Q_{i}\{\begin{array}{ll}R_{n} t_{i}O 1\end{array}\}\{\begin{array}{l}V1\end{array}\}$,$i=1,$$\ldots,$

$M$.

Let
$P_{i}=\{\begin{array}{l}p_{i}^{1}p_{i}^{2}p_{i}^{3}\end{array}\}=Q_{i}\{\begin{array}{ll}R_{i} t_{i}O 1\end{array}\}$,which is referred to as the normalized camem matrix, or simply as the camemmatrix. The coordinates of the imageis then given as$v_{i}=(x_{i}, y_{i})^{T}$$x_{i}= \frac{p_{i}^{1}U}{p_{i}^{3}U}$,$y_{i}= \frac{p_{i}^{2}U}{p_{i}^{3}U}$,

where $U=(V^{T}, 1)^{T}$ , if there is no noise. As mentioned above, however, this is
not the case in practice, and we need to determine the coordinates $(X, Y, Z)$ of
V so as to minimize the $2D$ residual ermr, which is defined as follows in termsofnorm:$L_{p}$$\gamma_{i}=\Vert\frac{p_{i}^{1}U}{p_{i}^{3}U}-x_{i},\frac{p_{i}^{2}U}{p_{i}^{3}U}-y_{i}\Vert_{p}$,$i=1,$$\ldots,$$M$.

If we measure the magnitude ofnorm, then trianin
gulation reduces to an optimization problem with $(X, Y, Z)$ , the first threecomponents of , as the variables:$\gamma=(\gamma_{1}, \ldots, \gamma_{M})^{T}$$L_{q}$$U$

minimizesubject to$\Vert\gamma\Vert_{q}$$\Vert\frac{p_{i}^{1}U}{p_{i}^{3}U}-x_{i},$$\frac{p_{i}^{2}U}{p_{i}^{3}U}-y_{i}\Vert_{p}=\gamma_{i}$, $i=1,$$\ldots$, M.

(2)

72
Table 1: Known and unknown parameters of computer vision problems.
$\frac{problem\neq camems\neq objectsknownunknown}{triangu1ationM1P_{i},v_{i1}V_{1}}$camera-resectioning
structure-and-motion

1$N$$V_{j},$ $v_{1j}$

$M$

$N$

$B_{n},$$v_{ij}$$P_{1}$$t_{i},$$V_{j}$

Other problems: Many other computer vision problems can also be formulated into optimization problems similar to (2) except that the number ofobjects is usually more than one.
$M$ , are given in the $3D$
Suppose $N$ points $V_{j}=(X_{j}, Y_{j}, Z_{j})^{T},$ $j=1,$
space. Let $v_{ij}=(x_{ij}, y_{ij})^{T}$ denote the image of
captured by camera . Also
let$\ldots,$$i$$V_{j}$$\gamma_{ij}=\Vert\frac{p_{i}^{1}U_{j}}{p_{i}^{3}U_{j}}-x_{ij},$

where$U_{j}=(V_{j}^{T}, 1)^{T}$

minimizesubject to$\frac{p_{i}^{2}U_{j}}{p_{i}^{3}U_{j}}-y_{ij}\Vert_{p}$,$i=1,$$\ldots,$$M;j=1,$ $\ldots N$ ,. In this case, the problem is written as follows$\Vert\gamma\Vert_{q}$$\Vert\frac{p_{i}^{1}U_{j}}{p_{i}^{3}U_{j}}-x_{ij},$$\frac{p_{i}^{2}U_{j}}{p_{i}^{3}U_{j}}-y_{ij}\Vert_{p}=\gamma_{ij}$,$\{\begin{array}{ll}i=1, \ldots, M (3)j=1, \ldots N, \end{array}$where. This seems a straightforward extension of (1),but describes various types of problems in computer vision depending on whatparameters are known or unknown. Table 1 shows three major examples; e.g.,
the number of variables amounts to $3(M+N)$ in (3) associated with a structureand-motion problem while that is only three in the case of triangulation.
$\gamma=(\gamma_{11}, \ldots, \gamma_{MN})^{T}$

3.

Solution approaches

In [4], Kahl and Hartley proposed a bisection algorithm to solve (3) with $p=2$
and $q=\infty$ . After illustrating their algorithm, we introduce here a practical
approximation approach to (3).

Bisection approach: In the usual applications of computer vision, we may
assume thatfor all . We can therefore rewrite (3) into the following$p_{i}^{3}U_{j}>0$$i,$ $j$73$\overline{\frac{A1gorithm1Bisectionalgorithmforthecasewherep=2andq=\infty}{Require:aninterva1[\gamma_{\ell},\gamma_{u}]knowntocontaintheoptima1va1uer^{*}and}}$atolerance
repeat$\epsilon>0$.

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

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

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

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

when $p=2$ and $q=\infty$ :minimizesubject to$\gamma$$\Vert p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j},$ $p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}\Vert_{2}\leq\gamma p_{i}^{3}U_{j}$,$\{\begin{array}{l}i=1, \ldots, Mj=1, \ldots N.\end{array}$(4)
As is shown in Algorithm 1, the bisection algorithm solves (4) by checking
repeatedly if the following system of inequalities is feasible for a fixed in a$r$given interval$[\gamma_{\ell}, \gamma_{u}]$:$\Vert p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j},$ $p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}\Vert_{2}\leq\gamma p_{i}^{3}U_{j},$(5)$\{\begin{array}{l}i=1, \ldots, Mj=1, \ldots N.\end{array}$Since the right-hand-sides turn into constants, the feasibility of (5) can bechecked by solving an SOCP problem with (5) as constraints.
As a natural extension of this Kahl and Hartley’s approach, we can apply the bisection algorithm to (3) with another combination of and . For
example, if we choose $p=1$ and $q=\infty$ , then (3) is rewritten as$p$minimizesubject to$q$$\gamma$$|p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j}|+|p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}|\leq\gamma p_{i}^{3}U_{j}$,$\{\begin{array}{l}i=1, \ldots, Mj=1, \ldots N.\end{array}$(6)

Note that the inequality$|p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j}|+|p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}|\leq\gamma p_{i}^{3}U_{j}$74
is equivalent to a set of inequalities$-(p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j})-(p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j})\leq\gamma p_{i}^{3}U_{j}-(U_{j}-xp^{3}U_{j})+(-y_{ij}p^{3})\leq\gamma p^{3}U_{j}(U_{j}U_{ji}(p^{1}U-xp^{3}U)+(U_{j}p_{i}^{3}U)\leq\gamma p^{3}U_{j}p^{2}-y_{ij}p_{i}^{3})\leq\gamma p^{3}U_{j}$$\}$(7).
In other words, the constraints of (6) can be thought of as linear constraints. Hence, Algorith11l 1 can solve (6) by checking theonce the value of is
feasibility of (7) for all $i,j$ , instead of (5). This can be done by solving an LP
problem with (7) for all $i,j$ as constraints. Similarly, we can solve (3) with
$p=\infty$ and $q=\infty$ using the bisection algorithm.
$fix\epsilon_{\text{ノ}}^{\backslash }(1$$\gamma$

Approximation approach: The reason why (3) cannot be solved directly
as an SOCP or LP problem is that the product of two variables appears in theright-hand-sides of the constraints when (3) is rewritten asminimize
subject to$\Vert\gamma\Vert_{q}$$\Vert p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j},$ $p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}\Vert_{p}=\gamma_{ij}p_{i}^{3}U_{j}$,$\{\begin{array}{l}i=1, \ldots, Mj=1, \ldots N.\end{array}$(8)Introducing new variables$\delta_{ij}$and letting$\delta_{ij}=\gamma_{ij}p_{i}^{3}U_{j}$,$i=1,$$\ldots,$$M;j=1,$$\ldots,$$N$,

we haveminimizesubject to$\Vert\delta\Vert_{q}$$\Vert p_{i}^{1}U_{j}-x_{ij}p_{i}^{3}U_{j},$ $p_{i}^{2}U_{j}-y_{ij}p_{i}^{3}U_{j}\Vert_{p}=\delta_{ij}$,$\{\begin{array}{l}i=1, \ldots, Mj=1, \ldots N,\end{array}$(9). It is easy to see that (9) reduces to an SOCP
problem if $p=2$ and $q=\infty$ , and is an LP problem if $p=1,$
and $q=\infty$ .
Furthermore, (9) can serve as a good approximation for (8) (and hence (3)).

where$\delta=(\delta_{11}, \ldots, \delta_{MN})^{T}$$\infty$

Proposition 1 If the optimal value of (8) vanishes, then any optimal solutionof (8), together with $\delta=O$ , solves (9) to optimality. Conversely, if the optimalvalue of (9) vanishes, then any optimal solution of (9), together with $\gamma=O$ ,solves (8) to optimality.
$\square$75
Table 2: CPU time (in seconds).
$\frac{p\backslash q\infty 12}{bisection\infty 0.0660.1170.729}$

approximation

0.003 0.003 0.075
0.009 0.004 0.080$\infty$

1

Table 3: $3D$ residual error (innorm).

$L_{2}$$\frac{p\backslash q\infty 12}{bisection\infty 0.0170.0180.016}$

approximation

0.036 0.036 0.035
0.014 0.014 0.013$\infty$1

Let us discuss the case where $q=\infty$ a little more closely. Let
anddenote the values offor optimal solutions of (8) and(9), respectively, and define$(\overline{\gamma}_{ij}, \overline{p_{i}^{3}U_{j}})$$(\overline{\delta}_{ij},\overline{p_{i}^{3}U_{j}})\wedge$$(\gamma_{ij}, p_{i}^{3}U_{j})$へ$\gamma$where$=\delta/p_{k}^{3}U_{l}$$\overline{\gamma}=\overline{\gamma}_{k\ell}\in\max_{i,j}\{\overline{\gamma}_{ij}\}$and,$\overline{\delta}=\overline{\gamma}\overline{p_{r}^{3}U_{s}}$,へ.

$\hat{\delta}=\hat{\delta}_{rs}\in\max_{i,j}\{\overline{\delta}_{ij}\}$

Proposition 2 When $q=\infty$ , the following relations hold:へ$\overline{\gamma}\leq$$\gamma\leq L\overline{\gamma}$

へ

,

where

$L= \max_{i,j}\{\overline{p_{i}^{3}U_{j}}\}/\min_{i,j}\{\overline{p_{i}^{3}U_{j}}\}$

4.

Numerical results,$\delta\leq\overline{\delta}\leq L\hat{\delta}$.

$\square$Lastly, we report some numerical results obtained with MATLAB (version 7.9,$R2009b)[5]$ .
The problem used as a benchmark is triangulation. The number of cameras
was fixed at $N=25$ , and each camera matrix was of the form:$P_{i}=\{\begin{array}{llll}1 0 0 00 1 0 00 0 1 z_{i}\end{array}\}\{\begin{array}{llll}1 0 0 00 cos\phi_{i} -sin\phi_{i} 00 sin\phi_{i} cos\phi_{i} 00 0 0 1\end{array}\}\{\begin{array}{llll}cos\theta_{i} 0 -sin\theta_{i} 00 1 0 0sin\theta_{i} 0 cos\theta_{i} 00 0 0 1\end{array}\}$,76
where$[$0.0,$\frac{\pi}{6}]$and are uniform random numbers on the intervals [10.0, 100.0] and, respectively, and$z_{i}$$\phi_{i}$$\theta_{i}=\frac{2\pi}{n(i-1)}$.

Each component of the object point V was also generated uniformly at random
in the interval $[-5.0,5.0]$ . Associated SOCP and LP problems were solved
using SeDuMi [7] and GLPK[I], respectively. Table 2 summarizes the average
CPU time in seconds taken to solve ten instances for each. Table 3 shows
the average residual error between actual and calculated values of V in norm.
These results suggest that our approach is superior to Karl and Hartley‘s inboth efficiency and accuracy.
$p,$ $q$$L_{2}$References[1] GLPK(GNU Linear Programming Kit) $http.\cdot//www$ . gnu. $org/software/glpk/$.
[2] Hartley, R.I. and Sturm, P. , Triangulation In COMPUTER ANALYSIS
OF IMAGES AND PATTERNS, pages 190-197, 1995.

[3] Kahl, F. , Multiple view geometry and the-norm. In Int. Conf. Computer Vision,pages 1002-1009, Beijing, China, 2005.
$L_{\infty}$[4] Kahl, F. and Hartley, R.I. , Multiple View Geometry Under the-norm.
In IEEE Tmnsactions on Pattem Analysis and Machine Intelligence, pages
1603-1617, September 2008.
$L_{\infty}$[5] MATLAB $www.mathworks.com/pmducts/matlab/$.

[6] Olsson, C. and Kahl, F. , Generalized Convexity in Multiple View Geometry. In JO URNAL OF MA THEMA TICAL IMA GING AND VISION,
pages 35-51, September 2010.
[7] SeDuMi $http.\cdot//sedumi.ie.lehigh.edu/$.
[8] Seo, Y. and Hartley, R.I. , Sequential
Norm Minimization for Triangulation. In Computer Vision ACCV 2007, vol. 4844, pages 322-331,$L_{\infty}$-2007.