이 프로젝트는 Camera Calibration과 Disparity Map, Depth Map을 생성하는 프로젝트이다.

Camera Calibration

먼저, 체커보드를 이용해서 Camera Calibration을 수행한다. 6개의 체커보드 이미지가 있고, 모서리를 찾아서 이를 평면좌표(world point)로 만들어주는 homography matrix를 찾고, 그 matrix에서 intrinsic, extrinsic 파라미터를 찾는 과정이다. Zhang’s method를 사용하여 수행할 수 있다.

Homography Calculation

먼저 $s\tilde{q}=H\tilde{p}$를 만족하는 행렬 $H$를 찾아야 한다. ($\tilde{q}$는 world point, $\tilde{p}$는 체커보드 좌표)

$H = \left[\begin{matrix}h_{11} & h_{12} & h_{13}\\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33}\end{matrix}\right]$

으로 표현되고, $\min_H\sum_{j=1}^m\left|\left|q_j-\hat{q_j}\right|\right|^2$를 최소화하는 $H$를 찾으면 된다. ( $q_j = \left[\begin{matrix}u_j \\ v_j\end{matrix}\right]$ , $\hat{q_j} = \frac{1}{h_{31}X_j+h_{32}Y_j+h_{33}}\left[\begin{matrix}h_{11}X_j+h_{12}Y_j+h_{33}\\h_{21}X_j+h_{22}Y_j+h_{23}\end{matrix}\right]$ )

이를 정리하면,

$\min_H\sum_j\left|\left|\left[\begin{matrix}u_0(h_{31}X_j+h_{32}Y_j+h_{33})-(h_{11}X_j+h_{12}Y_j+h_{13})\\v_0(h_{31}X_j+h_{32}Y_j+h_{33})-(h_{21}X_j+h_{22}Y_j+h_{23})\end{matrix}\right]\right|\right|^2$

이 되는데, 각 항을 풀어보면

$\left[\begin{matrix}-X_j & -Y_j & -1 & 0 & 0 & 0 & u_jX_j & u_jY_j & u_j \\ 0 & 0 & 0 & -X_j & -Y_j & -1 & v_jX_j & v_jY_j & v_j\end{matrix}\right] \left[\begin{matrix}h_{11}\\ h_{12} \\h_{13} \\h_{21} \\h_{22} \\h_{23} \\h_{31} \\h_{32} \\h_{33}\end{matrix}\right]$

이라고 쓸 수 있다. 앞에 곱해진 행렬을 $L_j$라고 하면 이 문제는 $\min_{\mathbf{x}}||L\mathbf{x}||^2$ 문제가 된다. ($L$은 $L_j$를 세로로 늘어뜨린 행렬) 따라서 $L^TL$의 가장 작은 eigenvalue중 가장 작은 것에 대응하는 eigenvector를 찾으면 된다.

이 때의 $H$ 행렬은 한 이미지당 하나를 찾을 수 있으므로 이 프로젝트에서는 6개의 $H$ 행렬이 나오게 된다.

Intrinsic Parameter

우리는 $H$행렬을 알고 있고, $H = K\left[\begin{matrix}r_1 & r_2 & t\end{matrix}\right]$ 이다. 그리고 $r_1^Tr_1=r_2^Tr_2$ , $r_1^Tr_2=0$ 이라는 제약을 갖고 있으므로, 우리는 intrinsic matrix $K$를 알아낼 수 있다.

먼저 $H=\left[\begin{matrix}h_1 & h_2 & h_3\end{matrix}\right]$ 이라고 하자. $h_1=Kr_1$, $h_2=Kr_2$이기 때문에 $r_1=K^{-1}h_1$, $r_2=K^{-1}h_2$이다.

위 조건들을 대입해 보면

$h_1^TK^{-T}K^{-1}h_1=h_2^TK^{-T}k^{-1}h_2$ $h_1^TK^{-T}K^{-1}h_2=0$

라는 식을 만들어낼 수 있다.

$B=K^{-T}K^{-1}=\left[\begin{matrix}B_{11} & B_{12} & B_{13} \\B_{21} & B_{22} & B_{23} \\B_{31} & B_{32} & B_{33}\end{matrix}\right]$

이라고 하면,

$B=\left[\begin{matrix}\frac{1}{\alpha^2} & -\frac{\gamma}{\alpha^2\beta} & \frac{v_0\gamma-u_0\beta}{\alpha^2\beta} \\ -\frac{\gamma}{\alpha^2}{\beta} & -\frac{\gamma}{\alpha^2\beta}+\frac{1}{\beta^2} & -\frac{\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}\\ \frac{v_0\gamma-u_0\beta}{\alpha^2\beta} & -\frac{\gamma(v_0\gamma-u_0\beta)}{\alpha^2\beta^2}-\frac{v_0}{\beta^2} & \frac{(v_0\gamma-u_0\beta)^2}{\alpha^2\beta^2}+\frac{v_0^2}{\beta^2}+1\end{matrix}\right]$

이다. ( $K=\left[\begin{matrix}\alpha & \gamma & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1\end{matrix}\right]$ 이므로)

$B_{11}, B_{12}, B_{13}, B_{22}, B_{23}, B_{33}$만 알면 $K$의 변수들을 구할 수 있으므로 $b=\left[\begin{matrix}B_{11} & B_{12} & B_{13} & B_{22} & B_{23} & B_{33}\end{matrix}\right]^T$ 를 구하는 것을 목표로 한다.

