[Camera model] Camera calibration

3차원 물체를 카메라로 관측하고, 이미지로 투영되는 과정을 생각해보자. 

동일한 3차원 물체를 카메라로 촬영해서 나오는 결과 이미지가 모두 같을까?

절대 그렇지 않다.

 

 

동일한 3차원 물체를 관측해도 대포 카메라로 촬영하는 것과 고프로를 사용해서 촬영하면 결과 이미지가 다르다.

동일한 카메라를 사용한다 하더라도 사진을 어디서 찍느냐에 따라서 결과 이미지가 달라질 수 있다. 

 

 

이것을 수학적으로 모델링 한 것이 바로 camera parameter이다. 

camera parameter는 또 다시 2가지로 구분 될 수 있는데, 

카메라 자체에서 결정되는 intrinsic parameter와, 

카메라의 자세를 통해 결정되는 extrinsic parameter로 구분할 수 있다. 

 

 

위와 같은 camera parameter를 구하는 작업을 camera calibration이라고 부른다.

 

이번 글에서는 camera calibration의 대표적인 방법 2가지에 대해서 알아본다.

굉장히 풀어서 단계별로 짚어가려고 하니 나와 같은 초심자들이 이해하기 좋을 것 같다!

 

• Direct Linear Transform (DLT)
• Zhang's Method

 

 

---

 

 

Direct Linear Transform

Direct Linear Transform은 world coordinate의 3차원 점 $X$ 와 대응되는 image coordinate의 2차원 점 $x$의 대응쌍들을 통해서 projection matrix $P$를 구하는 알고리즘이다.

 

3D point : $X = (X_i, Y_i, Z_i)$

2D point : $x = (x_i, y_i)$

 

3D point → 2D point : $x = PX$

 

위에 나타난 projection 식은 아래와 같이 쓸 수 있다. 

$$x = PX$$

$$x = K [ R\ |\  t]X$$

$$x = K R [I\ |\ -C]X$$

 

여기서 $K$는 camera intrinsic parameter으로 $f_x, f_y, c_x, c_y, skew_c$ 총 5개의 parameter를 갖고, 

$R, t$는 각각 3개의 parameter를 갖는다. 

 

 

따라서 $P$는 총 11개(5 + 3 + 3)의 unknown 값을 갖고, 

이 11개의 unknown 값을 구하기 위해서는 11개의 식이 필요하게 된다. 

DLT에서는 이 11개의 식을 3D point와 2D point간의 correspondence로부터 뽑아내어 사용한다. 

 

(※ 만약 calibrated camera라면 intrinsic을 이미 알고 있기 때문에, extrinsic에 관련된 6개의 unknown 값을 가지게 되고, 이 unknown 값을 알아내기 위해서 6개의 식을 필요로 하게 된다.)

 

그러면 이 식들을 어떻게 뽑아내는지 알아보도록 하자.

 

 

먼저 한 쌍의 3D point - 2D point 대응쌍을 나타내는 식을 살펴보자. 

 

$$\begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{14} \\ p_{21} & p_{22} & p_{23} & p_{24} \\ p_{31} & p_{32} & p_{33} & p_{34} \end{bmatrix} \begin{bmatrix} X \\ Y \\ Z \\ 1 \end{bmatrix}$$

 

여기서 2D point의 좌표는 아래와 같이 계산되어진다.

$$x = \frac{p_{11}X + p_{12}Y + p_{13}Z + p_{14}}{p_{31}X + p_{32}Y + p_{33}Z + p_{34}}$$

 

$$y = \frac{p_{21}X + p_{22}Y + p_{23}Z + p_{24}}{p_{31}X + p_{32}Y + p_{33}Z + p_{34}}$$

 

위와 같이 우리는 하나의 대응쌍을 알면 $x, y$ 각각에 대한 2개의 식을 뽑아낼 수 있다는 사실을 알 수 있다. 

여기서 우리가 사용할 DLT 방법에서는 11개의 식이 필요했으므로,

최소 6개의 대응쌍만 알면, 12개의 식을 얻을 수 있으므로, 우리는 DLT를 통해 $P$ matrix를 구할 수 있다.

 

 

 

먼저, 하나의 대응쌍을 ${x_i, X_i}$라고 정의하자.

그러면 $x_i = PX_i$ 라는 식이 성립하게 되고, $P$는 3x4 matrix 이므로 아래와 같이 쓸 수 있다.

 

