Zhang의 방법으로 카메라 캘리브레이션하기

카메라 캘리브레이션은 3D 복원부터 증강 현실에 이르기까지 많은 컴퓨터 비전 응용에서 기본적이며 필수적인 단계다. 이는 카메라의 내부 파라미터(초점 거리, 주점, 렌즈 왜곡 파라미터 등)와 외부 파라미터(월드 좌표계에 대한 회전과 이동)를 추정하는 과정을 포함한다. 다양한 캘리브레이션 기법 중에서, Zhang의 방법은 평면 체커보드 패턴을 여러 시점에서 관찰하기만 하면 되어 단순하고, 실용성 측면에서도 꽤 효과적이다.

Zhang 방법의 이론 이해하기

Zhang의 방법은 2D 평면 사영 변환인 호모그래피(homography) 개념을 활용한다. (체커보드 같은) 평면 물체를 카메라로 볼 때, 평면 위 점들의 3차원 위치와 이미지에서의 2D 좌표 위치 사이의 관계는 호모그래피 행렬로 기술할 수 있다. \(\mathbf{p}_w = [x, y, z]^T\)를 월드 좌표계에서 체커보드 평면 위 점들의 3D 좌표, \(\mathbf{p}_i = [u, v]^T\)를 그에 대응하는 2D 이미지 좌표라 하자. \(\mathbf{p}_w\)와 \(\mathbf{p}_i\)의 관계는 다음과 같이 표현할 수 있다.

\[s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} = K [R | \mathbf{t}] \begin{bmatrix} x \\ y \\ z \\ 1 \end{bmatrix}\]

여기서 \(s\)는 임의의 스케일 인자, \(K\)는 내부 카메라 행렬, \(R\)은 회전 행렬, \(\mathbf{t}\)는 이동 벡터다. 내부 카메라 행렬 \(K\)는 보통 다음과 같이 정의된다.

\[K = \begin{bmatrix} f_x & \gamma & c_x \\ 0 & f_y & c_y \\ 0 & 0 & 1 \end{bmatrix}\]

여기서 \(f_x\)와 \(f_y\)는 \(x\)축과 \(y\)축 방향의 픽셀 단위 초점 거리, \((c_x, c_y)\)는 주점의 위치 (광축과 이미지 평면의 교점), \(\gamma\)는 비틀림(skew) 파라미터(현대 카메라에서는 보통 0으로 가정)다.

랭크 부족(rank-deficiency) 문제를 피하고 단일 시점에서 카메라 내부 및 외부 파라미터를 유일하게 계산하기 위해, Zhang의 방법은 체커보드 평면의 월드 좌표에 대해 \(z = 0\)을 가정한다. \(z = 0\) 평면 위의 점들에 대해 투영 방정식을 단순화할 수 있다.

\[\begin{align} s \begin{bmatrix} u \\ v \\ 1 \end{bmatrix} &= K [\mathbf{r}_1 | \mathbf{r}_2 | \mathbf{r}_3 | \mathbf{t}] \begin{bmatrix} x \\ y \\ 0 \\ 1 \end{bmatrix} \\&= K [\mathbf{r}_1 | \mathbf{r}_2 | \mathbf{t}] \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} \\&= H \begin{bmatrix} x \\ y \\ 1 \end{bmatrix} \end{align}\]

여기서 \(\mathbf{r}_1, \mathbf{r}_2, \mathbf{r}_3\)은 회전 행렬 \(R\)의 열 벡터다. 행렬 \(H = K [\mathbf{r}_1 \vert \mathbf{r}_2 \vert \mathbf{t}]\)는 \(3 \times 3\) 호모그래피 행렬이다. 이 호모그래피는 평면 위의 2D 월드 좌표 \((x, y)\)와 2D 이미지 좌표 \((u, v)\) 사이의 관계를 정의한다.

선형 방정식을 세우기 위해 다음과 같이 두자.

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

그리고 위 방정식을 전개한다.

\[\begin{align} s u &= x h_{11} + y h_{12} + h_{13}\\ s v &= x h_{21} + y h_{22} + h_{23}\\ s &= x h_{31} + y h_{32} + h_{33}\\ \end{align}\]

스케일 인자 \(s\)를 소거하고 직접 선형 변환(DLT; direct linear transformation)을 사용하면, 선형 방정식 \(A \mathbf{h} = 0\)을 얻는다. 여기서

