[개념 정리] Bundle Adjustment

Bundle Adjustment란?

이미지들을 입력으로 '각 이미지의 camera parameter'와 '공유되는 3d point위치'를 예측하는 알고리즘입니다. 주로 SfM과 SLAM분야에서 사용됩니다.

압축해서 말하자면, 주황색 화살표 길이가 짧아지도록 빨간색점들을 이동시키는 알고리즘입니다. 상세 내용에 대해 길게 풀어 적어보겠습니다.

 

 

 

기본 아이디어

camera paramter로 3d point를 image로 projection했을 때, 2d point들의 color차이를 최소화하는 방식으로 작동합니다. 위 그림에서 3d point $X_i$과 camera parameter $C_j$은 예측값이고, i번째 point가 j번째 camera로 projection(=π)된 2d point $u_{ij}$는 입력(=관측)값입니다. 

각 이미지의 width, height(=n,m)에 대해 위 수식을 최소화하는 $X_i$와 $C_j$를 찾습니다. 초기값이 전혀 주어지지 않는 상태에선 해를 찾는데 시간이 오래 소요되므로, 어느 정도 값이 초기화 된 형태가 일반적입니다. 초기값에 대해선 SfM 글을 참고바랍니다.

 

 

 

모델 매개변수

모델 매개변수는 예측하고자 하는 값이며, 3d point와 camera parameter가 있습니다. 3d point는 단순히 x,y,z로 구성되며, camera parameter의 경우 카메라 모델에 대한 이해가 필요합니다. camera parameter $C_n$은 2d에 관한 intrinsic parameter, 3d에 관한 extrinsic parameter로 나뉘어집니다.

intrinsic은 2d translation, 2d scaling(=f=focal length),  2d shear(=s)으로 구성되어 있으며, 일반적인 경우 2d translation은 이미지 width, height를 나누기 2하면($x_0$ = w/2 ; $y_0$=h/2) 되고, shear s와 focal length f를 예측합니다.

 

extrinsic은 3d translation, 3d rotation으로 구성되어 있으며 모든 parameter 예측이 필요로 합니다.

카메라 3d 위치를 나타내는 3d translation은 $t_n$ 이고, 카메라의 회전 정도를 나타내는 3d rotation은 $r_{ij}$으로써 3x3 matrix로 구성됩니다. 3d point(=$x$)에 대해 matrix $P$로 matrix연산을 해주면 2d pixel(=$y$)로 projection(=π)됩니다($y = Px$). 보다 상세한 설명은 intrinsic/extrinsic parameter글 참고 바랍니다.

 

 

 

Bundle Adjustment 알고리즘

앞에서 언급한 수식은 비선형 최소 제곱문제(Nonlinear Least Squares Problem)입니다. $$\min \sum_{i=1}^n \sum_{j=1}^m (u_{ij} - \pi(C_j, X_i))^2$$ 이 문제는 Gradient-Descent, Gauss-Newton, Levenberg-Marquardt 방법으로 해결 가능하며, 이 중에서 Gauss-Newton(이전글) 방식으로 전개해보려합니다. $$ J^T J \Delta \beta = J^T r$$ 위 수식에 대입해서 설명하자면, $ r_{ij} = u_{ij} - \pi(C_j, X_i) = u_{ij} - \hat{u_{ij}} $이고, $i$는 3d point index, $j$는 camera index, $J$는 r을 각 모델 매개변수로 편미분한 Jacobian matrix입니다. $ \Delta \beta $는 모델 매개변수에 대한 업데이트 변화량이며, 최종적으로 계산하고자하는 패러미터입니다. Jacobian은 predicted되는 $\hat{u_{ij}}$ 에 대해서 parameter $x=[c_{j},p_{i}]$에 관한 수식으로 편미분됩니다. $$ J_{ij} = \frac{\partial r_{ij}}{\partial x_k} = \frac{\partial \hat{u}_{ij}}{\partial x_k} $$ $c$는 camera parameter, $p$는 3d point입니다. Jacobian의 경우, 편미분되는 모델 매개변수 갯수(=g) 만큼 열이 되고, 함수의 갯수(=d)만큼이 행이 됩니다. 여기서 함수의 갯수는 summation되는 횟수($d=i*j$)입니다. $$ J = \left[ \frac{\partial f}{\partial x_1} \cdots \frac{\partial f}{\partial x_g} \right] = \begin{bmatrix} \nabla f_1 \\ \vdots \\ \nabla f_d \end{bmatrix} = \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_g} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_d}{\partial x_1} & \cdots & \frac{\partial f_d}{\partial x_g} \end{bmatrix} $$ Jacobian 예시를 들어보겠습니다. 카메라가 3개, 3D point가 4개인 경우의 예시입니다. 아래는 Jacobian 수식입니다. 각 카메라에 대한 모델 매개변수들(focal, shear, translation 3x1, rotation 3x3)을 모두 표시하면 많으므로, 카메라마다 한개 변수(=$A_{i,j}$)로만 표시했습니다.

위와 같은 형태를 sparse matrix라고 합니다. $A_{i,j}$는 실제로 카메라 매개변수 갯수만큼 열이 길어지는 형태가 됩니다. 모델 매개변수 섹션 마지막줄에 언급된 $y=Px$ 다항식에 대해, 편미분하면 위 matrix를 채울 수 있습니다.

 

다음으로 Gauss-Newton의 좌측식을 계산합니다.

$$ U_j = \sum_{i=1}^{4} A_{ij}^T A_{ij}, \quad V_i = \sum_{j=1}^{3} B_{ij}^T B_{ij}, \quad W_{ij} = A_{ij}^T B_{ij} $$

