Search
Duplicate

Computer Vision/ 4. Stereo Systems

1 Introduction

이전 노트에서 장면에 대해 추가 viewpoint를 더하면 해당 장면에 대한 지식을 크게 향상시킬 수 있음을 다루었다. 3d 장면에 관한 정보를 추출하지 않고도 한 이미지 평면의 점을 다른 이미지 평면의 점에 관련지을 수 있도록 epipolar geometry 설정에 초점을 맞추었다. 이번 강의에서는 여러 2d 이미지에서 3d 장면에 관한 정보를 복구하는 방법에 대해 논의한다.

2 Triangulation

여러 view geometry에서 가장 근본적인 문제 중 하나는 triangulation(삼각측량) 문제이다. 이는 3d 점의 위치를 2개 이상의 이미지에 투영된 것으로부터 결정하는 과정이다.
두 view에서의 triangulation 문제에서는 각각 알려진 카메라 intrinsic 파라미터 KKKK'를 가진 2대의 카메라를 갖는다. 또한 이 카메라들의 상대적인 orientation과 offset R,TR, T를 갖는다. 3d에 점 PP가 있다고 가정하자. 이것은 두 카메라의 이미지에서 각각 pppp'에 해당한다. PP의 위치가 현재 알려지지 않지만 이미지에서 pppp'의 정확한 위치를 측정할 수 있다. K,K,R,TK, K', R, T가 알려져 있으므로 카메라 중심 O1,O2O_1, O_2와 이미지 위치 p,pp, p'에 의해 정의되는 두 시선 \ell\ell'을 계산할 수 있다. 따라서 PP\ell\ell'의 교점으로 계산될 수 있다.
이 절차가 간단하고 수학적으로 타당해 보이지만, 실제로는 잘 작동하지 않는다. 왜냐하면 현실에서 관측치 pppp'가 noisy이고 camera calibration 파라미터가 정밀하지 않기 때문에 \ell\ell'의 교차점을 찾는 것이 문제가 있다. 대부분의 경우에서 두 선이 교차하지 않으므로 교점이 전혀 존재하지 않을 수 있다.

2.1 A linear method for triangulation

이 섹션에서 ray 사이의 교차 점의 부재를 해결하는 간단한 linear triangulation 방법을 설명한다. 서로 대응하는 두 이미지의 점이 p=MP=(x,y,1)p = MP = (x, y, 1)p=MP=(x,y,1)p' = M'P = (x', y', 1)로 주어진다. cross product의 정의에 따라 p×(MP)=0p \times (MP) = 0이다. cross product에 의해 생성된 등식을 명시적으로 사용하여 3가지 제약조건을 형성할 수 있다.
x(M3P)(M1P)=0y(M3P)(M2P)=0x(M2P)y(M1P)=0(2.1)\begin{aligned} x(M_3P) - (M_1P) &= 0 \\ y(M_3P) - (M_2P) &= 0 \\ x(M_2P) - y(M_1P) &= 0 \end{aligned} \tag{2.1}
여기서 MiM_i는 행렬 MMii-번째 행이다. 유사한 제약조건을 pp'MM'에 대해 형식화할 수 있다. 두 이미지의 제약조건을 사용하여 AP=0AP = 0 형식의 선형 방정식을 형식화할 수 있다. 여기서
A=[xM3M1yM3M2xM3M1yM3M2](2.2)A = \begin{bmatrix} xM_3 - M_1 \\ yM_3 - M_2 \\ x'M_3' - M_1' \\ y'M_3' - M_2' \end{bmatrix} \tag{2.2}
이 등식은 점 PP의 best linear 추정치를 찾기 위해 SVD를 사용하여 해결될 수 있다. 이 방법의 또 다른 흥미로운 측면은 여러 view에서 triangulating도 잘 다룰 수 있다는 것이다. 그렇게 하기 위해 새로운 view에 의해 추가된 제약조건에 해당하는 행을 AA에 추가하면 된다.
그러나 이 방법은 projective-invariant이 아니기 때문에 projective reconstruction에 적합하지 않다. 예컨대 카메라 행렬 M,MM, M'을 projective transformation MH1,MH1MH^{-1}, M'H^{-1}에 의해 영향 받는 것으로 교체한다고 가정하자. 그러면 선형 방정식의 행렬 AAAH1AH^{-1}이 된다. 그러므로 이전 추정치 AP=0AP= 0의 해 PP는 변환된 문제 (AH1)(HP)=0(AH^{-1})(HP) = 0의 해 HPHP에 해당한다. SVD는 P=1\|P\| = 1의 제약조건을 해결하지만 이것은 projective transformation HH에 대해 불변이 아니다. 그러므로 이 방법은 간단하지만 종종 triangulation 문제의 최적 해가 아니다.

2.2 A nonlinear method for triangulation

