필립 클라인의 저서 코딩 더 매트릭스 6장 기저

교환 보조정리를 살펴봅니다.
기저 변경과 원근감 렌더링을 활용해 이미지의 원근감을 수정해 봅니다.
일차종속과 일차독립, 기저, 교환 정리와 관련된 프로시저를 작성해 봅니다.

교환 보조정리

$S$ 는 벡터들의 집합이라 하고 $A$ 는 $S$ 의 부분집합이라 할 때 $z$ 는 $Span\text{ }S$ 의 벡터이고 $A \cup \{z\}$ 는 일차독립이라 하겠습니다. 그러면 $Span\text{ }S = Span\text{ }(\{z\} \cup S - \{w\})$ 를 만족하는 $w \in S - A$ 가 존재합니다. $Span\text{ }(\{z\} \cup S - \{w\})$ 가 $Span\text{ }S$ 내의 $w$ 를 제외하고 $z$ 를 포함시키는 모양새라 교환정리라 부릅니다.

교환정리를 증명해 보겠습니다. $S = \{v_1, ..., v_n, w_1, ..., w_m\}$ , $A = \{v_1, ..., v_n\}$ 라 할때 $Span\text{ }S$ 내에 있는 $z$ 는 다음과 같이 쓸 수 있습니다.

z = \alpha_1 v_1 + ... + \alpha_n v_n + \beta_1 w_1 + ... + \beta_m w_m

여기서 $A \cup \{z\}$ 는 일차독립이므로 $\beta_1, ..., \beta_m$ 중 적어도 하나는 0이 아닙니다. 0이 아닌 계수 $\beta_j$ 가 있을때 아래와 같이 식을 다시 써보겠습니다.

\begin{aligned} w_j = & \frac{1}{\beta_j} z + (-\frac{\alpha_1}{\beta_j} v_1) + ... + (-\frac{\alpha_n}{\beta_j} v_n) + (-\frac{\beta_1}{\beta_j}) w_1 + ... + (-\frac{\beta_{j-1}}{\beta_j}) w_{j-1} \\ & + (-\frac{\beta_{j+1}}{\beta_j}) w_{j+1} + ... + (-\frac{\beta_m}{\beta_j}) w_m \end{aligned}

벡터 $w_j \in S$ 는 다른 벡터들의 선형결합으로 표현 가능하므로 Superfluous-Vector 정리에 의해 $Span\text{ }(\{z\} \cup S - \{w\}) = Span\text{ }(\{z\} \cup S) = Span\text{ }S$ 입니다.

원근감 수정

원근감 렌더링에서는 카메라 좌표상 점들을 이미지 센서 어레이에 픽셀로 표현했었습니다. 이번에는 원근감을 가진 이미지에 새로운 기저를 적용해 원근감을 수정해보겠습니다. 원근감 수정의 직관적 이해를 위해 결과 이미지를 먼저 보겠습니다.

원본

변경된 기저로 표현한 이미지

원본이미지는 카메라 좌표 기저에 대한 표현에서 건물의 한 면의 좌표 기저에 대한 표현으로 변환됐습니다.

원근감 수정은 다음 두 단계로 완성됩니다.

이미지상 픽셀들을 카메라 좌표 표현에서 건물 한 면에 대한 좌표 표현으로 변환하는 함수 $h$
건물 한 면에 대한 좌표 표현을 이미지 센서 어레이 상에 매핑하는 함수 $g$

따라서 구하고자 하는 최종 결과는 함수 $f = g \circ h$ 에 의해 완성됩니다.

첫번째 단계는 기저변경을 뜻합니다. 원본 그림에 건물 한 면의 기저를 표시해 놓겠습니다.

원본의 카메라 좌표계를 $a_1,a_2,a_3$ 라 하고 이미지 상의 한 점 $p$ 의 카메라 좌표계에 대한 좌표 표현을 $[x_1,x_2,x_3]$ 라 하겠습니다. 건물 한 면의 좌표계는 $c_1,c_2,c_3$ 라 하고 점 $p$ 의 건물 한 면의 좌표계에 대한 좌표 표현은 $[y_1,y_2,y_3]$ 라 하겠습니다. 같은 점 $p$ 를 서로 다른 좌표계로 표현하므로 다음처럼 쓸 수 있습니다.

\begin{bmatrix} a_1 & a_2 & a_3 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} c_1 & c_2 & c_3 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}