$$x_i = \begin{bmatrix} p_{11} & p_{12} & p_{13} & p_{14} \\ p_{21} & p_{22} & p_{23} & p_{24} \\ p_{31} & p_{32} & p_{33} & p_{34} \end{bmatrix} X_i$$

 

여기서 $P$ matrix의 각 row를 원소로 갖는 column vector $P_1,\ P_2,\ P_3$를 정의해보자. 

그러면 위의 식은 아래와 같이 표현될 수 있다.

 

$$x_i = \begin{bmatrix} \color{Red}p_{11} & \color{Red}p_{12} & \color{Red}p_{13} & \color{Red}p_{14} \\ \color{Green}p_{21} & \color{Green}p_{22} & \color{Green}p_{23} & \color{Green}p_{24} \\ \color{Blue}p_{31} & \color{Blue}p_{32} & \color{Blue}p_{33} & \color{Blue}p_{34} \end{bmatrix} X_i$$

 

$$x_i = \begin{bmatrix} \color{Red} P_1^T \\ \color{Green} P_2^T \\ \color{Blue} P_3^T \end{bmatrix} X_i$$

 

$$\color{Red} P_1 = \begin{bmatrix} p_{11} \\ p_{12} \\ p_{13} \\ p_{14} \end{bmatrix}$$

$$\color{Green} P_2 = \begin{bmatrix} p_{21} \\ p_{22} \\ p_{23} \\ p_{24} \end{bmatrix}$$

$$\color{Blue} P_3 = \begin{bmatrix} p_{31} \\ p_{32} \\ p_{33} \\ p_{34} \end{bmatrix}$$

 

 

이렇게 정의한 $P_1,\ P_2,\ P_3$를 사용해보자. 

앞에서 정의했었던 3D point - 2D point 대응쌍이 주어졌을 때 얻을 수 있는 2개의 식을 복기 해보자. 

$x_i, X_i$는 3D point - 2D point 대응쌍 중 $i$번째 대응쌍을 의미한다.

 

$$x= \frac{p_{11}X + p_{12}Y + p_{13}Z + p_{14}}{p_{31}X + p_{32}Y + p_{33}Z + p_{34}}$$

 

$$y = \frac{p_{21}X + p_{22}Y + p_{23}Z + p_{24}}{p_{31}X + p_{32}Y + p_{33}Z + p_{34}}$$

 

이 식은 $P_1,\ P_2,\ P_3$를 사용해서 아래와 같이 간단하게 표기할 수 있다. 

 

$$x_i = \frac{P_1^TX_i}{P_3^TX_i}$$

 

$$y_i = \frac{P_2^TX_i}{P_3^TX_i}$$

 

 

여기서 우리가 구하고자 하는게 무엇인지 한번 리마인드 하자.

우리는 카메라 파라미터가 포함된 $P$ matrix를 구하고 싶다.

 

따라서 우리는 위의 $x_i, y_i$에 대한 식을 아래와 같이 $P_1, P_2, P_3$에 대해서 정리한다.

 

위의 식에서 양변에 공통으로 들어간 분모인 $P_3^T X_i$를 곱해주고 우변의 값을 모두 좌변으로 넘겨서 정리한다.

그러면 아래와 같은 식을 얻을 수 있다.

 

$$-P_1^TX_i +\quad \quad \quad \ \ + x_iP_3^TX_i = 0$$

 

$$\quad \quad \quad \ \ + -P_2^TX_i + y_iP_3^TX_i = 0$$

 

 

여기서 우리가 구하고자 하는 $P_1, P_2, P_3$에 붙어있는 transpose를 양변에 모두 취해서 아래와 같이 없애준다. 

(Note that $(AB)^T = B^T A^T$)

 

$$-X_i^TP_1 + \quad \quad \quad \ \ + x_i X_i^TP_3 = 0$$

 

$$\quad \quad \quad \ \ + -X_i^TP_2 + y_iX_i^TP_3 = 0$$

 

 

그리고 위의 두 식을 Matrix form으로 아래와 같이 묶어줄 수 있다. 

 

$$\begin{bmatrix} -X_i^T & 0 & x_iX_i^T \\ 0 & -X_i^T & y_iX_i^T \end{bmatrix}\begin{bmatrix} P_1 \\ P_2 \\ P_3 \end{bmatrix} = 0$$

 

 