대신 현실 세계 시나리오에 대한 triangulation 문제는 종종 수학적으로 최소화 문제를 해결하는 것으로 특성화 된다.
minP^MP^p2+MP^p2(2.3)\min_{\hat{P}}\|M\hat{P} - p\|^2 + \|M'\hat{P}-p'\|^2 \tag{2.3}
위의 방정식에서 우리는 두 이미지에서 P^\hat{P}의 reprojection error의 best least-square을 찾음으로써 PP를 가장 잘 근사하는 3d 공간의 P^\hat{P}를 찾는다. 이미지에서 3d 점에 대한 reprojection error는 그 점을 이미지에 projection 한 것과 이미지 평면에서 관측된 대응점 사이의 거리이다. 그림 2의 예시에서 MM은 3d 공간에서 이미지 1에 대한 projective transformation이므로 P^\hat{P}를 이미지 1에 project한 점은 MP^M\hat{P}이다. 이미지 1에서 P^\hat{P}에 대응하는 관측치는 pp이다. 따라서 이미지 1에 대한 reprojection error는 거리 MP^p\|M\hat{P} - p\|이다. 방정식 2.3에서 발견한 전체 reprojection error는 모든 이미지에 걸친 reprojection error의 합이다. 2개 이상의 이미지가 있는 경우에 목적 함수에 거리 항을 더 추가하기만 하면 된다.
minP^iMP^ipi2(2.4)\min_{\hat{P}} \sum_i \|M\hat{P}_i - p_i\|^2 \tag{2.4}
실제에서 문제에 대한 좋은 근사치를 얻을 수 있는 다양한 정교한 최적화 방법이 존재한다. 그러나 이 수업의 범위에서 우리는 이러한 기법들 중 하나인 nonlinear least squares에 대한 Gauss-Newton algorithm에만 초점을 맞춘다. 일반적인 nonlinear least squares 문제는 다음을 최소화하는 xRnx \in \mathbb{R}^n을 찾는 것이다.
r(x)2=i=1mri(x)2(2.5)\|r(x)\|^2 = \sum_{i=1}^m r_i(x)^2 \tag{2.5}
여기서 rr은 어떤 함수 ff, 입력 xx와 관찰 yy에 대해 r(x)=f(x)yr(x) = f(x) - y와 같은 임의의 residual function r:RnRmr : \mathbb{R}^n \to \mathbb{R}^m이다. 함수 ff가 선형일 때 nonlinear least squares 문제는 regular linear least square 문제로 축소된다. 그러나 일반적으로 우리의 카메라 행렬이 affine이 아님을 떠올려라. 왜냐하면 이미지 평면에 대한 projection이 종종 homogeneous 좌표로 나누기 때문이다. 이미지에 대한 projection은 일반적으로 비선형이다.
eie_i2×12 \times 1 벡터 ei=MiP^pie_i = M_i\hat{P} - p_i(원문에는 ei=MP^ipie_i = M\hat{P}_i - p_i로 되어 있지만 오타인 듯) 설정하면 최적화 문제를 다음과 같이 재형식화 할 수 있다.
minP^iei(P^)2(2.6)\min_{\hat{P}} \sum_i e_i(\hat{P})^2 \tag{2.6}
이것은 완벽하게 nonlinear least squares 문제로 표현될 수 있다.
이 노트에서 Gauss-Newton algorithm을 사용하여 이 nonlinear least squares 문제에 대한 근사 해를 찾는 방법을 다룬다. 우선 이전 선형 방법으로 계산될 수 있는 3d 점 P^\hat{P}에 대한 합리적인 추정치가 있다고 가정한다. Guass-Newton algorithm의 핵심 통찰은 reprojection error를 최소화하는 더 나은 추정치 방향으로 현재 추정치를 보정하여 업데이트하는 것이다. 각 단계에서 추정치 P^\hat{P}를 어떤 δP:P^=P^+δP\delta_P : \hat{P} = \hat{P} + \delta_P만큼 업데이트한다.
그러나 업데이트 파라미터 δP\delta_P를 어떻게 선택할 수 있을까? Gauss-Newton 알고리즘의 핵심 통찰은 현재 추정치 P^\hat{P} 근처에서 residual 함수를 선형화하는 것이다. 우리 문제의 경우 이것은 점 PP의 residual error ee를 다음과 같이 생각할 수 있다는 것을 의미한다.
e(P^+δP)e(P^)+ePδP(2.7)e(\hat{P} + \delta_P) \approx e(\hat{P}) + {\partial e \over \partial P} \delta_P \tag{2.7}
결과적으로 최소화 문제를 다음으로 변환한다.
minδPePδP(e(P^))2(2.8)\min_{\delta_P} \|{\partial e \over \partial P} \delta_P - (-e(\hat{P}))\|^2 \tag{2.8}
이와 같이 residual을 형식화하면 표준 linear least squares 문제의 형식을 취하는 것을 볼 수 있다. NN개 이미지에 대한 triangulation 문제에 대해 linear least squares 해는 다음이 된다.
δP=(JJ)1Je(2.9)\delta_P = -(J^\top J)^{-1}J^\top e \tag{2.9}
여기서
e=[e1eN]=[p1M1P^pnMnP^](2.10)e = \begin{bmatrix} e_1 \\ \vdots \\ e_N \end{bmatrix} = \begin{bmatrix} p_1 - M_1\hat{P} \\ \vdots \\ p_n - M_n \hat{P} \end{bmatrix} \tag{2.10}
그리고
J=[e1P^1e1P^2e1P^3eNP^1eNP^2eNP^3](2.11)J = \begin{bmatrix} {\partial e_1 \over \partial \hat{P}_1} & {\partial e_1 \over \partial \hat{P}_2} & {\partial e_1 \over \partial \hat{P}_3} \\ & \vdots \\ {\partial e_N \over \partial \hat{P}_1} & {\partial e_N \over \partial \hat{P}_2} & {\partial e_N \over \partial \hat{P}_3} \end{bmatrix} \tag{2.11}
특정한 이미지의 residual error 벡터 eie_i는 이미지 평면에서 두 차원이 존재하기 때문에 2×12 \times 1 벡터이다. 따라서 triangulation의 가장 단순한 두 개의 카메라 경우 (N=2N=2), residual vector ee2N×1=4×12N \times 1 = 4 \times 1 벡터가 되고 야코비안 JJ2N×3=4×32 N \times 3 = 4 \times 3 행렬이 된다. 추가 이미지는 ee 벡터와 JJ 행렬에 해당하는 행을 추가하여 처리되므로 이 방법이 multiple view를 seamlessly 다룰 수 있음을 알 수 있다. 업데이트 δP\delta_P를 계산한 후에 고정된 단계 수만큼 또는 수치적으로 수렴할 때까지 이 과정을 반복하면 된다. Gauss-Newton 알고리즘의 가장 중요한 속성 중 하나는 residual function가 우리의 추정치 근처에서 선형이라는 가정이 수렴을 보장하지 않는다는 것이다. 따라서 실제로는 추정치에 대한 업데이트 횟수에 상한을 두는 것이 항상 유용하다.