$A = \begin{bmatrix}a_1 & a_2 & a_3\end{bmatrix}$ , $C = \begin{bmatrix}c_1 & c_2 & c_3\end{bmatrix}$ 라 하면 함수 $y \to Cy$ 는 $\Bbb{R}^3$ 에서 $\Bbb{R}^3$ 로의 가역함수이므로 $C^{-1}$ 이 존재합니다. 따라서 건물 한 면에 대한 좌표 표현 $[y_1,y_2,y_3]$ 는 식 $y = C^{-1}Ax$ 로 구할 수 있습니다. 여기서 $C^{-1}A$ 를 $H$ 라 하면 아래와 같이 쓸 수 있습니다.

H \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix}

이제 3 X 3 행렬 $H$ 의 엔트리들을 구해야 합니다.

\begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} h_{11} & h_{12} & h_{13} \\ h_{21} & h_{22} & h_{23} \\ h_{31} & h_{32} & h_{33} \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}

위 행렬방정식을 정리하면 다음을 얻습니다.

\begin{aligned} y_1 = h_{1i} \cdot x \\ y_2 = h_{2i} \cdot x \\ y_3 = h_{3i} \cdot x \end{aligned}

한편 건물의 한 면에 대한 좌표 표현 $[y_1,y_2,y_3]$ 를 이미지 상에 매핑하는 함수 $g$ 는 $g([y_1,y_2,y_3]) = [w_1,w_2,1]$ 과 같이 쓸 수 있습니다. 따라서 $w_1 = \frac{y_1}{y_3}$ , $w_2 = \frac{y_2}{y_3}$ 입니다. 식의 양변에 $y_3$ 을 곱하면 다음을 얻습니다.

\begin{aligned} w_1 y_3 = y_1 \\ w_2 y_3 = y_2 \end{aligned}

이 식에 앞서 정리한 $y_1 = h_{1i} \cdot x$ , $y_2 = h_{2i} \cdot x$ , $y_3 = h_{3i} \cdot x$ 를 대입하면 다음을 얻습니다.

\begin{aligned} w_1(h_{3i} \cdot x) - h_{1i} \cdot x = 0 \\ w_2(h_{3i} \cdot x) - h_{2i} \cdot x = 0 \end{aligned}

$x = [x_1,x_2,x_3]$ 는 카메라 좌표에 의한 좌표 표현인 동시에 원본 이미지 픽셀의 좌표 표현이고, $[w_1,w_2,1]$ 은 새로운 기저에 대한 좌표 표현이니 다음과 같이 이미지 상 픽셀 $x$ 를 각각 새로운 기저에 대한 좌표 표현으로 매핑할 수 있습니다.

이제 이미지 위의 4개에 점에 대해 각각을 $[x_1,x_2,1]$ 와 $[w_1,w_2,1]$ 에 대입해 식 2개씩 8개의 식을 얻을 수 있습니다. 이를 행렬로 표현하면 다음과 같습니다.

\begin{bmatrix} -x_{a1} & -x_{a2} & -x_{a3} & 0 & 0 & 0 & w_1 x_{a1} & w_1 x_{a2} & w_1 x_{a3} \\ \\ 0 & 0 & 0 & -x_{a1} & -x_{a2} & -x_{a3} & w_2 x_{a1} & w_2 x_{a2} & w_1 x_{a3} \\ \\ -x_{b1} & -x_{b2} & -x_{b3} & 0 & 0 & 0 & w_1 x_{b1} & w_1 x_{b2} & w_1 x_{b3} \\ \\ 0 & 0 & 0 & -x_{b1} & -x_{b2} & -x_{b3} & w_2 x_{b1} & w_2 x_{b2} & w_1 x_{b3} \\ \\ -x_{c1} & -x_{c2} & -x_{c3} & 0 & 0 & 0 & w_1 x_{c1} & w_1 x_{c2} & w_1 x_{c3} \\ \\ 0 & 0 & 0 & -x_{c1} & -x_{c2} & -x_{c3} & w_2 x_{c1} & w_2 x_{c2} & w_1 x_{c3} \\ \\ -x_{d1} & -x_{d2} & -x_{d3} & 0 & 0 & 0 & w_1 x_{d1} & w_1 x_{d2} & w_1 x_{d3} \\ \\ 0 & 0 & 0 & -x_{d1} & -x_{d2} & -x_{d3} & w_2 x_{d1} & w_2 x_{d2} & w_1 x_{d3} \\ \\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} h_{11} \\ h_{12} \\ h_{13} \\ h_{21} \\ h_{22} \\ h_{23} \\ h_{31} \\ h_{32} \\ h_{33} \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{bmatrix}