아래는 Gauss-Newton의 우측식에 해당합니다. 실제론 matrix형태이지만, 연산과정을 설명하기 위한 형태로 표기되었습니다. 위에 언급된 $J$를 Transpose해서 $r$을 곱해보면 의미가 와닿으실 겁니다.

다시 Gauss-Newton식을 정리하면, 아래 수식이 됩니다. $$ J^T J \Delta \beta = J^T r$$ $$\begin{bmatrix}
U & W \\
W^T & V
\end{bmatrix}
\begin{bmatrix}
\Delta c \\
\Delta p
\end{bmatrix}
=
\begin{bmatrix}
r_c \\
r_p
\end{bmatrix}$$

수식이 복잡해보이지만, 실제론 Matrix를 이용한 곱셈 덧셈연산을 사용해서, 처음에 언급한 최소제곱식을 최소화하도록 모델 매개변수들의 변화량  $\Delta c$와 $\Delta p$를 계산 할 수있는 방정식이 완성됬습니다. 

 

 

 

Reduced Camera System 기법

위 수식은 모델 매개변수가 적을 경우(ex, 카메라수와  3d point가 적은 경우)엔 충분히 사용 가능합니다. 모델 매개변수가 많아질 경우엔, 연산량을 줄여주도록 선형대수학의 Schur Complement(슈어 보수) 개념을 같이 사용하게 됩니다.

카메라의 갯수가 3d point 개수보다 훨씬 적은 경우, 3d point 변화량 $Delta p$를 계산하지 않고도 카메라 패러미터 변화량 $\Delta c$ 를 먼저 계산한 후에, $\Delta p$를 계산하는 전략을 취할 수 있습니다. 조금전에 언급한 수식에, $[I - WV^{-1}]^T$를 양변 왼쪽에 곱해서 $\Delta p$를 제거해줍니다. $$ \left[ \begin{array}{cc}
U - WV^{-1}W^T & 0 
\end{array} \right]
\left[ \begin{array}{c}
\Delta c \\
\Delta p
\end{array} \right] = 
\left[ \begin{array}{c}
I - WV^{-1}
\end{array} \right]^T
\left[ \begin{array}{c}
r_c \\
r_p
\end{array} \right]$$ $$(U - WV^{-1}W^T) \Delta c = r_c - WV^{-1} r_p$$ 마지막 수식을 Reduced Camera System이라고 칭한다고 하며, $ S = U - WV^{-1}W^T $부분이 V에 대한 슈어 보수라고 합니다. 이를 통해, 계산해야 할 $\Delta p$가 사라지면서, 계산 복잡도가 줄어들게 됩니다. S는 단순 matrix 연산입니다. 앞에서 언급드렸던 matrix 연산의 연장선이므로, 여기서부터는 복잡한 세부수식을 생략하겠습니다.

 

위 수식으로 현 시점의 미지수인 $\Delta c$ 를 계산 할 수 있습니다. 이제  $\Delta p$를 계산 할 차례입니다. 아래는 앞에 언급한 수식입니다.

빨간색 부분만 갖고와서, $\Delta p$에 관한 수식으로 만들 수 있고, $\Delta p$를 계산 할 수 있습니다.

$$ W^T \Delta c + V \Delta p = r_c $$

$$ \Delta p = V^{-1} (r_p - W^T \Delta c) $$ 여기까지가 $\Delta p$와 $\Delta c$를 분리해서 계산해서 연산복잡도를 낮추는 방법에 관한 내용입니다. 시간복잡도는 $O((c+p)^3)$에서 $O(c^3+cp)$로 줄어들게 됩니다.

비선형 최소 제곱 최적화 알고리즘 중 Gauss Newton을 사용한 방식을 설명드렸고, Levenber-Marquardt으로 할 경우에는 U,V의 대각성분을 augment하면 된다고 합니다. 그 외 추가로 자료를 찾다보니 $\Delta p$부터 구하고 $\Delta c$를 구하는 경우도 있습니다. 3d point 개수가 카메라의 갯수보다 적은 경우이지 않을까합니다.  

 

 

 

Conjugate Gradient기법

위에서 언급한 기법은 관측된 데이터(ex. feature matching에 에러)에 nosie가 많을 경우, 수렴 속도가 느려지게 됩니다. 이를 위해 또 다른 최적화 알고리즘인 Conjugate Gradient를 사용합니다. $Ax = b$ 에서 A matrix가 매우 크고 0이 많고 에러가 수반되는 경우 iterative방법으로 x를 계산하는 방법론입니다.

$$ \min q(x) = \frac{1}{2} x^T A x - b^T x $$ 1차 방정식을 2차 방정식으로 만들어 함수 $q$를 최소화하는 x를 계산합니다.

요약하자면, $b-Ax$를 오차값(s)으로 해서, loop를 통해 s가 지정한 threshold값까지 작아질 때까지 x를 갱신하면서 반복하는 알고리즘입니다. threshold를 넘기면 되므로, 정확한 값이 아닌 근사치를 계산 할 수 있습니다. $$ J^T J \Delta \beta = J^T r$$ 앞에서 언급된 위 수식을 아래 수식으로 만듭니다. 여기서 $\eta_j$는 0.1입니다.
$$ \| J^T J \Delta \beta + J^T \| \leq \eta_j \| J^T \| $$

Conjugate Graident기법에 대해서는 기본 아이디어만 전달하며, 글을 마치겠습니다.

 

 

 

출처

Bundle Adjustment Revisited, 2019