Search
Duplicate

AI/ Image Rectification

Image Rectification은 두 이미지 간의 대응점을 찾기 위해 시차 계산을 단순화하는 것이다. 이것은 두 카메라가 수평으로 배치된 경우 두 이미지의 행을 일치시키는 것이 되고, 두 카메라가 수직으로 배치된 경우 두 이미지의 열을 일치시키는 것에 해당한다.
카메라가 이동하거나 회전한다면 epipolar geometry가 변경되기 때문에 Image Rectification을 매번 다시 수행해야 한다.

epipolar geometry를 이용한 Image Rectification

epipolar geometry framework에서 두 카메라가 평행하다고 가정하자. 이때 두 카메라가 동일한 KK를 가졌다고 하면, 상대적 회전이 없고 (R=IR= I), xx축을 따라 translation만 존재하므로 T=(Tx,0,0)T = (T_x, 0, 0) Essential matrix EE는 다음과 같이 정의할 수 있다.
E=[T×]R=[00000−Tx0Tx0]E = [T_\times]R = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -T_x \\ 0 & T_x & 0 \end{bmatrix}
EE가 알려지면, 이미지 평면의 점과 연관된 epipolar line의 방향을 찾을 수 있다. 점 p′p'와 연관된 epipolar line ℓ\ell의 방향은 다음과 같이 계산된다.
ℓ=Ep′=[00000−Tx0Tx0][u′v′1]=[0−TxTxv′]\ell = Ep' = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & -T_x \\ 0 & T_x & 0 \end{bmatrix} \begin{bmatrix} u' \\ v' \\ 1 \end{bmatrix} = \begin{bmatrix} 0 \\ -T_x \\ T_x v' \end{bmatrix}
이에 대해 epipolar 제약조건 p⊤Ep′=0p^\top Ep' = 0을 적용하면 아래와 같고,
p⊤Ep′=[uv1][0−TxTxv′]=0⋅u+−Txv+Txv′=0p^\top Ep' = \begin{bmatrix} u & v & 1 \end{bmatrix} \begin{bmatrix} 0 \\ -T_x \\ T_x v' \end{bmatrix} = 0 \cdot u + -T_x v + T_xv' = 0
v=v′v=v'가 되어 pp와 p′p'가 동일한 vv-좌표를 공유한다는 사실을 알 수 있다. 결과적으로 대응점 사이에 매우 간단한 관계가 존재한다. 그러므로 주어진 두 이미지를 평행으로 만드는 과정인 rectification은 stereo matching을 위한 이미지에서 대응점 사이의 관계를 식별할 때 유용하다.

epipole 추정

이미지 쌍을 rectifying 하는데 두 카메라 행렬 K,K′K, K' 이나 그들 사이의 상대적 변환 R,TR, T에 대한 지식이 필요하지 않다. 대신 Normalized Eight Point Algorithm에 의해 추정된 Fundamental matrix를 사용할 수 있다. Fundamental matrix를 얻으면 각 대응점 pip_i와 pi′p_i'에 대한 epipolar line ℓi\ell_i와 ℓi′\ell_i'를 계산할 수 있다.
모든 epipolar line이 epipole을 교차하므로 epipolar line을 구하면, epipole을 ee과 e′e'를 추정할 수 있다. 그러나 실제로는 noisy 측정 때문에 epipolar line이 단일점에서 교차하지 않는다. 따라서 epipole을 찾기 위해 least square를 사용해야 한다.
각 epipolar line은 벡터 ℓ\ell로 표현할 수 있고, 직선 위의 모든 점은 {x∣ℓ⊤x=0}\{x|\ell^\top x = 0\}를 만족한다. 각 epipolar line을 직선의 방정식의 homogeneous 표현인 ℓi=[ℓi,1ℓi,2ℓi,3]⊤\ell_i = [\ell_{i,1} \quad \ell_{i,2} \quad \ell_{i,3}]^\top로 정의한 다음 선형 시스템의 방정식을 세우고 SVD를 사용하여 epipole ee를 찾을 수 있다.
[ℓ1⊤⋮ℓn⊤]e=0\begin{bmatrix} \ell_1^\top \\ \vdots \\ \ell_n^\top \end{bmatrix} e = 0
그러나 이렇게 찾은 e,e′e, e'는 무한대의 점이 아니다. 따라서 e,e′e, e'를 horizontal axis를 따라 무한대로 매핑하는 별도의 homography H1,H2H_1, H_2를 찾아야 한다.