3 Affine structure from motion

이전 섹션의 마지막 부분에서 장면의 두 view를 너머 3d 장면에 관한 정보를 얻을 수 있음을 시사했다. 이제 두 대의 카메라의 geometry를 multiple 카메라로 확장하는 방법을 탐구한다. multiple view에서 관찰 점을 결합하여 structure from motion으로 알려진대로 장면의 3d 구조와 카메라의 파라미터 모두를 동시에 결정할 수 있다.
여기서 structure from motion 문제를 형식적으로 유도한다. mm개의 카메라에 대해 intrinsic과 extrinsic 파라미터 모두를 인코딩하는 카메라 변환 MiM_i가 있다고 가정하자. XjX_j를 장면의 nn개 3d 포인트 중 하나라고 하자. 각 3d point는 multiple 카메라에서 볼 수 있으며, xijx_{ij}XjX_j를 projective transformation MiM_i를 사용하여 카메라 ii의 이미지로 projection한 것이다. structure from motion의 목표는 모든 관찰 xijx_{ij}에서 장면의 structure(nn개 3d 점 XjX_j)와 카메라의 motion(mm개 projection 행렬 MiM_i) 모두를 복구하는 것이다.

3.1 The affine structure from motion problem

일반적인 structure from motion 문제를 다루기 전에 우선 카메라가 affine 또는 weak perspective라고 가정하는 더 간단한 문제부터 시작한다. 궁극적으로 perspective scaling 연산이 없어지면 이 문제에 대한 수학적 유도가 더 쉬워진다.
이전에 perspective와 weak perspective 경우에 대해 위의 등식을 유도했다. full perspective 모델에서 카메라 행렬이 다음과 같이 정의되는 것을 떠올려라.
M=[Abv1](3.1)M = \begin{bmatrix} A & b \\ v & 1 \end{bmatrix} \tag{3.1}
여기서 vv는 어떤 non-zero 1×31 \times 3 벡터이다. 반면에 weak perspective model의 경우 v=0v = 0이다. 이 속성은 MXMX의 homogeneous 좌표를 1과 동일하게 만든다.
x=MX=[m1m20001][X1X2X31]=[m1Xm2X1](3.2)x = MX = \begin{bmatrix} m_1 \\ m_2 \\ 0 \quad 0 \quad 0 \quad 1 \end{bmatrix} \begin{bmatrix} X_1 \\ X_2 \\ X_3 \\ 1 \end{bmatrix} = \begin{bmatrix} m_1X \\ m_2X \\ 1 \end{bmatrix} \tag{3.2}
결과적으로 homogeneous에서 유클리드 좌표로 이동하면서 projective transformation의 비선형성은 사라지고 weak perspective transformation은 더 magnifier(확대기) 역할을 한다. 우리는 projection을 다음과 같이 더 간결하게 나타낼 수 있다.
[m1Xm2X]=[Ab]X=AX+b(3.3)\begin{bmatrix} m_1X \\ m_2X \end{bmatrix} = \begin{bmatrix} A & b \end{bmatrix}X = AX + b\tag{3.3}
그리고 모든 카메라 행렬을 Maffine=[Ab]M_\text{affine} = [A \quad b] 형식으로 나타낼 수 있다. 따라서 우리는 이제 affine 카메라 모델을 사용하여 3d에서 점 XjX_j과 각 affine 카메라에서의 대응 관측치간의 관계(예컨대 카메라 ii에서 xijx_{ij})를 표현할 수 있다.
structure from motion 문제로 돌아가서, 우리는 mnmn개 관찰에서 mm개 행렬 MiM_inn개 월드 좌표 벡터 XjX_j, 총 8m+3n8m + 3n개의 미지수를 추정해야 한다. 각 관측치는 카메라당 2개의 제약을 생성하므로 8m+3n8m + 3n개의 미지수에 대해 2mn2mn 방정식이 존재한다. 우리는 이 방정식을 사용하여 각 이미지에서 대응점이 필요한 최소 개수의 하한을 알 수 있다. 예컨대 m=2m=2개의 카메라를 가지면, 3d에서 최소 n=16n=16개의 점이 필요하다. 그러나 각 이미지에서 라벨링된 대응점을 충분하면, 이 문제를 어떻게 해결할 수 있는가?

3.2 The Tomasi and Kanade factorization method

