Triangulation

Triangulation

2021, Feb 15    


  • 참조 : https://www.cs.cmu.edu/~16385/s17/Slides/



목차



Triangulation 개념


Drawing


  • 위 표의 3가지 항목인 Pose Estimation, Triangulation, Reconstruction 중 이번 글에서 다룰 내용은 Triangulation이며 이 태스크는 도표의 분류와 같이 Pose Estimation 또는 Reconstruction과 차이가 있음을 알 수 있습니다. 그 내용에 대하여 좀 더 자세하게 살펴보겠습니다.


Drawing


  • Triangulation은 2개의 이미지에서 같은 지점을 가리키는 2D point 셋과 카메라 파라미터가 있을 때, 2D point의 실제 3D point를 추정하는 작업을 의미합니다.
  • 2D Point는 보통 Feature Extraction & Matching 작업 등과 같은 방식을 통하여 얻기 때문에 일부 오차가 있음을 가정합니다. 또한 Triangulation에서는 이미지의 왜곡 보정이 반영되어 있음을 가정합니다.
  • 이번 글에서는 다음과 같이 용어 및 기호를 통일하여 사용하도록 하겠습니다.


  • Given a set of (noisy) matched points: x=(u,v),x=(u,v)Given a set of (noisy) matched points: x=(u,v),x=(u,v)(1)
  • Camera projection matrices: P,PCamera projection matrices: P,P(2)
  • Estimated 3D points: XT=[XYZ]Estimated 3D points: XT=[XYZ](3)


  • 식 (1)에서 x,xx,x 각각은 이미지 상에서 서로 같은 지점의 이미지 좌표를 의미합니다.
  • 식 (2)에서 PPCamera Extrinsic MatrixCamera Intrinsic Matrix를 모두 반영한 Camera Projection Matrix를 의미합니다. 상세 내용은 글 상단의 카메라 캘리브레이션 링크를 참조해 보시기 바랍니다.
  • 식 (3)의 XX 는 식 (1)의 이미지 2D points인 x,xx,x 가 이미지에 투영되기 이전의 3D points를 의미하며 최종 추정하고자 하는 값이 됩니다.


  • x=PXx=PX(4)


  • 식 (4)에서 3D pointsXXCamera Projection Matrix를 통해 이미지에 투영하면 xx 가 됩니다.
  • 이 때, XX 를 알고 싶으나 XX 의 후보가 되는 값은 무수히 많이 있습니다. 관련 내용은 아래 링크를 참조하시면 됩니다.
    • 동차좌표계 : https://gaussian37.github.io/vision-concept-homogeneous_coordinate/


Drawing


  • 위 그림과 같이 빨간색 선인 Ray에 대응되는 모든 점들이 XX 가 될 수 있으므로 유일한 해를 결정할 수 없습니다. scale αα 값에 따라 빨간색 선의 어떤 점으로도 대응될 수 있기 때문에 1개의 카메라 만을 이용하면 무수히 많은 해를 가지게 됩니다.