현재는 하나의 대응쌍에 대해서 위와 같은 (2x12) (12x1) = 0 꼴의 Matrix form 방정식이 만들어졌다. 

6개의 대응쌍에 대해서 모두 위와 같은 작업을 거쳐주게 되면 Matrix form이 아래와 같이 (12x12) (12x1) = 0 꼴이 된다.

 

$$\begin{bmatrix} -X_1^T & 0 & x_1X_1^T \\ 0 & -X_1^T & y_1X_1^T \\ -X_2^T & 0 & x_2X_2^T \\ 0 & -X_2^T & y_2X_2^T \\ ... & ... & ... \\ -X_6^T & 0 & x_6X_6^T \\ 0 & -X_6^T & y_6X_6^T\end{bmatrix}\begin{bmatrix} P_1 \\ P_2 \\ P_3 \end{bmatrix} = 0$$

 

 

위의 식은 linear algebra에서 단골 손님으로 나오는 $Ax = 0$의 형태이다. 

위에서 Coefficient matrix를 $M$이라고 하고, 우리가 구하고자 했던 $P_1, P_2, P_3$으로 이루어진 행렬을 $P$라고 하자.

 

그러면 우리는 $M$에 Singular Value Decomposition을 적용해서 $P$를 구할 수 있다. 

 

 

$P$ Matrix를 구했지만 우리는 아직 끝나지 않았다. 

왜냐하면 $P$ matrix는 intrinsic과 extrinsic이 모두 곱해진 형태이므로, 각각에 대한 정보를 알 수 없기 때문이다. 

 

따라서 이제는 $P$의 decomposition을 통해 intrinsic 과 extrinsic을 각각 알아낼 것이다. 

먼저 $P$의 형태를 살펴보자. 

$P$는 intrinsic과 extrinsic의 곱이므로 아래와 같이 표현할 수 있다. 

 

$$P = K[R|t]$$

 

$$P = KR[I | -C]$$

 

$$P = [KR | -KRC]$$

 

이 때, 편의상 아래와 같이 $H, h$ matrix를 정의한다.

 

$$H = KR$$

$$h = -KRC$$

 

$H$는 3x3 matrix이고, $h$는 3x1 matrix이다.

$P$를 사전에 구했기 때문에 이를 3x3, 3x1로 나누면 $H, h$ 의 값을 알 수 있다. 

 

이제 우리는 $H, h$ 값을 통해 우리는 $K, R, C$의 값을 각각 구할 것이다. 

 

먼저 camera center인 $C$는 쉽게 구할 수 있다.

왜냐하면 아래와 같이 $C$를 $H, h$로 쉽게 표현할 수 있기 때문이다. 

 

$$h = -KRC = -HC$$

$$C = -H^{-1}h$$

 

 

이제는 intrinsic matrix $K$와 rotation matrix $R$을 구할 것이다.

우리는 $H$(정확히는 $H^{-1}$)에 QR decomposition을 적용하여 각각의 값을 구할 수 있다. 

 

QR decomposition을 통해 얻어지는 $Q$ matrix는 rotation matrix와 대응되고, $R$ matrix는 intrinsic과 대응된다. 

하지만 $H = KR$이기 때문에 이 상태로 그냥 QR decomposition을 적용하게 되면 $H = QR$과 비교했을 때

$Q$의 위치에 intrinsic matrix가, $R$의 위치에는 rotation matrix가 대응되게 된다. 

 

이 문제를 해결해주기 위해서, 우리는 아래와 같이 H의 inverse matrix인 $H^{-1}$에 QR decomposition을 적용한다.

 

$H = KR$

 

$H^{-1} = (KR)^{-1} = R^{-1}K^{-1} = \color{Red} R^T \color{Green} K^{-1}$ --- (1)

 

$H^{-1} = \color{Red} Q \color{Green} R$ --- (2)

 

따라서 QR decomposition의 결과로 나온 (2)번식의 $Q, R$ 값으로 

(1)번식의 $K, R$을 구해줄 수 있다. 

 

이렇게 해서 $K, R, C$를 모두 구할 수 있다. 

 

 

여기서 하나 고려해야 할 부분이 있다. 

$H$는 homogeneous matrix였기 때문에 $K$도 homogeneous matrix이다. 

 

