[MIT 머신비전 19강] 회전·이동·배율을 한 번에: 절대 방향의 닫힌 해와 RANSAC | 레일리 몫·이상치·회전 뽑기

MIT OpenCourseWare 6.801 19강 원본 영상 (영어)

한 줄 요약: 두 점 구름 사이의 이동은 무게중심 차이, 배율은 두 구름 크기의 비, 회전은 4×4 행렬의 최대 고유벡터로 모두 닫힌 해가 나옵니다. 다만 최소제곱은 튀는 점 몇 개에 크게 휘둘리므로, 실전에서는 RANSAC으로 정상 점만 골라 낸 뒤 풉니다.

사원수 ↔ 회전 행렬

지난 글에서 회전을 단위 사원수로 나타냈습니다. 실제 계산에서는 사원수와 회전 행렬을 자주 오가야 하므로 바꾸는 법부터 봅니다.

사원수 → 회전 행렬. \(\mathring{q}\,\mathring{r}\,\mathring{q}^*\)를 행렬 곱으로 풀어 쓰면 3×3 회전 행렬이 나옵니다.

\[R = \begin{pmatrix} q_0^2 + q_x^2 - q_y^2 - q_z^2 & 2(q_xq_y - q_0q_z) & 2(q_xq_z + q_0q_y) \\ 2(q_yq_x + q_0q_z) & q_0^2 - q_x^2 + q_y^2 - q_z^2 & 2(q_yq_z - q_0q_x) \\ 2(q_zq_x - q_0q_y) & 2(q_zq_y + q_0q_x) & q_0^2 - q_x^2 - q_y^2 + q_z^2 \end{pmatrix}\]

회전 각만 알고 싶다면 대각 성분의 합(trace)을 보면 됩니다. \(\text{tr}(R) = 3q_0^2 - (q_x^2 + q_y^2 + q_z^2) = 2\cos\theta + 1\)이므로

\[\cos\theta = \tfrac{1}{2}\big(\text{tr}(R) - 1\big)\]

그런데 이 식은 \(\theta\)가 0이나 180° 근처에서 부정확합니다. 코사인 그래프가 그 근처에서 평평해서, 값이 조금만 흔들려도 각도가 크게 바뀌기 때문입니다. 그럴 때는 \(\sin\)에 비례하는 비대각 성분을 쓰는 편이 낫습니다. 한쪽이 둔할 때 다른 쪽이 예민한 짝입니다.

회전 행렬 → 사원수. 대각 성분을 더하고 빼면 네 성분의 제곱이 각각 나옵니다.

\[\begin{aligned} 1 + r_{11} + r_{22} + r_{33} &= 4q_0^2 \\ 1 + r_{11} - r_{22} - r_{33} &= 4q_x^2 \\ 1 - r_{11} + r_{22} - r_{33} &= 4q_y^2 \\ 1 - r_{11} - r_{22} + r_{33} &= 4q_z^2 \end{aligned}\]

제곱근을 그냥 다 취하면 부호 조합이 너무 많아집니다. 그래서 넷 중 가장 큰 값 하나만 골라 그 성분을 구하고(숫자 정확도가 가장 좋음), 나머지는 비대각 성분의 합·차로 구합니다.

\[r_{32} - r_{23} = 4q_0q_x,\quad r_{13} - r_{31} = 4q_0q_y,\quad r_{21} - r_{12} = 4q_0q_z,\quad r_{21} + r_{12} = 4q_xq_y,\ \ldots\]

손으로 계산해 보면: \(z\)축 둘레 90° 회전 행렬은 \(R = \begin{pmatrix} 0 & -1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{pmatrix}\)입니다.

1. 대각 조합: \(4q_0^2 = 1 + 0 + 0 + 1 = 2\), \(4q_x^2 = 1 + 0 - 0 - 1 = 0\), \(4q_y^2 = 0\), \(4q_z^2 = 1 - 0 - 0 + 1 = 2\)
2. 가장 큰 쪽(공동 1등 \(q_0\))을 양수로 골라 \(q_0 = \sqrt{0.5} \approx 0.707\)
3. \(r_{21} - r_{12} = 1 - (-1) = 2 = 4q_0q_z\)이므로 \(q_z = 2/(4 \times 0.707) \approx 0.707\)

결과 \(\mathring{q} = (0.707, 0, 0, 0.707)\)은 지난 글에서 쓴 “z축 90°” 사원수와 정확히 같습니다.

배율까지 넣기