이 파트에서 우리는 affine structure from motion 문제를 해결하기 위한 Tomasi와 Kanade의 fatorization 방법을 개략적으로 설명한다. 이 방법은 2가지 주요 단계로 구성된다. data centering step과 actual factorization step.
data centering 단계에서 시작하자. 이 단계에서 주요 아이디어는 원점에 데이터를 center하는 것이다. 이것을 위해 각 이미지 ii에 대해 각 이미지 점 xijx_{ij}에 대해 그들의 centroid xˉi\bar{x}_i를 빼서 새로운 좌표 x^ij\hat{x}_{ij}를 재정의한다.
x^ij=xijxˉi=xij1nj=1nxij(3.4)\hat{x}_{ij} = x_{ij} - \bar{x}_i = x_{ij} - {1\over n} \sum_{j=1}^n x_{ij} \tag{3.4}
affine structure from motion 문제가 이미지 포인트 xijx_{ij}, 카메라 행렬 변수 AiA_ibib_i 그리고 3d 점 XjX_j 사이의 관계가 다음과 같음을 떠올려라.
xij=AiXj+bi(3.5)x_{ij} = A_i X_j + b_i \tag{3.5}
이 centering 단계 이후에, 방정식 3.4의 centered 이미지 점 x^ij\hat{x}_{ij}의 정의와 방정식 3.5에서 affine 표현식을 결합할 수 있다.
x^ij=xij1nk=1nxik=AiXj1nk=1nAiXk=Ai(Xj1nk=1nXk)=Ai(XjXˉ)=AiX^j(3.6)\begin{aligned} \hat{x}_{ij} &= x_{ij} - {1\over n} \sum_{k=1}^n x_{ik} \\ &= A_iX_j - {1\over n}\sum_{k=1}^n A_iX_k \\ &= A_i(X_j - {1\over n}\sum_{k=1}^n X_k) \\ &= A_i(X_j-\bar{X}) \\ &= A_i \hat{X}_j \end{aligned} \tag{3.6}
방정식 3.6에서 볼 수 있듯이, world reference system의 원점을 centroid Xˉ\bar{X}로 translate 하면 이미지 점의 centered 좌표 x^ij\hat{x}_{ij}와 3d 점의 centered 좌표 X^ij\hat{X}_{ij}는 단일 2×32 \times 3 행렬 AiA_i에 의해서만 연관된다. 궁극적으로 factorization 방법의 centering 단계를 통해 3d 구조와 multiple 이미지에서의 관측 점 간의 관계를 간편한 행렬 곱 표현으로 만들 수 있다.
그러나 행렬 곱 x^ij=AiX^j\hat{x}_{ij} = A_i \hat{X}_j에서 등식의 좌변의 값에만 접근할 수 있음에 유의하라. 따라서 어떻게든 motion 행렬 AiA_i와 구조 XjX_j를 인수분해 해야 한다. mm개 카메라에 대한 모든 관찰을 사용하여 nn개의 관찰으로 이루어진 measurement matrix DD를 구축할 수 있다. (각 x^ij\hat{x}_{ij} 항목이 2×12\times1 벡터임을 기억하라)
D=[x^11x^12...x^1nx^21x^22...x^2nx^m1x^m2...x^mn](3.7)D = \begin{bmatrix} \hat{x}_{11} & \hat{x}_{12} & ... & \hat{x}_{1n} \\ \hat{x}_{21} & \hat{x}_{22} & ... & \hat{x}_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ \hat{x}_{m1} & \hat{x}_{m2} & ... & \hat{x}_{mn} \end{bmatrix} \tag{3.7}
이제 affine 가정 때문에 DD2m×32m \times 3 motion matrix MM(이것은 카메라 행렬 A1,...,AmA_1, ..., A_m로 구성됨)과 3×n3 \times n structure matrix SS(이것은 3d 점 X1,...,XnX_1, ..., X_n으로 구성됨)의 곱으로 표현될 수 있음을 떠올려라. 중요한 사실은 DD가 최대 차원이 33인 두 행렬의 곱이므로 rank(D)=3\text{rank}(D) = 3을 사용할 수 있다는 것이다.
DDMMSS로 인수분해하여, SVD D=UΣVD = U\Sigma V^\top를 사용할 수 있다. rank(D)=3\text{rank}(D) = 3이므로 Σ\Sigma에는 3개의 non-zero singular value σ1,σ2,σ3\sigma_1, \sigma_2, \sigma_3만 존재한다. 따라서 표현식을 더 축소하여 다음의 분해를 얻을 수 있다.
D=UΣV=[u1...un][σ1000...00σ200...000σ30...00000...0][v1vn]=[u1u2u3][σ1000σ2000σ3][v1v2v3]=U3Σ3V3(3.8)\begin{aligned} D &= U\Sigma V^\top \\ &= \begin{bmatrix} u_1 & ... & u_n \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 & 0 &0 & ... & 0 \\ 0 & \sigma_2 & 0 & 0 & ... & 0 \\ 0 & 0 & \sigma_3 & 0 & ... & 0 \\ &&&&\ddots \\ 0 & 0 & 0 & 0 & ... & 0 \end{bmatrix} \begin{bmatrix} v_1^\top \\ \vdots \\ v_n^\top \end{bmatrix} \\ &= \begin{bmatrix} u_1 & u_2 & u_3 \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 &\sigma_3 \end{bmatrix} \begin{bmatrix} v_1^\top \\ v_2^\top \\ v_3^\top \end{bmatrix} \\ &= U_3\Sigma_3V_3^\top \end{aligned} \tag{3.8}
이 decomposition에서 Σ3\Sigma_3은 non-zero singular 값으로 이루어진 대각 행렬이며, U3U_3V3V_3^\top은 각각 UUVV^\top의 3개 열과 행을 각각 취하여 얻어진다. 불행히도 실제에서 측정 noise와 affine 카메라 근사치 때문에 rank(D)>3\text{rank}(D) > 3이다. 그러나 rank(D)>3\text{rank}(D) > 3일 때 U3W3V3U_3W_3V_3^\top이 여전히 Frobenius norm의 관점에서 MSMS의 가장 좋은 rank-3 근사이다.
자세히 보면, 행렬 곱 Σ3V3\Sigma_3V_3^\top3×n3\times n 행렬로 structure 행렬 SS와 정확히 동일한 크기이고, 유사하게 UU2m×32m \times 3 행렬로 motion 행렬 MM과 동일한 크기이다. SVD decomposition의 성분을 MMSS에 연관하는 방법은 affine structure from motion 문제에 대한 물리적, 기하학적으로 타당한 해를 제공하지만 이 선택이 유일한 해는 아니다. 예컨대 motion 행렬을 M=U3Σ3M = U_3\Sigma_3으로 structure 행렬을 S=V3S = V_3^\top로 설정해도 관찰 행렬 DD는 동일하다. 그렇다면 어떤 factorization을 선택해야 하는가? Tomasi와 Kanade는 그들의 논문에서 견고한 factorization의 선택은 M=U3Σ3M = U_3\sqrt{\Sigma_3}S=Σ3V3S = \sqrt{\Sigma_3}V_3^\top라고 결론 내렸다.