행렬의 마지막 행에 $[1,0,...,0]$ 이 추가되었습니다. $\alpha Hx = [\alpha y_1, \alpha y_2, \alpha y_3]$ 일때 $g(\alpha Hx) = [\frac{\alpha y_1}{\alpha y_3}, \frac{\alpha y_2}{\alpha y_3}, \frac{\alpha y_3}{\alpha y_3}]$ 이므로 이는 이미지로부터 건물 한 면의 크기를 정할 수 없음을 의미합니다. 따라서 행렬 $H$ 의 엔트리 $h_{11}$ 를 1로 만드는 스케일링 행 $[1,0,...,0]$ 을 추가해 $H$ 의 엔트리들을 얻습니다. 위 과정을 numpy를 이용한 코드로 옮겨보겠습니다. 이미지를 행렬로 바꾸고 브라우저에 그리기 위해서 소스파일에서 image_mat_util을 다운로드 받습니다.

from image_mat_util import file2mat, mat2display
from vec import Vec
import numpy as np

(X_pts, colors) = file2mat('../assets/perspective.png')

# 좌표쌍으로부터 2개의 식을 추출
def np_make_equation(x1,x2,w1,w2):
    x3 = 1
    u = [ -x1, -x2, -x3, 0, 0, 0, w1 * x1, w1 * x2, w1 * x3  ]
    v = [ 0, 0, 0, -x1, -x2, -x3, w2 * x1, w2 * x2, w2 * x3  ]
    return (u, v)

# 스케일링 행
def np_scale_equation():
    return [1, 0, 0, 0, 0, 0, 0, 0, 0]

left_top = np_make_equation(127, 66, 0, 0)
right_top = np_make_equation(160, 71, 1, 0)
left_bottom = np_make_equation(126, 107, 0, 1)
right_bottom = np_make_equation(160, 102, 1, 1)

L = np.matrix([
    left_top[0], left_top[1],
	right_top[0], right_top[1],
    left_bottom[0], left_bottom[1],
    right_bottom[0], right_bottom[1],
    np_scale_equation()
])

b = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1])

H_np = np.linalg.solve(L, b)
H_np = H_np.reshape(3,3)

변수 H_np는 numpy의 array이므로 이를 Mat 인스턴스로 변환하겠습니다.

from mat import Mat
from matutil import zero_mat

def np2mat(H):
    R = ['y1', 'y2', 'y3']
    C =  [ 'x', 'y', 'u' ]
    M = zero_mat((set(R), set(C)))
    for i, r in zip(range(len(H)), R):
        for j, c in zip(range(len(H[0])), C):
            M[r,c] = H[i][j]
    return M

H = np2mat(H_np)
print(H)
Y = H * X_pts

>>>
            u      x      y
      ---------------------
 y1  |   -129      1 0.0244
 y2  |  -51.1 -0.166   1.09
 y3  |     83 -0.323 0.0244

두 단계 중 첫번째 단계, 기저변경을 완성했습니다. 이제 건물 한 면에 대한 좌표 표현들을 모두 이미지 평면 상에 매핑하는 프로시저 mat_move2board을 작성하겠습니다.

from matutil import mat2coldict, coldict2mat
from vec import Vec

def mat_move2board(M):
    m = {}
    coldict = mat2coldict(M)
    for k, item in coldict.items():
        m[k] = Vec(item.D, { d: item[d] / item['y3'] for d in item.D })
    return coldict2mat(m)

Y_board = mat_move2board(Y)
mat2display(Y_board, colors, ('y1', 'y2', 'y3'), scale=100, xmin=None, ymin=None)

기저 관련 프로시저

벡터리스트의 한 벡터가 일차종속인지 판단하는 프로시저 is_superfluous

from solver import solve
from matutil import coldict2mat
from vecutil import zero_vec

def is_superfluous(L, i):
    if len(L) == 1:
        return False
    M = coldict2mat(L[0:i] + L[i + 1:])
    u = solve(M, L[i])
    return (L[i] - M * u).is_almost_zero()

a0 = Vec({ 'a', 'b', 'c', 'd' }, { 'a': 1 })
a1 = Vec({ 'a', 'b', 'c', 'd' }, { 'b': 1 })
a2 = Vec({ 'a', 'b', 'c', 'd' }, { 'c': 1 })
a3 = Vec({ 'a', 'b', 'c', 'd' }, { 'a': 1, 'c': 1 })

# a1은 다른 벡터들 [a0, a2, a3]의 선형결합으로 표현할 수 없음
assert is_superfluous([a0, a1, a2, a3], 1) == False

# a3은 a0와 a2의 선형결합으로 표현할 수 있음
assert is_superfluous([a0, a1, a2, a3], 3) == True