\[A = \begin{bmatrix} -x & -y & -1 & 0 & 0 & 0 & xu & yu & u \\ 0 & 0 & 0 & -x & -y & -1 & xv & yv & v \\ \end{bmatrix}\]

이고 \(\mathbf{h} = [h_{11}, h_{12}, \ldots, h_{33}]^T\)이다.

Zhang 방법의 핵심은 회전 행렬의 성질을 이용하는 것이다. \(\mathbf{r}_1\)과 \(\mathbf{r}_2\)는 직교하는 단위 벡터이므로 다음 조건을 만족한다.

\[\begin{align} \mathbf{r}_1^T \mathbf{r}_2 &= 0 \\ ||\mathbf{r}_1||^2 &= 1 \\ ||\mathbf{r}_2||^2 &= 1 \\ \end{align}\]

호모그래피의 정의로부터 \(\mathbf{r}_1 = K^{-1}\mathbf{h}_1\), \(\mathbf{r}_2 = K^{-1}\mathbf{h}_2\)(\(\mathbf{h}_1\)과 \(\mathbf{h}_2\)는 호모그래피 행렬 \(H = [\mathbf{h}_1 \vert \mathbf{h}_2 \vert \mathbf{h}_3]\)의 첫 두 열)를 위에 대입하면 다음을 얻는다.

\[\begin{align} \mathbf{h}_1^T K^{-T} K^{-1} \mathbf{h}_2 &= 0 \\ \mathbf{h}_1^T K^{-T} K^{-1} \mathbf{h}_1 &= \mathbf{h}_2^T K^{-T} K^{-1} \mathbf{h}_2 = 1 \end{align} \tag{1}\label{eq:homography}\]

여기서 \(B = K^{-T} K^{-1}\)이라 하자. 이 행렬 \(B\)는 대칭이고 양의 정부호(positive definite)다. 따라서 \(B\)는 여섯 개의 파라미터로 정의할 수 있다. 위 두 식은 관찰된 각 체커보드 이미지에 대해 내부 파라미터에 대한 두 개의 제약을 제공한다. 적어도 세 가지 다른 방향에서 이미지를 촬영하면, \(B\)의 원소를 풀기 위한 선형 방정식 시스템을 세울 수 있다(6개의 파라미터와 6개의 방정식). \(B\)가 결정되면, 내부 행렬 \(K\)는 콜레스키 분해(Cholesky decomposition)를 통해 구할 수 있다.

\(K\)를 얻은 후, 각 시점의 외부 파라미터(\(R\)과 \(\mathbf{t}\))는 호모그래피 행렬의 열을 이용해 계산할 수 있다.

\[\begin{align} \mathbf{r}_1 &= \lambda K^{-1} \mathbf{h}_1 \\ \mathbf{r}_2 &= \lambda K^{-1} \mathbf{h}_2 \\ \mathbf{r}_3 &= \mathbf{r}_1 \times \mathbf{r}_2 \\ \mathbf{t} &= \lambda K^{-1} \mathbf{h}_3 \\ \end{align}\]

여기서 \(\lambda = 1 / \vert\vert K^{-1} \mathbf{h}_1 \vert\vert\)는 \(\mathbf{r}_1\)과 \(\mathbf{r}_2\)가 단위 벡터가 되도록 하는 스케일링 인자다.

마지막으로, 추정된 파라미터를 정밀하게 계산하고 렌즈 왜곡(반경 방향 및 접선 방향)을 고려하기 위해 비선형 최적화(예: Levenberg-Marquardt)를 수행한다. 이 최적화는 재투영 오차(reprojection error), 즉 관찰된 이미지 점과 현재 카메라 파라미터로 이미지 평면에 투영된 점 사이의 거리 제곱의 합을 최소화한다.

MATLAB 구현

제공된 MATLAB 코드는 Zhang 방법의 실용적인 구현을 보여준다. 주요 함수들을 살펴보자.

체커보드 점 검출

MATLAB의 Computer Vision Toolbox의 일부인 detectCheckerboardPoints 함수는 입력 이미지 집합에서 체커보드 코너를 자동으로 검출하는 데 사용된다. 이 코너들의 픽셀 좌표 imagePoints와 체커보드 격자의 크기 boardSize를 반환한다.

호모그래피 계산

computeHomography 함수는 체커보드 평면 위의 3D 월드 점(P_w)을 2D 이미지 점(P_i)으로 매핑하는 호모그래피 행렬 H를 계산한다. 호모그래피 행렬의 원소들을 담은 벡터 h에 대한 선형 방정식 시스템 A*h = 0을 세운 뒤, 특이값 분해(SVD)를 이용해 h를 푼다.