3.3 Ambiguity in reconstruction

그럼에도 불구하고, D=MSD = MS의 factorization 선택에 본질적인 모호성을 발견할 수 있다. 이는 임의의 가역인 3×33 \times 3 행렬 AA를 decomposition에 추가할 수 있기 때문이다.
D=MAA1S=(MA)(A1S)(3.9)D = MAA^{-1}S = (MA)(A^{-1}S) \tag{3.9}
이것은 motion MM에서 얻은 카메라 행렬과 structure SS 에서 얻은 3d 점이 공통 행렬 AA에 의해 곱해진 것으로 결정됨을 뜻한다. 그러므로 우리 해는 under-determined 이고 이 affine ambiguity를 해결하기 위해 추가 제약조건이 필요하다. reconstruction이 affine ambiguity를 가질 때, 평행성은 보존되지만 metric scale은 알 수 없음을 의미한다.
reconstruction을 위한 또 다른 중요한 ambiguities 클래스는 similarity ambiguity이다. 이것은 reconstruction이 similarity transform(rotation, translation, scaling)까지만 정확할 때 발생한다. similarity ambiguity만 갖는 reconstruction을 metric reconstruction이라 한다. 이 ambiguity는 카메라가 intrinsically calibrated이더라도 존재한다. 좋은 소식은 calibrated 카메라의 경우 similarity ambiguity는 유일한 ambiguity이라는 것이다.
이미지에서 장면의 절대 scale을 복구하는 방법이 없다는 사실은 꽤 직관적이다. 추가 가정(예: 그림에서 집의 높이를 알고 있음) 또는 더 많은 데이터를 통합하지 않는한, 객체의 scale, 절대 위치와 canonical orientation은 항상 알려지지 않는다. 이는 일부 특성이 다른 것에 대해 compensate(보상하다)할 수 있기 때문이다. 예컨대 동일한 이미지를 얻으려면 객체를 뒤로 이동시키고 그에 따라 scale만 하면 된다. similarity ambiguity을 제거하는 한 가지 예는 camera calibration 절차에서 발생했는데, 여기서 world reference system에 대한 calibration 점의 위치를 알고 있다고 가정했다. 이를 통해 3d 구조의 metric scale을 학습하기 위해 체커보드의 정사각형의 크기를 알 수 있었다.

4 Perspective structure from motion

motion 문제에서 단순화된 affine structure를 공부한 후에 projective 카메라 MiM_i에 대한 일반적인 경우를 고려한다. projective 카메라의 일반적인 경우, 각 카메라 행렬 MiM_i는 scale까지 정의되어 11개 자유도를 포함한다.
Mi=[a11a12a13b1a21a22a23b2a31a32a331](4.1)M_i = \begin{bmatrix} a_{11} & a_{12} & a_{13} & b_1 \\ a_{21} & a_{22} & a_{23} & b_2 \\ a_{31} & a_{32} & a_{33} & 1 \end{bmatrix} \tag{4.1}
게다가 affine 경우에서 해가 affine transformation까지만 결정될 수 있었던 것처럼 일반적인 경우에 structure와 motion에 대한 해는 projective transformation까지만 결정될 수 있다. structure matrix에 inverse transformation H1H^{-1}를 적용하는 한, motion matrix에 임의의 4×44 \times 4 projective transformation HH를 적용할 수 있다. 이미지 평면에서 결과 관측치는 여전히 동일할 것이다.
affine 경우와 유사하게 motion 문제에서 일반적인 structure를 mnmn개 관찰 xijx_{ij}에서 mm개 motion 행렬 MiM_inn개 3d 점 XjX_j 모두를 추정하여 설정할 수 있다. 카메라와 포인트는 4×44\times 4 projective transformation(15개 파라미터)까지 scale까지만 복구될 수 있기 때문에 2mn2mn개의 등식에서 11m+3n1511m + 3n - 15개의 미지수를 갖는다. 이러한 사실에서 미지수를 해결하기 위해 필요한 view와 관찰의 수를 결정할 수 있다.