homography 추정

Homography H2H_2 추정

우선 homography H2H_2을 찾는 것에서 시작하면, 두 번째 이미지에 대해 homogeneous 좌표에서 중심이 (0,0,1)(0, 0, 1)이 되도록 translate 하는 translation matrix를 정의한다.
T=[10−width201−height2001]T = \begin{bmatrix} 1 & 0 & -{\text{width} \over 2} \\ 0 & 1 & -{\text{height} \over 2} \\ 0 & 0 & 1 \end{bmatrix}
translation을 적용한 후에 epipole을 horizontal axis 상에 어떤 점 (f,0,1)(f, 0, 1)에서 배치하도록 rotation을 적용한다. 평행이동 된 epipole Te′Te'가 homogeneous 좌표 (e1′,e2′,1)(e_1', e_2', 1)에 위치한다면 적용되는 rotation은 다음과 같다.
R=[αe1′e1′2+e2′2αe2′e1′2+e2′20−αe2′e1′2+e2′2αe1′e1′2+e2′20001]R = \begin{bmatrix} \alpha{e_1' \over \sqrt{e_1'^2 + e_2'^2}} & \alpha{e_2' \over \sqrt{e_1'^2 + e_2'^2}} & 0 \\ -\alpha{e_2' \over \sqrt{e_1'^2 + e_2'^2}} & \alpha{e_1' \over \sqrt{e_1'^2 + e_2'^2}} & 0 \\ 0 & 0 & 1 \end{bmatrix}
여기서 e1′≥0e_1' \ge 0이면 α=1\alpha = 1이고 그렇지 않으면 α=−1\alpha = -1이다.
rotation을 적용한 후에 (f,0,1)(f, 0, 1)의 임의의 점을 horizontal axis 상의 무한대 점 (f,0,0)(f, 0, 0)으로 가져오려면 다음의 변환만 적용하면 된다.
G=[100010−1f01]G = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -{1\over f} & 0 & 1 \end{bmatrix}
이 변환을 적용한 후에 마침내 무한대에 있는 epipole을 갖게 된다.
최종적으로 T−1T^{-1}을 이용해서 homogeneous 좌표를 regular 이미지 공간으로 다시 translate 한다. 따라서 두 번째 이미지를 rectify 하기 위해 적용하는 homography H2H_2는 다음 형식이 된다.
H2=T−1GRTH_2 = T^{-1}GRT

Homography H1H_1 추정