원래 $K$의 (3,3) 원소 값은 1인데 (intrinsic matrix form),

현재 $K$는 homogeneous matrix이기 때문에 1이 아닐수도 있다. 

따라서 $K$의 (3,3) 원소를 1로 만들어 주기 위해서 현재 구한 $K$ matrix를 normalize 하는 과정이 필요하다. 

 

$$K = \frac{1}{K_{33}}K$$ 

 

 

 

---

 

 

 

Zhang's Method

 

Zhang's method는 camera의 extrinsic을 제외한 intrinsic parameter만을 간편하게 구하는 방법이다. 

Intrinsic parameter를 구하는 것이기 때문에 총 5개의 unknown을 구하면 된다. 

따라서 필요한 식의 수도 상대적으로 적다.

 

Zhang's method는 calibration을 위햇 아래와 같은 checkerboard pattern을 사용한다. 

(아마도 OpenCV로 calibration을 했다면 이런 checkerboard pattern을 사용해봤던 기억이 있을 것이다.)

 

 

Zhang's method는 앞에서 살펴보았던 DLT 방법과 논리는 동일한데, 약간 잔머리를 사용했다고 보면 된다. 

 

DLT에서 우리는 3차원 3D 점에 대한 정보를 알고, 해당 점과 대응되는 2D 점의 대응쌍들을 통해 calibration을 하였다. 

Zhang's method에서는 3차원 3D 점들을 checkerboard에서 검출한다.

 

이게 그래서 뭐가 좋은데?

 

Checkerboard 위에서 검출된 점들은 '한 평면에 놓인 점들'이라는 constraint를 갖는다. 

여기서 우리는 이 checkerboard가 world coordinate의 XY-plane 위에 놓여있다고 가정한다. 

(World coordinate는 우리가 기준을 정해주기 나름이다)

 

Checkerboard가 XY-plane에 놓여있으므로 checkerboard 위의 모든 3D 점들의 Z좌표는 0이 된다. 

이런 단순한 가정 하나가 어떤 계산적인 측면에서 이득을 가져오는지 투영식을 통해 알아보자.

 

DLT에서 사용했었던 일반적인 투영식은 아래와 같다. 

 

$$\begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & skew_c & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} R_{11} & R_{12} & R_{13} & t_{1} \\ R_{21} & R_{22} & R_{23} & t_{2} \\ R_{31} & R_{32} & R_{33} & t_{3} \end{bmatrix} \begin{bmatrix} X \\ Y \\ Z\\ 1 \end{bmatrix}$$

 

 

다만 우리는 checkerboard 위의 모든 점들이 Z = 0의 값을 가진다고 가정했기 때문에, 

3D 점의 Z 값은 항상 0의 값을 가지게 되고, Z와 직접적으로 곱해지는 $p_13, p_23, p_33$도 그 값이 무엇이던 간에 Z와 곱해진 이후 결국 0이 된다. 

따라서 위 식을 아래와 같이 간소화 할 수 있다.

 

$$\begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & skew_c & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} R_{11} & R_{12} & 0 & t_{1} \\ R_{21} & R_{22} & 0 & t_{2} \\ R_{31} & R_{32} & 0 & t_{3} \end{bmatrix} \begin{bmatrix} X \\ Y \\ 0 \\ 1 \end{bmatrix}$$

 

$$\begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = \begin{bmatrix} f_x & skew_c & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} R_{11} & R_{12} & t_{1} \\ R_{21} & R_{22} & t_{2} \\ R_{31} & R_{32} & t_{3} \end{bmatrix} \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}$$

 

 

여기서 편의상 $H$ matrix를 아래와 같이 정의하겠다. 

 

$$H = \begin{bmatrix} f_x & skew_c & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}\begin{bmatrix} R_{11} & R_{12} & t_{1} \\ R_{21} & R_{22} & t_{2} \\ R_{31} & R_{32} & t_{3} \end{bmatrix}$$

 

 

그러면 위의 식은 아래와 같이 $H$를 통해 표현할 수 있다. 

 

$$\begin{bmatrix} x \\ y \\ 1 \end{bmatrix} = H \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}$$

 

 

이제는 DLT에서와 마찬가지로 3D - 2D correspondence를 이용해서 $H$ matrix를 구할 것이다. 

(이 과정은 DLT에서의 과정과 동일한 논리이므로 간략히 수식만 나열하겠다. )

 