function H = computeHomography(P_w, P_i)
    numPoints = size(P_w,1);
    A = zeros(2*numPoints, 9);

    for i = 1:numPoints
        X = P_w(i,1);
        Y = P_w(i,2);
        u = P_i(i,1);
        v = P_i(i,2);

        A(2*i-1,:) = [-X, -Y, -1,  0,  0,  0, X*u, Y*u, u];
        A(2*i,  :) = [ 0,  0,  0, -X, -Y, -1, X*v, Y*v, v];
    end

    [~,~,V] = svd(A);
    H = reshape(V(:,end), [3,3])'; % the smallest singular value
end

내부 파라미터 계산

메인 스크립트의 이 부분은 호모그래피 행렬을 기반으로 행렬 B를 풀기 위한 선형 방정식 시스템을 세운다.

식 \eqref{eq:homography}에서 \(B = K^{-T} K^{-1}\)로 두었다. 행렬 \(B\)는 대칭이고 양의 정부호이므로, \(B\)의 독립적인 원소 수는 여섯 개다. 따라서 다음과 같이 둔다.

\[B = \begin{bmatrix} b_1 & b_2 & b_4 \\ b_2 & b_3 & b_5 \\ b_4 & b_5 & b_6 \\ \end{bmatrix}.\]

\(\mathbf{h}_1 = [h_{11}, h_{21}, h_{31}]^T\)와 \(\mathbf{h}_2 = [h_{12}, h_{22}, h_{32}]^T\)를 이용해, 식 \eqref{eq:homography}의 제약을 \(b_*\)에 대해 전개할 수 있다.

\[\begin{align} & \mathbf{h}_i^T B \mathbf{h}_j \\ &= (h_{1i} h_{1j}) b_1 \\&+ (h_{1i} h_{2j} + h_{2i} h_{1j}) b_2 \\&+ (h_{2i} h_{2j}) b_3 \\&+ (h_{3i} h_{1j} + h_{1i} h_{3j}) b_4 \\&+ (h_{2i} h_{3j} + h_{3i} h_{2j}) b_5 \\&+ (h_{3i} h_{3j}) b_6, \end{align}\]

여기서 \(i,j \in \{1,2\}\)이다. 그러면 각 호모그래피 행렬에 대해 두 개의 제약이 충족되어야 하므로, 총 \(N\)개의 호모그래피에 대해 \(V \mathbf{b} = \mathbf{0}\)을 만족한다. 여기서 \(V \in \mathbb{R}^{2N \times 6}\)이고 \(\mathbf{b} = [b_1, b_2, \ldots, b_6]^T\)는 \(V\)의 영공간(null space)이다.

v() 인라인 함수는 위 방정식들을 세우는 데에 사용하였다. SVD를 통해 b(B의 벡터화된 형태)를 얻으면 B를 복원하고, 콜레스키 분해를 이용해 내부 행렬 K를 유도한다.

% b11 b12 b22 b13 b23 b33
v = @(H,i,j) [
    H(1,i)*H(1,j), ...
    H(1,i)*H(2,j) + H(2,i)*H(1,j), ...
    H(2,i)*H(2,j), ...
    H(3,i)*H(1,j) + H(1,i)*H(3,j), ...
    H(3,i)*H(2,j) + H(2,i)*H(3,j), ...
    H(3,i)*H(3,j)
];

for i = 1:numImages
    H{i} = computeHomography(objectPoints, imagePoints(:,:,i));

    v1 = v(H{i}, 1, 2);
    v2 = v(H{i}, 1, 1) - v(H{i}, 2, 2);

    V = [V; v1; v2];
end

[~,~,Vmat] = svd(V);
b = Vmat(:,end);
B = [b(1), b(2), b(4); b(2), b(3), b(5); b(4), b(5), b(6)];
if det(B) < 0
    B = -B;
end
Kinv = chol(B,'upper');
K = inv(Kinv/Kinv(3,3));

외부 파라미터 계산

초기 내부 행렬 K로, 각 시점의 회전 행렬 R과 이동 벡터 t를 호모그래피 행렬로부터 계산한다.

R = {};
t = {};
for i = 1:numImages
    r1 = K\H{i}(:,1);
    r2 = K\H{i}(:,2);

    n = norm(r1);
    r1 = r1 / norm(r1);
    r2 = r2 / norm(r2);

    r3 = cross(r1, r2);
    r3 = r3 / norm(r3);

    R{i} = [r1 r2 r3];
    t{i} = K\H{i}(:,3);
    t{i} = t{i} / n;

    if t{i}(3) < 0 % the plane lies behind the camera
        t{i} = -t{i};
        R{i} = -R{i};
    end