위성 사진으로 지형을 복원할 때처럼, 두 좌표계의 단위(배율)가 다를 수도 있습니다. 무게중심을 뺀 좌표로 쓰면

\[\mathbf{r}'_r = s\,R(\mathbf{r}'_l)\]

\(s > 0\)이 배율입니다. 계산에 쓸 세 값을 정의해 둡니다.

\[s_r = \sum_i \|\mathbf{r}'_{r,i}\|^2, \qquad s_l = \sum_i \|\mathbf{r}'_{l,i}\|^2, \qquad D = \sum_i \mathbf{r}'_{r,i} \cdot R(\mathbf{r}'_{l,i})\]

\(s_r\), \(s_l\)은 각 점 구름이 얼마나 퍼져 있는지(크기)를, \(D\)는 돌린 왼쪽과 오른쪽이 얼마나 같은 방향인지를 나타냅니다.

한쪽 방향으로 풀면 오차 \(\sum\|\mathbf{r}'_r - sR(\mathbf{r}'_l)\|^2 = s_r - 2sD + s^2 s_l\)을 \(s\)로 미분해 \(s = D/s_l\)을 얻습니다. 그런데 문제가 있습니다. 같은 문제를 오른쪽에서 왼쪽으로 풀어 얻은 배율의 역수가 이 값과 다릅니다.

한쪽 방향 최소제곱은 배율을 줄이면 점들이 서로 가까워져 오차가 작아 보이는 점을 “이용”합니다. 그래서 어느 쪽에서 푸느냐에 따라 답이 갈립니다.

대칭 버전은 오차를 양쪽에 나눠 줍니다.

\[e_i = \frac{1}{\sqrt{s}}\,\mathbf{r}'_{r,i} - \sqrt{s}\,R(\mathbf{r}'_{l,i}), \qquad \sum\|e_i\|^2 = \frac{1}{s}s_r - 2D + s\,s_l\]

미분해서 0으로 놓으면

\[-\frac{s_r}{s^2} + s_l = 0 \quad\Longrightarrow\quad s = \sqrt{\frac{s_r}{s_l}}\]

방향을 바꾸면 \(s_r\)과 \(s_l\)이 자리만 바꾸므로 정확히 역수가 됩니다. 게다가 이 식에는 \(D\)도, \(R\)도 없습니다. 배율은 그냥 두 점 구름의 크기 비율이고, 대응도 회전도 몰라도 구할 수 있습니다.

정리하면 이동, 배율, 회전이 서로 분리됩니다. 배율과 이동은 간단한 공식, 남는 것은 회전뿐입니다.

회전: 레일리 몫으로 보기

회전은 지난 강의처럼 “길이 1인 \(\mathring{q}\) 중에서 \(\mathring{q}^T N \mathring{q}\)를 최대로”입니다. 이번에는 길이 조건을 나눗셈으로 처리하는 방법을 봅니다. 목적 함수를 \(\mathring{q}\)의 길이 제곱으로 나누면

\[J(\mathring{q}) = \frac{\mathring{q}^T N \mathring{q}}{\mathring{q}^T \mathring{q}}\]

이 값은 \(\mathring{q}\)를 몇 배로 늘리든 변하지 않습니다. 그래서 조건 없이 마음대로 미분할 수 있습니다. 이런 꼴을 레일리 몫(Rayleigh quotient)이라고 합니다. 미분해서 0으로 놓으면

\[N\mathring{q} = \underbrace{\frac{\mathring{q}^T N \mathring{q}}{\mathring{q}^T \mathring{q}}}_{\lambda}\,\mathring{q}\]

\(N\)을 곱해도 방향이 그대로인 벡터, 즉 고유벡터이고, 목적 함수의 값이 곧 고유값 \(\lambda\)입니다. 최대화 문제이니 가장 큰 고유값의 고유벡터를 고르면 됩니다. 라그랑주 승수로 푼 지난 강의와 같은 답입니다.

회전 행렬로 직접 풀려면 복잡한 제약 최적화나 특이값 분해가 필요한데, 사원수로는 4×4 고유값 문제 하나로 끝납니다.

짝은 몇 개나 필요할까

미지수는 이동 3 + 회전 3 = 6개(배율까지 7개)입니다. 대응 점 한 쌍은 3D 좌표가 같아야 하니 조건 3개를 줍니다. 그럼 2쌍이면 충분할까요?

두 점 사이의 거리는 원래 정해져 있으니, 두 번째 쌍은 사실상 조건 2개만 더합니다.

