
Triangulation
2021, Feb 15
- 참조 : https://www.cs.cmu.edu/~16385/s17/Slides/
- 사전 지식 : Direct Linear Transformation
- 사전 지식 : 특이값 분해 (Singular Value Decomposition)
- 사전 지식 : 카메라 모델 및 카메라 캘리브레이션의 이해와 Python 실습
목차
Triangulation 개념

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

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,P′Camera projection matrices: P,P′(2)
- Estimated 3D points: XT=[XYZ]Estimated 3D points: XT=[XYZ](3)
- 식 (1)에서 x,x′x,x′ 각각은 이미지 상에서 서로 같은 지점의 이미지 좌표를 의미합니다.
- 식 (2)에서 PP 는
Camera Extrinsic Matrix
와Camera Intrinsic Matrix
를 모두 반영한Camera Projection Matrix
를 의미합니다. 상세 내용은 글 상단의카메라 캘리브레이션
링크를 참조해 보시기 바랍니다. - 식 (3)의 XX 는 식 (1)의 이미지 2D points인 x,x′x,x′ 가 이미지에 투영되기 이전의 3D points를 의미하며 최종 추정하고자 하는 값이 됩니다.
- x=PXx=PX(4)
- 식 (4)에서
3D points
인 XX 를Camera Projection Matrix
를 통해 이미지에 투영하면 xx 가 됩니다. - 이 때, XX 를 알고 싶으나 XX 의 후보가 되는 값은 무수히 많이 있습니다. 관련 내용은 아래 링크를 참조하시면 됩니다.
- 동차좌표계 : https://gaussian37.github.io/vision-concept-homogeneous_coordinate/

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

- 반면 이상적인 상황에서 2개의 카메라에서 생성한
Ray
의 교점을 이용하여 XX 의 좌표를 구하는 방법을 사용해 볼 수 있습니다. 이 경우에는 다음과 같은 식을 이용합니다.
- {x=PXx′=P′X{x=PXx′=P′X(5)
- 앞에서 가정한 바와 같이 x,x′x,x′ 의
matched points
는 노이즈를 포함할 수 밖에 없음을 가정하였습니다. 따라서 식(5)에서 완벽히 일치하는 해 XX 를 찾기는 어렵습니다. 대신에 가장 적합한 해를 찾는 최적화 방법을 이용해야 합니다. - 한 개의 점 X={X,Y,Z}X={X,Y,Z} 에 대하여 찾고자 하는 값은 3개이므로 최소 3개식을 이용해야 X,Y,ZX,Y,Z 를 추정할 수 있습니다.
- 식을 구하기 위하여 2개의
Ray
αxαx 와 PXPX 가 서로 평행하다는 조건을 이용하도록 하겠습니다. 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
αx 와 PX 가 평행하기 때문에 다음 식을 만족합니다.scale factor
α 는 아래 식에 영향이 없으므로 제거 할 수 있습니다.
- x×PX=0

- aT=[a1a2a3]
- bT=[b1b2b3]
- a×b=[a2b3−a3b2a3b1−a1b3a1b2−a2b1]
- 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]=[vPT3X−PT2XPT1X−uPT3XuPT2X−vPT1X]=[000]
- 식 (15)에서 3행은 1행과 2행의 선형 결합으로 이루어져 있습니다. (1행에 x 를 곱하고 2행에 y 를 곱하여 더하면 3행이 도출됩니다.) 따라서 1행과 2행의 식만 이용할 수 있습니다. 즉, 식 (15)를 통해 2개의
2D point
와3D point
간의 대응 관계를 구할 수 있습니다.
- [vPT3X−PT2XPT1X−uPT3X]=[00]
- [vPT3−PT2PT1−uPT3]X=[00]
- 같은 방식으로 x′×P′X=0 을 이용한 후 식(17)의 좌변 행렬에 누적시키면 다음과 같이 식을 구성할 수 있습니다.
- [vPT3−PT2PT1−uPT3v′P′T3−P′T2P′T1−u′P′T3]X=[0000]
- 식 (18)의 좌변을 간단히 다음과 같이 나타낼 수 있습니다.
- AX=0
- 식 (19)에서 X 이
trivial solution
이 아닌 해를 구하려면 특이값 분해 (Singular Value Decomposition)를 이용하여 구할 수 있습니다. 상세 내용은 아래 링크에서 확인할 수 있습니다.- 링크 : SVD 활용
Python 실습
- 아래 예제에서는 앞에서 배운 내용을 바탕으로
Triangulation
을 적용하는 코드를 다루었습니다.Triangulated 3D point
와X_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. ].