Drawing


  • 반면 이상적인 상황에서 2개의 카메라에서 생성한 Ray의 교점을 이용하여 XX 의 좌표를 구하는 방법을 사용해 볼 수 있습니다. 이 경우에는 다음과 같은 식을 이용합니다.


  • {x=PXx=PX{x=PXx=PX(5)


  • 앞에서 가정한 바와 같이 x,xx,xmatched points는 노이즈를 포함할 수 밖에 없음을 가정하였습니다. 따라서 식(5)에서 완벽히 일치하는 해 XX 를 찾기는 어렵습니다. 대신에 가장 적합한 해를 찾는 최적화 방법을 이용해야 합니다.
  • 한 개의 점 X={X,Y,Z}X={X,Y,Z} 에 대하여 찾고자 하는 값은 3개이므로 최소 3개식을 이용해야 X,Y,ZX,Y,Z 를 추정할 수 있습니다.


  • 식을 구하기 위하여 2개의 Ray αxαxPXPX 가 서로 평행하다는 조건을 이용하도록 하겠습니다. 2개의 Ray가 평행하다면 cross product가 0이 되기 때문에 이 값을 이용하면 식을 구할 수 있습니다.


  • 먼저 아래와 같이 αx=PXαx=PX 를 전개해 보도록 하겠습니다.


  • x=PXx=PX(6)
  • α[uv1]=[p11p12p13p14p21p22p23p24p31p32p33p34][XYZ1]=[PT1PT2PT3][XYZ1]αuv1=p11p12p13p14p21p22p23p24p31p32p33p34⎢ ⎢ ⎢XYZ1⎥ ⎥ ⎥=⎢ ⎢PT1PT2PT3⎥ ⎥⎢ ⎢ ⎢XYZ1⎥ ⎥ ⎥(7)
  • {PT1=[p11p12p13p14]PT2=[p21p22p23p24]PT3=[p31p32p33p34]


  • 2개의 Ray αxPX 가 평행하기 때문에 다음 식을 만족합니다. scale factor α 는 아래 식에 영향이 없으므로 제거 할 수 있습니다.


  • x×PX=0


Drawing


  • aT=[a1a2a3]
  • bT=[b1b2b3]
  • a×b=[a2b3a3b2a3b1a1b3a1b2a2b1]
  • Cross product of two vectors in the same direction is zero: a×a=0


  • 식 (7)을 이용하여 전개하면 다음과 같습니다.


  • α[uv1]=[PT1PT2PT3][XYZ1]=[PT1XPT2XPT3X]
  • α[uv1]=[PT1XPT2XPT3X]


  • 앞에서 설명한 바와 같이 α 를 제거하겠습니다. 식 (14)에서 좌변의 αscale에 해당하므로 cross proudct 식을 적용할 때, normalized scale인 1 로 두어도 cross product 식 전개에는 영향이 없습니다. 따라서 다음과 같이 식을 전개할 수 있습니다.


  • [uv1]×[PT1XPT2XPT3X]=[vPT3XPT2XPT1XuPT3XuPT2XvPT1X]=[000]


  • 식 (15)에서 3행은 1행과 2행의 선형 결합으로 이루어져 있습니다. (1행에 x 를 곱하고 2행에 y 를 곱하여 더하면 3행이 도출됩니다.) 따라서 1행과 2행의 식만 이용할 수 있습니다. 즉, 식 (15)를 통해 2개의 2D point3D point 간의 대응 관계를 구할 수 있습니다.


  • [vPT3XPT2XPT1XuPT3X]=[00]
  • [vPT3PT2PT1uPT3]X=[00]


  • 같은 방식으로 x×PX=0 을 이용한 후 식(17)의 좌변 행렬에 누적시키면 다음과 같이 식을 구성할 수 있습니다.


  • [vPT3PT2PT1uPT3vPT3PT2PT1uPT3]X=[0000]


  • 식 (18)의 좌변을 간단히 다음과 같이 나타낼 수 있습니다.


  • AX=0



Python 실습


  • 아래 예제에서는 앞에서 배운 내용을 바탕으로 Triangulation을 적용하는 코드를 다루었습니다. Triangulated 3D pointX_true의 값이 동일한 것을 통해 개념에서 배운 내용이 유효한 것을 확인할 수 있습니다.


import numpy as np

def triangulate_point(x, P, x_prime, P_prime):
    u, v, _ = x
    u_prime, v_prime, _ = x_prime
    P_1 = P[0, :]
    P_2 = P[1, :]
    P_3 = P[2, :]

    P_prime_1 = P_prime[0, :]
    P_prime_2 = P_prime[1, :]
    P_prime_3 = P_prime[2, :]

    A = np.array([
        v * P_3 - P_2,
        P_1 - u * P_3,
        v_prime * P_prime_3 - P_prime_2,
        P_prime_1 - u_prime * P_prime_3
    ])

    # Perform SVD on matrix A
    _, _, Vt = np.linalg.svd(A)
    X = Vt[-1]

    # Convert from homogeneous to inhomogeneous coordinates
    X = X / X[-1]
    return X

# Define the true 3D point in homogeneous coordinates
X_true = np.array([45, -35, 150, 1])

# Define the camera matrices
P1 = np.array([[700, 120, 320, 80],
               [60, 650, 230, -50],
               [0.5, 0.3, 1, 0.1]])

P2 = np.array([[650, -100, 310, -140],
               [-80, 700, 240, 90],
               [0.4, -0.2, 1, 0.2]])

# Project the 3D point onto the image planes to get x and x'
x = P1 @ X_true
x = x / x[-1]  # Normalize to get inhomogeneous coordinates

x_prime = P2 @ X_true
x_prime = x_prime / x_prime[-1]  # Normalize to get inhomogeneous coordinates

# Triangulate the 3D point using the projections
X_triangulated = triangulate_point(x, P1, x_prime, P2)
print(f"Triangulated 3D point: {X_triangulated} Vs. X true: {X_true}.")
# Triangulated 3D point: [ 45. -35. 150.   1.] Vs. X true: [ 45 -35 150   1].

# Reproject the 3D point back to 2D to check
x_reprojected = P1 @ X_triangulated
x_reprojected = x_reprojected / x_reprojected[-1]

x_prime_reprojected = P2 @ X_triangulated
x_prime_reprojected = x_prime_reprojected / x_prime_reprojected[-1]

print(f"Original x: {x} Vs. Reprojected x: {x_reprojected}.")
print(f"Original x_prime: {x_prime} Vs. Reprojected x_prime: {x_prime_reprojected}.")
# Original x: [465.02159161  88.83405305   1.        ] Vs. Reprojected x: [465.02159161  88.83405305   1.        ].
# Original x_prime: [451.54109589  45.60502283   1.        ] Vs. Reprojected x_prime: [451.54109589  45.60502283   1.        ].