$$H = \begin{bmatrix} \color{Red} h_{11} & \color{Red} h_{12} & \color{Red} h_{13} \\ \color{Green} h_{21} & \color{Green} h_{22} & \color{Green} h_{23} \\ \color{Blue} h_{31} & \color{Blue} h_{32} & \color{Blue} h_{33} \end{bmatrix}$$

 

$$H = \begin{bmatrix} \color{Red} H_1^T \\ \color{Green} H_2^T \\ \color{Blue} H_3^T \end{bmatrix}$$

 

$$H_1 = \begin{bmatrix} \color{Red} h_{11} \\ \color{Red} h_{12} \\ \color{Red} h_{13} \end{bmatrix}$$

$$H_2 = \begin{bmatrix} \color{Green} h_{21} \\ \color{Green} h_{22} \\ \color{Green} h_{23} \end{bmatrix}$$

$$H_2 = \begin{bmatrix} \color{Blue} h_{31} \\ \color{Blue} h_{32} \\ \color{Blue} h_{33} \end{bmatrix}$$

 

 

DLT와 마찬가지로 하나의 3D-2D 대응 쌍을 알면 아래와 같이 2개의 식을 얻을 수 있다.

 

$$x= \frac{h_{11}X + h_{12}Y + h_{13}}{h_{31}X + h_{32}Y + h_{33}}$$

$$y = \frac{h_{21}X + h_{22}Y + h_{23}}{h_{31}X + h_{32}Y + h_{33}}$$

 

 

이를 앞에서 정의한 $H_1, H_2, H_3$를 사용해서 아래와 같이 간략하게 표현할 수 있다. 

(이 때 $x_i X_i$는 각각 2D, 3D 점에 해당한다)

 

$$x_i= \frac{H_1^TX_i}{H_3^TX_i}$$

$$y_i = \frac{H_2^TX_i}{H_3^TX_i}$$

 

 

위 식에서 분모를 양변에 곱하고 좌변으로 모든 식을 넘겨서 정리하면 아래와 같이 된다.

 

$$-H_1^TX_i + \quad \quad \ \  + x_i H_3^TX_i = 0$$

$$\quad \quad \ \ -H_2^TX_i + y_iH_3^TX_i = 0$$

 

 

이를 matrix form으로 정리하면 아래와 같다.

 

$$\begin{bmatrix} -X_i^T & 0 & x_iX_i^T \\ 0 & -X_i^T & y_iX_i^T \end{bmatrix}\begin{bmatrix} H_1 \\ H_2 \\ H_3 \end{bmatrix} = 0$$

 

$H$ matrix는 3x3 matrix이지만 homogeneous matrix이므로 8 DoF를 가진다. 

따라서 8개의 식만 있어도 $H$를 구할 수 있다. 

8개의 식을 얻기 위해서는 4개의 대응쌍을 필요로 한다. 

 

4개의 대응쌍에 대해서 위의 matrix form을 완성해주면 아래와 같이 (8x9) (9x1) = 0나타난다. 

 

$$\begin{bmatrix} -X_1^T & 0 & x_1X_1^T \\ 0 & -X_1^T & y_1X_1^T \\ -X_2^T & 0 & x_2X_2^T \\ 0 & -X_2^T & y_2X_2^T \\ ... & ... & ... \\ -X_4^T & 0 & x_4X_4^T \\ 0 & -X_4^T & y_4X_4^T\end{bmatrix}\begin{bmatrix} H_1 \\ H_2 \\ H_3 \end{bmatrix} = 0$$

 

 

위의 식은 DLT에서 살펴보았던 것과 같이 $Ax = 0$ 꼴이므로 coefficient matrix에 SVD를 적용해주면 우리가 구하고자 하는 $H$ matrix를 얻을 수 있다. 

 

 

여기까지는 DLT 방법과 동일했다. 

하지만, Zhang's method에서는 Z = 0라는 가정을 통해 rotation matrix를 간소화 했기 때문에 $H$를 QR decomposition했을 때, $Q, R$ matrix를 대응시켜줄 수 없다. 

 

따라서 QR decomposition을 적용할 수 없고, DLT와는 다른 방법으로 구해야 한다. 

 

(추후 업데이트)

 

 

Implementation

 

Checkerboard image

https://markhedleyjones.com/projects/calibration-checkerboard-collection

 

 

Algorithm