• 2쌍: 두 점을 잇는 직선 둘레로 돌려도 두 쌍이 다 맞습니다. 조건이 5개뿐입니다.
• 3쌍: 한 줄에 있지만 않으면 자세가 하나로 정해집니다. 조건 9개로 오히려 남습니다.
• 일반 선형 변환(회전 대신 아무 3×3 행렬 + 이동)까지 허용하면 미지수 12개라 4쌍이 필요합니다. 하지만 라이다 두 대의 관계처럼 실제 문제는 대부분 회전 + 이동이면 충분합니다.

언제 실패할까: 점들이 한 평면에 있을 때

대응이 부족하면 \(N\)이 특이 행렬이 되어 고유값 0이 생깁니다. 더 흥미로운 경우는 점들이 한 평면 위에 있을 때입니다.

\(N\)의 고유값은 4차 방정식 \(\lambda^4 + c_3\lambda^3 + c_2\lambda^2 + c_1\lambda + c_0 = 0\)의 근이고, 지난 강의에서 본 대로 \(c_3 = 0\), \(c_1 = -8\det M\)입니다(\(M = \sum_i \mathbf{r}'_{l,i}\mathbf{r}'^T_{r,i}\)). 점들이 법선 \(\hat{n}\)인 평면 위에 있으면 모든 점이 \(\mathbf{r}' \cdot \hat{n} = 0\)이라

\[M\hat{n} = \sum_i \mathbf{r}'_{l,i}\,(\mathbf{r}'_{r,i}\cdot\hat{n}) = \mathbf{0}\]

\(M\)이 특이 행렬이 되고 \(c_1 = 0\)입니다. 방정식이 \(\lambda^4 + c_2\lambda^2 + c_0 = 0\)으로 줄어, \(\lambda^2\)에 대한 2차 방정식이 됩니다.

양쪽 점 구름이 모두 평면이면 오히려 더 쉽게 풀 수 있습니다.

1. 평면끼리 포개기: 두 평면의 법선 \(\hat{n}_1\), \(\hat{n}_2\)로 축과 각도를 바로 얻습니다.

   \[\hat{\omega} = \frac{\hat{n}_1 \times \hat{n}_2}{\|\hat{n}_1 \times \hat{n}_2\|}, \qquad \theta = \text{atan2}\big(\|\hat{n}_1 \times \hat{n}_2\|,\; \hat{n}_1\cdot\hat{n}_2\big)\]

2. 평면 안에서 돌리기: 포갠 뒤에는 2D 회전 하나만 남으니 2D 최소제곱으로 풉니다.

튀는 점 하나가 다 망친다: RANSAC

최소제곱은 잡음이 정규분포(종 모양)일 때 가장 좋은 추정입니다. 하지만 잘못 짝지어진 점(이상치) 몇 개만 섞여도 크게 휘둘립니다. 오차를 제곱하니까 멀리 튄 점의 목소리가 지나치게 커지기 때문입니다. 제곱 대신 절댓값을 쓰면 덜 휘둘리지만 닫힌 해가 없습니다.

대안이 RANSAC(RANdom SAmple Consensus, 무작위 표본 합의)입니다.

RANSAC은 “가장 많은 점이 동의하는 답”을 찾은 뒤, 동의한 점들로만 다시 최소제곱을 합니다.

1. 무작위로 최소 개수만 뽑습니다. 직선이면 2개, 절대 방향이면 3쌍입니다.
2. 그 점들로 답을 만듭니다.
3. 나머지 점으로 검사합니다. 답에서 \(\epsilon\) 이내(폭 \(2\epsilon\)인 띠 안)에 드는 점을 정상치로 셉니다.
4. 여러 번 반복해 정상치가 가장 많은 답을 고르고, 그 정상치들로만 최소제곱을 다시 합니다.

\(\epsilon\)은 보통 측정 잡음의 크기를 알고 정합니다. 강의에서는 RANSAC을 14강의 허프 변환으로도 해석합니다. 뽑을 때마다 나온 답을 매개변수 공간의 칸에 한 표씩 던지고, 표가 가장 많이 쌓인 칸을 고르는 것과 같다는 겁니다.

회전을 고르게 뽑기

닫힌 해가 없는 방향 문제에서는 회전을 이것저것 시험해 봐야 할 때가 있습니다. 그러려면 회전의 공간, 즉 4D 단위구에서 점을 고르게 뽑아야 합니다. 3D 구로 먼저 연습해 봅시다.