4.1 The algebraic approach

이제 algebraic(대수) 접근을 다룬다. 이것은 2대의 카메라에서 structure from motion 문제를 해결하기 위해 Fundamental 행렬 FF의 개념을 활용한다. 그림 5와 같이 대수 접근의 주요 아이디어는 perspective transformation HH까지만 계산할 수 있는 2대 카메라 행렬 M1M_1M2M_2를 계산하는 것이다. 각 MiM_i가 perspective transformation HH까지만 계산할 수 있기 때문에 첫 번째 카메라 projection 행렬 M1H1M_1H^{-1}이 canonical인 HH를 항상 고려할 수 있다. 물론 동일한 변환이 두 번째 카메라에도 적용될 수 있다.
M1H1=[I0]M2H1=[Ab](4.2)\begin{aligned} M_1H^{-1} &= [I \quad 0] \\ M_2H^{-1} &= [A \quad b ] \end{aligned} \tag{4.2}
이 작업을 수행하기 위해 먼저 이전 코스에서 다루었던 eight point algorithm을 사용하여 Fundamental 행렬 FF를 계산해야 한다. 이제 FF를 사용하여 projective 카메라 행렬 M1M_1M2M_2를 추정한다. 이 추정을 하기 위해 PP를 이미지의 대응 관측치 pppp'에 대한 대응 3d 점이라고 정의한다. H1H^{-1}을 두 카메라 projection 행렬에 적용하기 때문에 HH도 structure에 적용해야 한다. 따라서 P~=HP\tilde{P} = HP가 된다. 그러므로 픽셀 좌표 pppp'를 다음처럼 변환된 구조에 연관시킬 수 있다.
p=M1P=M1H1HP=[I0]P~p=M2P=M2H1HP=[Ab]P~(4.3)\begin{aligned} p &= M_1P = M_1H^{-1}HP = [I|0]\tilde{P} \\ p' &= M_2P = M_2H^{-1}HP = [A|b]\tilde{P} \end{aligned} \tag{4.3}
두 이미지 대응점 pppp' 사이의 흥미로운 속성이 일부 창의적인 대체에 의해 나타난다.
p=[Ab]P~=A[I0]P~+b=Ap+b(4.4)\begin{aligned} p &= [A|b]\tilde{P} \\ &= A[I|0]\tilde{P} + b \\ &= Ap + b \end{aligned} \tag{4.4}
방정식 4.4를 사용하여 pp'bb 사이에 cross product를 작성할 수 있다.
p×b=(Ap+b)×b=Ap×b(4.5)p' \times b = (Ap + b) \times b = Ap \times b \tag{4.5}
cross product의 정의에 따라 p×bp'\times bpp'에 수직이다. 그러므로 다음과 같이 작성할 수 있다.
0=p(p×b)=p(Ap×b)=p(b×Ap)=p[b]×Ap(4.6)\begin{aligned} 0 &= p'^\top(p'\times b) \\ &= p'^\top(Ap \times b) \\ &= p'^\top \cdot(b \times Ap) \\ &= p'^\top[b] _\times Ap \end{aligned} \tag{4.6}
이 제약조건을 보면 Fundamental 행렬 pFp=0p'^\top Fp = 0의 일반적인 정의가 떠오를 것이다. F=[b]×AF = [b]_\times A로 설정하면 AAbb를 추출하는 것은 decomposition 문제로 간단히 분해된다.
bb를 결정하는 것에서 시작하자. 다시 cross product의 정의에 의해 FbFb를 다음과 같이 간단히 작성할 수 있다.
Fb=[b]×Ab=(b×A)b=0(4.7)Fb = [b]_\times Ab = (b \times A)b = 0 \tag{4.7}
FF가 singular이므로 bb는 SVD를 사용하여 b=1\|b\| = 1Fb=0Fb = 0의 least square 해로 계산될 수 있다.
bb을 알면 AA를 계산할 수 있다. A=[b]×FA = -[b]_\times F로 설정하면 이 정의가 F=[b]×AF = [b]_\times A를 만족하는 것을 확인할 수 있다.
[b×]A=[b×][b×]F=(bbb2I)F=bbF+b2F=0+1F=F(4.8)\begin{aligned} [b_\times]A' &= -[b_\times][b_\times]F \\ &= (bb^\top - |b|^2I)F \\ &= bb^\top F + |b|^2F \\ &= 0 + 1\cdot F \\ &= F\end{aligned} \tag{4.8}
결과적으로 카메라 행렬 M1H1M_1H^{-1}M2H1M_2H^{-1}에 대한 두가지 표현식을 결정할 수 있다.
M~1=[I0]M~2=[[b×]Fb](4.9)\begin{aligned} \tilde{M}_1 &= [I \quad 0] \\ \tilde{M}_2 &= [-[b_\times]F & b] \end{aligned} \tag{4.9}
이 섹션을 마무리하기 전에 bb에 대한 geometrical 해석을 제공한다. 우리는 bbFb=0Fb = 0을 만족한다는 것을 안다. 이전 코스 노트에서 유도한 epipolar 제약조건을 떠올려라. 이는 이미지에서 epipole이 Fundamental 행렬에 의해 변환될 때 0에 매핑되는 점이라는 것을 발견했다(즉 Fe2=0Fe_2 = 0Fe1=0F^\top e_1 = 0). 그러므로 bb는 epipole임을 볼 수 있다. 이를 통해 카메라 projection 행렬에 대한 새로운 방정식 집합(방정식 4.10)을 얻는다.
M~1=[I0]M~2=[[e×]Fe](4.10)\begin{aligned} \tilde{M}_1 &= [I \quad 0] \\ \tilde{M}_2 &= [-[e_\times]F & e] \end{aligned} \tag{4.10}