end

재투영 오차 계산

이 함수는 비선형 최적화의 핵심이다. 현재 카메라 파라미터(내부와 외부)를 받아 3D objectPoints를 이미지 평면에 투영한 뒤, 이 투영된 점과 관찰된 imagePoints 사이의 차이를 계산한다. 이 차이가 최적화 알고리즘이 최소화하려는 재투영 오차다.

function err = reprojectionError(params, objectPoints, imagePoints)
    [fx, fy, s, cx, cy, distCoeffs, R, t] = parse_params(params);
    
    K = [fx, s, cx; 0, fy, cy; 0, 0, 1];

    err = [];
    numPoints = size(objectPoints,1);
    numImages = (length(params)-10)/7;

    for i = 1:numImages
        imagePoint = imagePoints(:,:,i);
        X = [R{i} t{i}; 0 0 0 1] * [objectPoints, ones(size(objectPoints,1),1)]';
        projected = X(1:3,:) ./ X(3,:);
        projected = K * distort(projected, distCoeffs);
        projected = projected(1:2,:)';

        % compute reprojection error
        err = [err; imagePoint(:) - projected(:)];
    end
end

내부 및 외부 파라미터 최적화와 왜곡 파라미터 계산

이 함수들은 lsqnonlin(Levenberg-Marquardt 알고리즘)을 이용한 비선형 최적화를 활용한다. 초기 파라미터를 받아 재투영 오차를 최소화함으로써 반복적으로 최적화한다. variations 파라미터는 최적화 중 파라미터 탐색의 범위를 제어하여 단계적 최적화 과정을 가능하게 한다.

function params_opt = myEstimateCameraParameters(K, R, t, objectPoints, imagePoints)

    params_init = [K(1,1), K(2,2), K(1,2), K(1,3), K(2,3), zeros(1,5)]; % fx, fy, s, cx, cy, distCoeffs
    for i = 1:length(R)
        quat = rotm2quat(R{i});
        params_init = [params_init, quat, t{i}'];
    end

    options = optimoptions('lsqnonlin', 'Display', 'final-detailed', 'FunctionTolerance', 1e-10, 'StepTolerance', 1e-20, 'MaxIterations', 2000, 'Algorithm', 'levenberg-marquardt');

    variations = [0.2, 0.2, 0.2, 0.2, 0.2, Inf, Inf, 0, 0, 0, Inf*ones(1,7*length(R))];
    params_opt = optimize(params_init, objectPoints, imagePoints, variations, options);

    options = optimoptions('lsqnonlin', 'Display', 'final-detailed', 'FunctionTolerance', 1e-20, 'StepTolerance', 1e-20, 'MaxIterations', 2000, 'Algorithm', 'levenberg-marquardt');

    variations = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, Inf, Inf, 0, 0.2*ones(1,7*length(R))];
    params_opt = optimize(params_opt, objectPoints, imagePoints, variations, options);

    options = optimoptions('lsqnonlin', 'Display', 'final-detailed', 'FunctionTolerance', 1e-30, 'StepTolerance', 1e-20, 'MaxIterations', 2000, 'Algorithm', 'levenberg-marquardt');

    variations = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, Inf, 0.2*ones(1,7*length(R))];
    params_opt = optimize(params_opt, objectPoints, imagePoints, variations, options);
end

function params_opt = optimize(params_init, objectPoints, imagePoints, variations, options)

params_lower = params_init .* (1 - variations);
params_upper = params_init .* (1 + variations);

idx = params_init < 0;
tmp = params_upper(idx);
params_upper(idx) = params_lower(idx);
params_lower(idx) = tmp;

idx = params_init == 0;
params_upper(idx) = variations(idx);
params_lower(idx) = -variations(idx);

params_opt = lsqnonlin(@(params) reprojectionError(params, objectPoints, imagePoints), params_init, params_lower, params_upper, options);

초기 파라미터 추정을 위한 선형 해법과 왜곡 파라미터를 위한 비선형 최적화를 결합함으로써, Zhang의 방법은 카메라를 강건하고 정확하게 캘리브레이션하는 방법을 제공한다. 이는 컴퓨터 비전 작업에서 정밀한 측정과 정확한 공간 이해에 매우 중요하다.