# a0은 다른 벡터들 [a1, a3]의 선형결합으로 표현할 수 없음
assert is_superfluous([ a0, a1, a3 ], 0)  == False

소스파일에서 제공하는 solver.py의 프로시저 solve를 이용합니다. 행렬방정식 $Mu = v$ 가 해를 가지면 벡터 $v$ 는 $M$ 의 열벡터들의 선형결합으로 표현가능하므로 $v$ 는 $M$ 의 열벡터들에 대해 일차종속입니다. 해 $u$ 가 존재하면 L[i] - M * u의 값이 영벡터일 것입니다. 하지만 Python의 부동소수점 오차를 감안해 L[i] - M * u의 엔트리가 0이 아닌 0에 가까운 값일 수 있으므로 Vec의 인스턴스 메서드 is_almost_zero를 이용해 값을 검증합니다. is_almost_zero는 다음과 같습니다.

def is_almost_zero(self):
	s = 0
	for x in self.f.values():
		if isinstance(x, int) or isinstance(x, float):
			s += x*x
		elif isinstance(x, complex):
			y = abs(x)
			s += y*y
		else: return False
	return s < 1e-20

벡터리스트가 일차독립인지 판단하는 프로시저 is_independent

def is_independent(L):
    for i in range(len(L)):
        if is_superfluous(L, i):
            return False
    return True

assert is_independent([a0,a1,a2,a3]) == False
assert is_independent([a0,a1,a2]) == True
assert is_independent([a0,a2,a3]) == False
assert is_independent([a0,a1,a3]) == True

is_superfluous를 벡터들의 리스트 L의 반복에 서브루틴으로서 이용합니다.

주어진 벡터리스트들에 대한 기저를 찾는 프로시저 subset_basis

def subset_basis(T):
    subset = []
    for v in T:
        if not is_superfluous([v] + subset, 0):
            subset = [v] + subset
    return subset

print(subset_basis([a0,a3,a1,a2]))

>>>
[
	Vec({'d', 'c', 'b', 'a'},{'b': 1}),
	Vec({'d', 'c', 'b', 'a'},{'a': 1, 'c': 1}),
	Vec({'d', 'c', 'b', 'a'},{'a': 1})
]

전체 프로세스는 Grow 알고리즘을 따릅니다. is_superfluous프로시저를 이용해 기저들이 추가되고 있는 리스트 subset와 T의 한 벡터 v가 일차종속인지 판단후 일차독립이라면 subset에 추가합니다.

교환 보조정리의 $w$ 를 구하는 프로시저 exchange

def exchange(S, A, z):
    for i in range(len(S)):

		# S[i] = w 가 A의 생성에 속하면 더 이상 진행하지 않습니다.
        if not is_independent(A + [S[i]]):
            continue

		# M = Span (S + {z} - {w})
        M = S[0:i] + [z] + S[i+1:]

		# S의 모든 벡터가 M에 속하는지 검증합니다.
        for v in S:
            if not is_superfluous([v] + M, 0):
                break
        else:
            return S[i]

S = [ list2vec(v) for v in [ [ 0,0,5,3 ], [2,0,1,3], [0,0,1,0], [1,2,3,4] ] ]
A = [ list2vec(v) for v in [ [0,0,5,3], [2,0,1,3] ] ]
z = list2vec([ 0,2,1,1 ])
print(exchange(S, A, z))

>>>
 0 1 2 3
--------
 0 0 1 0

변수 $z = [0,2,1,1]$ 는 $S$ 의 선형결합으로 표현가능하고 $A \cup \{z\}$ 는 일차독립입니다.

-\frac{1}{2} \begin{bmatrix} 0 \\ 0 \\ 5 \\ 3 \end{bmatrix} -\frac{1}{2} \begin{bmatrix} 2 \\ 0 \\ 1 \\ 3 \end{bmatrix} + \begin{bmatrix} 0 \\ 0 \\ 1 \\ 0 \end{bmatrix} + \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} = \begin{bmatrix} 0 \\ 2 \\ 1 \\ 1 \end{bmatrix}

프로시저는 오직 $S - A$ 에 속하는 S[i]에 대해서만 6번째 줄 이후 코드를 실행합니다. 변수 M은 $Span\text{ }(S \cup \{z\} - \{w\})$ 을 뜻합니다. M에 대해 S의 모든 벡터가 생성에 속하는지 판단합니다. 만약 M이 모든 벡터를 포함한다면 S[i]를 $w$ 로서 반환합니다.

코딩 더 매트릭스 - 6장 기저 (2)

교환 보조정리, 기저 변경을 이용한 원근감 수정

2018년 12월 17일

교환 보조정리

원근감 수정

기저 관련 프로시저