4.2 Determining motion from the Essential matrix

대수적 접근으로 얻은 재구성을 개선하는 한 가지 유용한 방법은 calibrated 카메라를 사용하는 것이다. normalized 좌표에 대한 Fundamental 행렬의 특별한 경우인 Essential 행렬을 사용하여 카메라 행렬의 더 정확한 초기 추정을 추출할 수 있다. Essential 행렬 EE을 사용함으로써 카메라를 calibrate 하고 따라서 intrinsic 카메라 행렬 KK를 안다고 가정 한다. 우리는 normalized 이미지 좌표에서 직접 Essential 행렬 EE를 계산하거나 Fundamental 행렬 FF과 intrinsic 행렬 KK의 관계에서 계산할 수 있다.
E=KFK(4.11)E = K^\top FK \tag{4.11}
Essential 행렬이 calibrate 된 카메라를 가졌다고 가정하기 때문에 extrinsic 파라미터를 인코딩하는 5개의 자유도(두 카메라 사이의 rotation RR과 translation tt)만 갖는다는 것을 기억해야 한다. 다행히 이것은 우리의 motion 행렬을 생성하기 위해 추출하기를 원하는 정확한 정보이다. 우선 Essential 행렬 EE가 다음과 같이 표현될 수 있음을 떠올려라.
E=[t]×R(4.12)E = [t]_\times R \tag{4.12}
이와 같이, 아마도 EE를 두 성분으로 분해하는 전략을 발견할 수 있다. 우선 cross product 행렬 [t]×[t]_\times가 skew-symmetric이라는 것에 유의하라. 우리는 decomposition에서 사용할 2개의 행렬을 정의한다.
W=[010100001],Z=[010100000](4.13)W = \begin{bmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}, Z = \begin{bmatrix} 0 & 1 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \tag{4.13}
나중에 사용할 중요한 속성 하나는 부호까지 Z=diag(1,1,0)WZ = \text{diag}(1, 1, 0)W라는 것이다. 유사하게 부호까지 ZW=ZW=diag(1,1,0)ZW=ZW^\top = \text{diag}(1,1,0)라는 사실을 사용할 것이다.
eigenvalue decomposition의 결과로 일반적인 skew-symmetric 행렬의 block decomposition을 scale까지 생성할 수 있다. 따라서 [t]×[t]_\times를 다음과 같이 작성할 수 있다.
[t]×=UZU(4.14)[t]_\times = UZU^\top \tag{4.14}
여기서 UU는 어떤 orthogonal 행렬이다. 그러므로 decomposition을 다음과 같이 재작성할 수 있다.
E=Udiag(1,1,0)(WUR)(4.15)E = U \text{diag}(1, 1,0)(WU^\top R) \tag{4.15}
이 표현식을 조심스럽게 살피면, Singular Value Decomposition E=UΣVE = U\Sigma V^\top과 밀접하게 유사함을 볼 수 있다. 여기서 Σ\Sigma는 2개의 동일한 singular 값을 포함한다. EE를 scale까지 알고 E=Udiag(1,1,0)VE = U \text{diag}(1,1,0)V^\top 형식을 취한다고 가정하면, 다음과 같은 EE의 분해에 도달한다.
[t]×=UZUR=UWV or UWV(4.16)\begin{aligned} [t]_\times &= UZU^\top \\ R &= UWV^\top \text{ or } UW^\top V^\top \end{aligned} \tag{4.16}
주어진 분해가 유효함을 검사를 통해 증명할 수 있고 또한 다른 분해가 없음을 증명할 수 있다. [t]×[t]_\times의 형식은 그 left null space가 EE의 null space와 동일해야 한다는 사실에 의해 결정된다. unitray 행렬 UUVV가 주어지면 임의의 rotation RRUXVUXV^\top으로 decompose 될 수 있다. 여기서 XX는 또 다른 rotation 행렬이다. 이 값들을 대입하면 scale까지 ZX=diag(1,1,0)ZX = \text{diag}(1, 1, 0)가 된다. 따라서 XXWWWW^\top과 일치해야 한다.
EE의 분해가 행렬 UWVUWV^\topUWVUW^\top V^\top가 직교한다는 것만 보장한다. RR이 유효한 rotation이라는 것을 보장하려면 단순히 RR의 행렬식이 양수임을 확인하면 된다.
R=(detUWV)UWV or (detUWV)UWV(4.17)R = (\det UWV^\top)UWV^\top \text{ or } (\det UW^\top V^\top)UW^\top V^\top \tag{4.17}
rotation RR이 잠재적인 2개의 값을 취할 수 있는 것과 유사하게 translation 벡터 tt도 여러 값을 취할 수 있다. cross product의 정의에서 다음을 알 수 있다.
t×t=[t]×t=UZUt=0(4.18)t \times t = [t]_\times t = UZU^\top t = 0 \tag{4.18}
UU가 unitary라는 것을 알면, [t]×F=2\|[t]_\times \|_F = \sqrt{2}임을 알 수 있다. 그러므로 이 분해에서 tt에 대한 추정치는 위의 방정식과 EE가 scale까지 알려져 있다는 사실로부터, 다음을 의미한다.
t=±U[001]=±u3(4.19)t = \pm U \begin{bmatrix}0 \\ 0 \\ 1 \end{bmatrix} = \pm u_3 \tag{4.19}
여기서 u3u_3UU의 3번째 행이다. 검사에 의해 [t]×=UZU[t]_\times = UZU^\top를 벡터 vv로 부호까지 재형식화하여 동일한 결과를 얻을 수 있다는 것을 확인할 수도 있다.
그림 6에 나온대로, RRtt 모두에 대해 각각 2가지 옵션이 존재하기 때문에 잠재적인 4개의 R,tR, t 쌍이 존재한다. 직관적으로 4개 쌍은 특정한 방향에서 카메라를 회전시키거나 반대 방향에서 카메라를 회전시키는 모든 가능한 쌍과 특정한 방향 또는 반대 방향으로 translating 하는 옵션을 결합한 것을 포함한다. 그러므로 이상적인 조건하에 올바른 R,tR, t 쌍을 결정하기 위해 단 한 점만 triangulate하면 된다. 올바른 R,tR, t 쌍에 대해 triangulated 점 P^\hat{P}는 두 카메라의 앞에 존재하며, 이것은 두 카메라 reference system 측면에서 양의 zz-좌표를 갖는다는 것을 의미한다. 측정 노이즈 때문에 우리는 종종 단일 점만으로 triangulating 하지 않고 대신 많은 점을 triangulate 하고 두 카메라 앞에 가장 많은 점을 포함하는 R,tR, t 쌍을 올바른 쌍으로 결정한다.