두 개의 조건을 압축해서

$\left[\begin{matrix}v_{12}^T \\ (v_{11}-v_{22})^T\end{matrix}\right]b=0$

이라고 쓸 수 있다. ( $v_{kl}=\left[\begin{matrix}h_{1k}h_{1l} & h_{1k}h_{2l}+h_{2k}h_{1l} & h_{1k}h_{3l}+h_{3k}h_{1l} & h_{2k}h_{2l} & h_{2k}h_{3l}+h_{3k}h_{2l} & h_{3k}h_{3l} \end{matrix}\right]^T$ )

$V=\left[\begin{matrix}v_{12}^T \\ (v_{11}-v_{22})^T\end{matrix}\right]$ 라고 하면 각 $V$들은 이미지당 하나씩 있으므로, 우리는 총 6개의 $V$를 가지고 $B$를 찾을 수 있다. 이 $V$들을 세로로 합쳐서 $V’$행렬을 만들고, $\min_b\|V'b\|^2$ 문제를 풀면 $b$를 구할 수 있다. 이 때의 해는 $V’^TV’$의 eigenvalue중 가장 작은 것에 대응하는 eigenvector를 구하면 된다.

이렇게 구한 b를 이용해서 각각의 변수에

$v_0=\left(\frac{B_{12}B_{13}-B_{11}B_{23}}{B_{11}B_{22}-B_{12}^2}\right)$ $\lambda=B_{33}-\frac{B_{13}^2+v_0(B_{12}B_{13}-B_{11}B_{23})}{B_{11}}$ $\alpha=\sqrt{\frac{\lambda}{B_{11}}}$ $\beta=\sqrt{\frac{\lambda B_{11}}{B_{11}B_{22}-B_{12}^2}}$ $\gamma=-\frac{B_{12}\alpha^2\beta}{\lambda}$ $u_0=\frac{\gamma v_0}{\beta}-\frac{B_{13}\alpha^2}{\lambda}$

를 대입하면 intrinsic matrix $K=\left[\begin{matrix}\alpha & \gamma & u_0 \\ 0 & \beta & v_0 \\ 0 & 0 & 1\end{matrix}\right]$ 를 구할 수 있다.

Extrinsic Parameter

지금까지 구한 $H$와 $K$로 extrinsic parameter들을 구할 수 있다. 먼저 정규화를 위한 상수 $\lambda'=\frac{\frac{1}{\|K^{-1}h_1\|}+\frac{1}{\|K^{-1}h_2\|}}{2}$ 를 구한다. 그 후에는

$r_1=\lambda' K^{-1}h_1$ $r_2=\lambda' K^{-1}h_2$ $r_3=r_1\times r_2$ $t=\lambda'K^{-1}h_3$

을 이용해 extrinsic parameter들을 구한다.

여기서 구한 rotation matrix $R=\left[\begin{matrix}r_1 & r_2 & r_3\end{matrix}\right]$ 은 일반적으로 rotation matrix의 성질을 만족하지 않으므로, 보정을 해줘야 한다. $R=U\Sigma V^T$로 분해할 수 있는데, 여기서 $R’=UV^T$를 구하면 된다.

Maximum Likelihood Estimation (MLE)

우리가 위에서 구한 값은 그저 방정식을 푼 것이지, 물리적으로 아무 의미가 없는 숫자이다. 따라서 MLE를 통해 값을 보정해줄 필요가 있다. MLE추정은 $\min_{K, R_i, t_i}\sum_i\sum_j\|q_{ij}-\hat{q_{ij}}\|^2$ 를 통해 할 수 있다. (사실 이 식은 $H$를 찾을 때 썼던 식이다.)

코드에서는 intrinsic, extrinsic 파라미터들을 벡터로 만들어 넘겨서 lsqnonlin 함수를 사용하여 최적화했다.(Matlab)

Depth Map

Depth map을 구하기 위한 과정은 다음과 같다. (좌, 우에서 찍은 이미지가 주어진다.)

이미지를 흑백화시킨다.
cost function을 정해서 cost volume을 만든다.
cost aggregation을 한다.
최소 cost를 갖는 인덱스를 골라 disparity map을 만든다.
disparity map을 카메라 파라미터를 이용해 depth map으로 변환시킨다.

Preprocessing

단순히 두 이미지를 흑백화시킨다.

Cost Volume

좌측 이미지를 기준으로 한다. 좌측 이미지의 모든 픽셀에 대하여 오른쪽 이미지에서는 왼쪽으로 $d$만큼 떨어진 픽셀을 고른다. (좌측 이미지에서 $(10,5)$에 대해, $d=5$라면 우측 이미지에서는 $(5,5)$를 고르는 셈.) 좌우 이미지 모두 해당 픽셀 주변의 window크기만큼의 행렬을 골라 벡터화시킨 후, cost를 계산해 cost volume에 저장한다. 이 계산을 minDisparity부터 maxDisparity까지 수행해 cost volume에 저장한다. (이 프로젝트에서 window 크기는 15로 지정하였고, minDisparity는 11, maxDisparity는 140으로 지정하였다. cost 대신에 $NCC=\frac{\sum_i\sum_jA_{ij}B_{ij}}{\sqrt{\sum_i\sum_jA_{ij}^2}\sqrt{\sum_i\sum_j{B_{ij}^2}}}$ 를 사용하였다.)