위도·경도 격자는 극 근처에서 경선이 모여 촘촘해집니다. 정육면체는 모서리 쪽이 구에서 더 멀어 그쪽 방향 점이 많아집니다.

• 위도·경도를 고르게: 극에 몰립니다.
• 정육면체 안에서 고르게 뽑아 구에 투영: 정육면체 모서리 방향에 몰립니다.
• 구 밖에 떨어진 점은 버리고 투영: 고르게 퍼집니다. 원점에 너무 가까운 점도 숫자 문제로 버립니다. 4D에서도 똑같이 하면 단위 사원수를 고르게 뽑을 수 있습니다.

버리는 점이 아깝다면 정다면체(정사면체 4면, 정육면체 6면, 정팔면체 8면, 정십이면체 12면, 정이십면체 20면)처럼 구에 더 꼭 맞는 도형을 씁니다. 회전 자체를 규칙적으로 뽑으려면 정육면체의 회전 대칭을 이용합니다. 기본 회전 10개(제자리, 세 축 둘레 180°, 세 축 둘레 ±90°)를 사원수로 쓰고 서로 곱해 새 회전을 만듭니다.

\[\tfrac{1}{\sqrt{2}}(1, \hat{x}) \cdot \tfrac{1}{\sqrt{2}}(1, \hat{y}) = \tfrac{1}{2}\big(1 - \hat{x}\cdot\hat{y},\; \hat{x} + \hat{y} + \hat{x}\times\hat{y}\big) = \tfrac{1}{2}(1,\; \hat{x} + \hat{y} + \hat{z})\]

\(\cos(\theta/2) = 1/2\)이니 \(\theta = 120°\), 축은 \((1, 1, 1)/\sqrt{3}\)입니다. 정육면체의 꼭짓점을 관통하는 축 둘레 120° 회전으로, 기본 10개에 없던 새 회전입니다. 이렇게 곱해 나가면 정육면체 회전군의 24개 회전이 모두 나옵니다.

직접 해 보기

① 절대 방향 닫힌 해. 파란 점 구름(왼쪽 좌표계)과 주황 점 구름(오른쪽 좌표계)이 있고, 초록 고리는 4×4 행렬 \(N\)의 최대 고유벡터로 구한 회전과 이동으로 파란 점을 옮긴 위치입니다. 이상치 개수를 늘려 보세요.

이상치가 없을 때 회전 오차는 0.16°입니다. 이상치를 5개(16개 중) 넣으면 25.82°로 망가지고, 초록 고리가 주황 점에서 멀리 벗어납니다. “RANSAC 쓰기”를 켜면 정상치 11개를 골라 오차가 0.42°로 돌아옵니다. 절반인 8개가 이상치여도 RANSAC은 버팁니다.

② 직선 RANSAC. “한 번 뽑기”를 몇 번 눌러 보라색 두 점이 세운 직선과 그 직선의 정상치 수를 보세요. “50번 뽑기”를 누르면 가장 많은 표를 받은 직선이 남습니다.

처음 데이터에서 50번 뽑으면 최소제곱의 기울기 오차가 8.3°인 반면, RANSAC은 0.1°입니다. 이상치 비율을 0.7까지 올려도 대개 버티지만, 뽑아야 하는 횟수가 늘어납니다. 두 점이 모두 정상치일 확률이 \(0.3^2 = 0.09\)로 떨어지기 때문입니다.

정리

• 사원수 ↔ 회전 행렬 변환은 가장 큰 대각 조합을 골라 숫자 정확도를 지킵니다. 각도는 trace보다 비대각 성분이 안정적일 때가 있습니다.
• 배율은 대칭 버전 \(s = \sqrt{s_r/s_l}\)을 써야 방향을 바꿔도 일관됩니다. 대응도 회전도 필요 없습니다.
• 회전은 레일리 몫 \(\mathring{q}^TN\mathring{q}/\mathring{q}^T\mathring{q}\)의 최대화, 즉 \(N\)의 최대 고유벡터입니다.
• 대응은 3쌍 이상 필요하고, 점들이 한 평면에 있으면 \(\det M = 0\)이 되어 특별히 다뤄야 합니다.
• 이상치에는 RANSAC: 최소 개수를 뽑아 답을 세우고, 가장 많은 점이 동의하는 답을 고릅니다.

참고: MIT 6.801/6.866 Machine Vision (Fall 2020), Berthold Horn · Ryan Sander · Tadayuki Yoshitake 강의 노트 (MIT OpenCourseWare, CC BY-NC-SA 4.0)