5 An example structure from motion pipeline

상대적인 motion 행렬 MiM_i을 발견한 후에 이를 사용하여 점 XjX_j의 월드 좌표를 결정할 수 있다. 대수적 방법의 경우, 그런 점의 추정치는 perspective 변환 HH까지 정확하다. Essential 행렬에서 카메라 행렬을 추출할 때 추정치는 scale까지만 알 수 있다. 두 경우 모두 앞서 설명한 triangulartion 방법을 통해 추정된 카메라 행렬에서 3d 점을 계산할 수 있다.
multi-view 경우에 대한 확장은 쌍별 카메라를 연결하여 수행될 수 있다. 대응점이 충분하면 대수적 접근 또는 Essential 행렬을 사용하여 임의의 카메라 행렬에 3d 점의 해를 얻을 수 있다. 재구성된 3d 점은 카메라 쌍 사이에서 사용 가능한 대응점과 연관된다. 이러한 쌍별 해는 다음에 보게 될 bundle adjustment이라는 접근법으로 결합(최적화) 할 수 있다.

5.1 Bundle adjustment

지금까지 논의한 structure from motion 문제를 해결하기 위한 이전 방법들에는 주요한 한계가 존재한다. factorization 방법은 모든 점이 모든 이미지에서 visible이라고 가정한다. 이것은 occlusion이나 대응점 찾기 실패 때문에 이미지가 많거나 멀리 떨어져 촬영된 경우 일어날 가능성이 매우 없다. 마지막으로 대수적 접근은 카메라 체인으로 결합할 수 있는 쌍별 해를 생성하지만 모든 카메라와 3d 점을 사용하여 coherent(일관된) 최적화 재구성을 해결하지 않는다.
이러한 한계를 해결하기 위해 bundle adjustment를 소개한다. 이것은 structure from motion 문제를 해결하기 위한 비선형 방법이다. 최적화에서 우리는 reprojection error를 최소화하는데 초점을 맞춘다. 이것은 재구성된 점을 추정된 카메라로 projection할 때와 모든 카메라와 모든 점에 대한 대응 관측치 사이의 픽셀 거리를 의미한다. 이전에 triangulation을 위한 비선형 최적화 방법을 논의할 때 주로 2대 카메라 경우에 초점을 맞추었는데, 여기서는 자연스럽게 각 카메라가 두 카메라 사이의 모든 대응점을 보는 것으로 가정했다. 그러나 bundle adjustment는 여러 카메라를 다루므로 각 카메라가 볼 수 있는 관찰에 대해서만 reprojection error를 계산한다. 그렇지만 궁극적으로 이 최적화 방법은 triangulation에 대한 비선형 방법을 논의할 때 소개된 것과 매우 유사하다.
bundle adjustment의 비선형 최적화를 해결하기 위한 두 가지 일반적인 접근은 Gauss-Newton algorithm과 Levenberg-Marquardt algorithm이다. Gauss-Newton algorithm에 대해서는 이전 섹션에서 상세한 내용을 찾을 수 있고 Levenberg-Marquardt algorithm에 대해서는 Hartley와 Zisserman 교재에서 상세한 내용을 찾을 수 있다.
결론적으로 bundle adjustment는 우리가 조사한 다른 방법과 비교하여 중요한 이점과 한계를 갖는다. 이것은 많은 수의 view를 smoothly 다룰 수 있고 특정한 점이 모든 이미지에서 관찰되지 않은 경우도 다룰 수 있다는 점에서 특히 유용하다. 그러나 주요한 한계는 파라미터가 view의 수에 따라 증가함에 하므로 특히 큰 최소화 문제라는 것이다. 또한 이것은 비선형 최적화 기법에 의존하기 때문에 좋은 초기화 조건이 필요하다. 이러한 이유로 bundle adjustment는 종종 structure from motion 구현에서 마지막 단계로 사용되는 경우가 많다(즉 factorization 또는 algebraic 접근 이후에). factorization 또는 algebraic 접근이 최적화 문제에 대한 좋은 초기 해를 제공할 수 있기 때문이다.