H2H_2를 찾으면 첫 번째 이미지에 대한 H1H_1을 찾을 수 있다. 이것은 이미지의 대응점 사이의 제곱 거리의 합을 최소화하는 방식으로 찾을 수 있다.
arg min⁡H1∑i∥H1pi−H2pi′∥2\argmin_{H_1} \sum_i \|H_1p_i - H_2p_i'\|^2
유도 없이 H1H_1이 다음 형식을 갖는다는 것을 보인다.
H1=HAH2MH_1 = H_A H_2 M
여기서 Fundamental Matrix F=[e]×MF = [e]_\times M이고
HA=[a1a2a3010001]H_A = \begin{bmatrix} a_1 & a_2 & a_3 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}
여기서 (a1,a2,a3)(a_1, a_2, a_3)은 나중에 계산될 특정한 벡터 a\bold{a}의 성분을 구성한다.
MM을 알기 전에 임의의 3×33 \times 3 skew-symmetric matrix AA에 대해 A=A3A = A^3(scale까지)라는 것을 이용한다. 임의의 cross product matrix [e]×[e]_\times가 skew-symmetric이고 Fundamental matrix FF는 scale까지만 알 수 있으므로
F=[e]×M=[e]×[e]×[e]×M⏟F=[e]×[e]×FF = [e]_\times M = [e]_\times [e]_\times \underbrace{[e]_\times M}_{F} = [e]_\times [e]_\times F
[e]×=[e]×3[e]_\times = [e]_\times^3이므로 [e]×[e]×[e]×M=[e]×[e]×F[e]_\times[e]_\times[e]_\times M = [e]_\times[e]_\times F를 다음과 같이 정리할 수 있다.
M=[e]×FM = [e]_\times F
MM의 열에 ee의 임의의 스칼라 곱을 추가해도 F=[e]×MF = [e]_\times M은 여전히 scale까지 성립한다. 그러므로 임의의 벡터 vv에 대해 MM을 다음과 같이 더 일반적인 형식으로 정의할 수 있다.
M=[e]×F+ev⊤M = [e]_\times F + ev^\top
v⊤=[111]v^\top = [1 \quad 1 \quad 1]로 설정하면 MM을 매우 잘 정의할 수 있다.
H2H_2와 MM의 값을 알고 있으므로, H1H_1의 최소화 식에 대해 p^i=H2Mpi\hat{p}_i = H_2Mp_i와 p^i′=H2pi′\hat{p}_i' = H_2p_i'로 치환하면 최소화 문제는 다음이 된다.
arg min⁡H1∑i∥H1pi−H2pi′∥2=arg min⁡HAH2M∑i∥HAH2Mpi⏟p^i−H2pi′⏟p^i′∥2=arg min⁡HA∑i∥HAp^i−p^i′∥2\begin{aligned} \argmin_{H_1} \sum_i \|H_1p_i - H_2p_i'\|^2 &=\argmin_{H_AH_2M} \sum_i \|H_A \underbrace{H_2Mp_i}_{\hat{p}_i} - \underbrace{H_2p_i'}_{\hat{p}_i'}\|^2 \\ &= \argmin_{H_A} \sum_i \|H_A\hat{p}_i - \hat{p}_i'\|^2 \end{aligned}
p^i=(x^i,y^i,1)\hat{p}_i = (\hat{x}_i, \hat{y}_i, 1)과 p^i′=(x^i′,y^i′,1)\hat{p}_i' = (\hat{x}_i', \hat{y}_i', 1)를 설정하면 최소화 문제는 다음과 같이 바꿀 수 있다.
arg⁡min⁡a∑i(a1x^i+a2y^i+a3−x^i′)2+(y^i−y^i′)2\arg \min_\bold{a} \sum_i (a_1 \hat{x}_i + a_2 \hat{y}_i + a_3 - \hat{x}_i')^2 + (\hat{y}_i - \hat{y}_i')^2
y^i−y^i′\hat{y}_i - \hat{y}_i'이 상수 값이므로 최소화 문제는 추가로 다음으로 축소된다.
arg⁡min⁡a∑i(a1x^i+a2y^i+a3−x^i′)2\arg \min_\bold{a} \sum_i (a_1 \hat{x}_i + a_2 \hat{y}_i + a_3 - \hat{x}_i')^2
궁극적으로 이것은 a\bold{a}에 대한 least-square 문제 Wa=bW\bold{a} = b를 해결하는 것으로 분해될 수 있다.
W=[x^1y^11⋮x^ny^n1],b=[x^1′⋮x^n′]W = \begin{bmatrix} \hat{x}_1 & \hat{y}_1 & 1 \\ & \vdots \\ \hat{x}_n & \hat{y}_n & 1 \end{bmatrix}, b = \begin{bmatrix} \hat{x}_1' \\ \vdots \\ \hat{x}_n' \end{bmatrix}
최종적으로 a\bold{a}를 계산한 후에 HAH_A를 계산하고 마지막으로 H1H_1을 구할 수 있다.

Rectify 절차

최종적으로 이미지 Rectify 절차는 다음과 같다.
1.
Normalized Eight Point Algorithm을 이용하여 Fundamental matrix 계산
2.
Fundamental matrix 이용하여 각 대응점 pip_i와 pi′p_i'에 대한 epipolar line ℓi\ell_i와 ℓi′\ell_i'를 계산
3.
각 epipolar line ℓ,ℓ′\ell, \ell'을 이용하여 선형 방정식을 세우고 SVD를 해결하여 epipole e,e′e, e' 계산
4.
e′e'를 무한대로 매핑하는 homography H2H_2 계산
a.
두 번째 이미지에 대해 homogeneous 좌표에서 중심이 (0,0,1)(0, 0, 1)이 되도록 translate하는 translation matrix TT 적용
b.
epipole을 horizontal axis 상에 어떤 점 (f,0,1)(f, 0, 1)에서 배치하도록 rotation 하는 rotation matrix RR 적용
c.
rotation 적용 후 (f,0,1)(f, 0, 1)의 점을 horizontal axis 상의 무한대 점 (f,0,0)(f, 0, 0)으로 매핑하는 행렬 GG 적용
d.
최종적으로 이미지 좌표로 되돌리도록 translation matrix의 역행렬 T−1T^{-1} 적용
5.
H2H_2를 이용하여 ee을 무한대로 매핑하는 homography H1H_1 계산
a.
이미지의 대응점 사이의 제곱 거리의 합을 최소화하는 식 정의. 여기서 H1=HAH2MH_1 = H_A H_2 M 로 정의된다.
b.
Fundamental Matrix FF를 이용하여 MM 계산. M=[e]×F+ev⊤M = [e]_\times F + ev^\top
c.
H2H_2와 MM을 이용하여 최소화 문제로 HAH_A 계산
d.
최종적으로 구해진 H2,M,HAH_2, M, H_A를 이용하여 H1H_1 계산

카메라 Intrinsic과 Extrinsic을 이용한 Image Rectification

만일 카메라의 Intrinsic과 Extrinsic과 왜곡이 알려져 있다면, Fundamental Matrix와 Epipolar Line을 찾을 필요 없이 Intrinsic과 Extrinsic을 이용하여 Rectify를 직접 수행할 수 있다.

Rotation, Projection, Disparity-to-depth mapping 행렬 계산

우선 rectify를 수행하기 위한 rotation, projection, disparity-to-depth mapping 행렬을 계산해야 한다.
K1,K2K_1, K_2를 두 카메라의 카메라 행렬, D1,D2D_1, D_2를 두 카메라의 distortion coefficients, R,TR, T를 두 카메라의 extrinsic(Rotation과 Translation)이라 할 때, 다음의 단계를 통해 rectification을 위한 rotation 행렬과 projection 행렬과 disparity-to-depth 매핑 행렬을 구할 수 있다.
1.
두 카메라의 이미지 중심점 계산
cx=cx1+cx22cy=cy1+cy22c_x = {c_{x1} + c_{x2} \over 2} \\ c_y = {c_{y1} + c_{y2} \over 2}
2.
두 카메라의 rotation matrix 계산
첫 번째 카메라의 회전 행렬 R1R_1을 단위 행렬 II로 설정한다.
R1=IR_1 = I
두 번째 카메라의 회전 행렬 R2R_2을 두 카메라 간의 회전 행렬 RR과 첫 번째 카메라의 회전 행렬 R1R_1을 곱한 것으로 설정한다.
R2=R⋅R1R_2 = R \cdot R_1
3.
두 카메라의 projection matrix 계산
첫 번째 카메라의 projection matrix P1P_1를 카메라 행렬 K1K_1과 2에서 구한 rotation 행렬 R1R_1을 이용하여 다음처럼 정의한다.
P1=K1⋅[R1∣0]P_1 = K_1 \cdot [R_1 |0]
두 번째 카메라의 projection matrix P2P_2를 카메라 행렬 K2K_2과 2에서 구한 rotation 행렬 R2R_2과 두 카메라 사이의 translation vector TT를 이용하여 다음처럼 정의한다.
P2=K2⋅[R2∣T]P_2 = K_2 \cdot [R_2 |T]
4.
Disparity-to-depth mapping matrix 계산
1에서 구한 새로운 이미지 중심점과 focal length ff와 두 카메라 사이의 baseline BB을 이용하여 Disparity-to-depth mapping 행렬 QQ을 다음처럼 정의한다.
Q=[100−cx010−cy000f00−1B0]Q = \begin{bmatrix} 1 & 0 & 0 & -c_x \\ 0 & 1 & 0 & -c_y \\ 0 & 0 & 0 & f \\ 0 & 0 & -{1\over B} & 0 \end{bmatrix}

Rectification 계산

다음으로 위에서 계산된 Rotation 행렬과 Projection 행렬을 이용하여 Rectification Map을 생성한다. 이것은 그 자체로 이미지는 아니고, 원본 이미지를 Rectified 이미지로 변환하는데 사용되는 매핑 행렬이다.
1.
왜곡 보정
카메라 intrinsic의 radial distortion 계수 k1,k2,k3k_1, k_2, k_3과 tangential distortion 계수 p1,p2p_1, p_2를 사용하여 다음과 같이 Radial Distortion 보정 xr,yrx_r, y_r와 Tangential Distortion 보정 xt,ytx_t, y_t을 수행한다.
xr=x(1+k1r2+k2r4+k3r6)yr=y(1+k1r2+k2r4+k3r6)xt=x+[2p1xy+p2(r2+2x2)]yt=y+[p1(r2+2y2)+2p2xy]\begin{aligned} x_r &= x(1 + k_1r^2 + k_2r^4 + k_3r^6) \\ y_r &= y(1 + k_1r^2 + k_2r^4 + k_3r^6) \\ x_t &= x + [2p_1xy + p_2(r^2 + 2x^2)] \\ y_t &= y + [p_1(r^2 + 2y^2) + 2p_2xy] \end{aligned}
여기서 r2=x2+y2r^2 = x^2 + y^2.
2개의 보정값을 구했으면 원래의 값 x,yx, y에 대해 다음처럼 보정된 값 xundistorted,yundistortedx_\text{undistorted}, y_\text{undistorted}을 구한다.
xundistorted=xr+xt−xyundistorted=yr+yt−y\begin{aligned} x_\text{undistorted} &= x_r + x_t - x \\ y_\text{undistorted} &= y_r + y_t - y \end{aligned}
2.
rectification 변환
왜곡 보정을 수행한 xu,yux_u, y_u에 대해 Rotation 행렬을 적용하여 새로운 좌표를 구한다.
[x′y′1]=R[xuyu1]\begin{bmatrix} x' \\ y' \\ 1\end{bmatrix} = R \begin{bmatrix} x_u \\ y_u \\ 1\end{bmatrix}
3.
projection 수행
rectification 변환을 수행한 결과에 대해 Projection 행렬을 적용하여 최종 rectified 좌표를 구한다.
[uvw]=P[x′y′1]\begin{bmatrix} u \\ v \\ w \end{bmatrix} = P \begin{bmatrix} x' \\ y' \\ 1\end{bmatrix}
최종적으로 uu와 vv를 ww로 나누어 정규화한다.
(u,v)=(uw,vw)(u, v) = \left({u \over w}, {v \over w} \right)
이렇게 얻은 u,vu, v는 기존 이미지의 x,yx, y 좌표에 해당하는 rectified 이미지의 픽셀 좌표가 된다. 그러나 이 값은 정수가 아니라 실수값이기 때문에 좌표값으로 그대로 사용할 수 없다. 따라서 다음과 같이 후처리를 통해 최종적으로 픽셀 값을 생성한다.

후처리 Interpolation

rectified 된 좌표의 값을 그대로 사용하지 않고 해당 좌표의 인접 픽셀과 보간하여 최종 Rectified Image를 생성한다.
1.
좌표 클리핑
rectified를 통해 얻은 u,vu, v 좌표가 원래 이미지의 영역을 벗어날 수 있기 때문에 다음과 같이 clipping한다.
x=max⁡(0,min⁡(u,cols−1))y=max⁡(0,min⁡(v,rows−1))\begin{aligned} x &= \max(0, \min(u, \text{cols}-1)) \\ y &= \max(0, \min(v, \text{rows}-1)) \end{aligned}
2.
정수와 소수 부분 분리
u,vu, v가 실수값이기 때문에 그대로 좌표로 사용할 수 없으므로 다음과 같이 정수와 소수 부분을 분리한다.
x1=⌊u⌋y1=⌊v⌋x2=x1+1y2=y1+1α=u−x1β=v−y1\begin{aligned} x_1 &= \lfloor u \rfloor \\ y_1 &= \lfloor v \rfloor \\ x_2 &= x_1 + 1 \\ y_2 &= y_1 + 1 \\ \alpha &= u - x_1 \\ \beta &= v - y_1 \end{aligned}
3.
4개 인접 픽셀 값 가져오기
분리한 정수 부분을 이용하여 인접 픽셀의 값을 가져온다. u,vu, v의 소수 부분이 0에 가까웠다면 I11I_{11}이 가장 가까운 픽셀에 해당하고, 소수 부분이 1에 가까웠다면 I22I_{22}가 가장 가까운 픽셀에 해당한다.
I11=src(y1,x1)I12=src(y2,x1)I21=src(y1,x2)I22=src(y2,x2)\begin{aligned} I_{11} &= \text{src}(y_1, x_1) \\ I_{12} &= \text{src}(y_2, x_1) \\ I_{21} &= \text{src}(y_1, x_2) \\ I_{22} &= \text{src}(y_2, x_2) \end{aligned}
4.
선형 보간
α,β\alpha, \beta는 소수 부분에 해당하므로, 4개의 인접 픽셀을 합하는데, 비중을 조절하는 역할을 한다.
I(u,v)=(1−α)(1−β)I11+α(1−β)I21+(1−α)βI12+αβI22I(u, v) = (1-\alpha)(1-\beta)I_{11} + \alpha(1-\beta)I_{21} + (1-\alpha)\beta I_{12} + \alpha\beta I_{22}

Sample Code

위의 절차는 opencv에서 stereoRectify(), initUndistortRectifyMap(), remap() 함수로 구현되어 있다. 아래 코드 참조. stereoRectify()는 두 카메라에 대해 한 번만 수행하지만, initUndistortRectifyMap()와 remap()은 두 카메라에 대해 각각 수행한다는 것에 유의.
#include <opencv2/opencv.hpp> #include <iostream> using namespace cv; using namespace std; int main() { // 이미지 캡처 (예제 이미지 사용) Mat img1 = imread("left_image.jpg", IMREAD_GRAYSCALE); Mat img2 = imread("right_image.jpg", IMREAD_GRAYSCALE); if (img1.empty() || img2.empty()) { cerr << "Could not open or find the images!\n"; return -1; } // 카메라 캘리브레이션 파라미터 (예제 데이터) Mat K1 = (Mat_<double>(3, 3) << 700, 0, 320, 0, 700, 240, 0, 0, 1); Mat K2 = (Mat_<double>(3, 3) << 700, 0, 320, 0, 700, 240, 0, 0, 1); Mat D1 = Mat::zeros(1, 5, CV_64F); // 왜곡 계수 Mat D2 = Mat::zeros(1, 5, CV_64F); // 왜곡 계수 Mat R = (Mat_<double>(3, 3) << 0.9998, -0.0176, 0.0045, 0.0176, 0.9998, -0.0045, -0.0045, 0.0045, 1); // 예제 회전 행렬 Mat T = (Mat_<double>(3, 1) << 0.1, 0, 0); // 두 카메라 간의 거리 (baseline) // Rectification 변환 행렬 계산 Mat R1, R2, P1, P2, Q; stereoRectify(K1, D1, K2, D2, img1.size, R, T, R1, R2, P1, P2, Q); // Rectification 맵 생성 // initUndistortRectifyMap()의 output은 2개인데, distortion이 보정된 x값 행렬과 y값 행렬에 해당한다. Mat map1x, map1y, map2x, map2y; initUndistortRectifyMap(K1, D1, R1, P1, imageSize, CV_32FC1, map1x, map1y); initUndistortRectifyMap(K2, D2, R2, P2, imageSize, CV_32FC1, map2x, map2y); // 이미지 보정 및 변환 Mat rectifiedImage1, rectifiedImage2; remap(img1, rectifiedImage1, map1x, map1y, INTER_LINEAR); remap(img2, rectifiedImage2, map2x, map2y, INTER_LINEAR); // 결과 출력 imshow("Rectified Left Image", rectifiedImage1); imshow("Rectified Right Image", rectifiedImage2); waitKey(0); return 0; }
C++
복사

참고

Made with