Generalized Hough Transform from scratch with python3
허프 변환(Hough Transformation)은 인라이어와 아웃라이어가 섞여있는 상황에서 인라이어 군집을 찾아내는 군집화 알고리즘이다. 허프변환의 이러한 특성을 활용하여 직선, 원 등과 같은 특정 도형을 찾아내는데 활용할 수 있는데, Generalized Hough Transformation(GHT)은 임의의 도형에 일반적으로 사용할 수 있도록 설계된 버전의 허프 변환이다.
허프 변환은 찾고자 하는 도형을 (도형의 특성을 반영한) 허프 공간으로 변환한 후, 누적 테이블(accumulator)을 만들어 투표(voting)하는 방식을 취한다. 이에 대한 자세한 설명을 위해, GHT를 사용하여 아래의 그림 (B)에서 (A)를 찾아내는 예시를 살펴보자. 먼저 (A)는 찾고자 하는 도형(우산)을 정의한 영상(reference image)이고, (B)는 (A)를 포함한 영상(target image)이다. GHT 수행을 위해선 먼저 찾고자 하는 도형인 (A)를 엣지의 방향(orientation)과 벡터로 구성된 허프공간으로 변화시켜야 하는데, GHT에서는 이를 특별히 'r-table'이라고도 한다.
r-table은 위에서 언급한 바와 같이 특정 도형에 대한 엣지의 방향과 벡터로 구성된 테이블인데, 달리 말하면 찾고자 하는 도형을, '해당 도형을 구성하고 있는 엣지의 방향과 벡터로 정의(표현)하는 테이블'이라고 생각하면 된다. 즉, 본 예시에서의 (A)를 엣지의 방향과 벡터로 표현하는 과정이 필요한데, 이는 (A)에 대한 (1) edge map, (2) oreintation map을 통해 표현할 수 있다. 먼저 edge map의 경우, canny 알고리즘을 사용하여 생성할 수 있고 그 결과는 다음과 같다. 이때 edge인 부분은 1(흰색), 아닌 부분은 0(검은색)으로 이진화한다.
엣지의 방향은 y, x 각각의 변화량을 사용하여 $atan(\frac{dy}{dx})$로 정의된다. 따라서 이를 위해 먼저 영상 (A)에 대하여 소벨(sobel) 등의 필터를 활용하여 $dy$, $dx$ 맵을 만들고 이를 사용하여 각 픽셀별 orientation을 생성할 수 있다. 단, r-table 생성을 위해선 각 orientation을 양자화하여 진행하는 것이 필요하다. 본 글에서는 8방향으로 양자화(0~7)하였다. (예 : -27.5도 ~ 27.5도 : 0, 27.5 ~ 72.5 : 1 등)
지금까지 (1) edge map, (2) orientation map 생성을 통해 r-table을 만들 준비를 완료했다. 다시 언급하자면, r-table은 orientation과 벡터로 해당 도형을 표현한 테이블이다. orientation의 경우, edge map에서 찾은 엣지의 위치 Ek(yk, xk)를 기준으로 orientation map에서도 동일한 위치인 Ok(yk, xk)에 위치한 orientation, $phi_k$를 가져오면 된다. 그렇다면 벡터는 무엇을 의미하는 것일까. 여기서 필요한 벡터는 위에서 찾은 엣지 Ek를 기준으로 정의되는 벡터를 의미한다. 벡터는 크기와 방향으로 정의되는 값인데, 그렇다면 해당 Ek를 기준으로 어떻게 방향과 크기를 정의해야하는 것인지 의문이 든다.
GHT에서는 해당 엣지 Ek를 기준으로 벡터를 정의하기 위해, '해당 도형의 임의의 중심'을 도입한다. 다시 말해 실제 도형의 중심이 아니라, 아무 점이나 임의로 중심으로 가정하고 해당 점을 기준으로 엣지에 대한 벡터를 정의한다. 예를 들어, 영상 (A)에서 임의로 중심 C(yc, xc)를 잡았다고 하면, 두 점 C와 Ek를 이용하여 1개의 벡터를 계산할 수 있다. 즉, 두 점 사이의 변량을 $\triangle y, \triangle x$라고 하면, Ek에서 C로의 방향 ok = $atan(\frac{\triangle y}{\triangle x})$가 되고, 크기 mk = $\sqrt(\triangle y^2 + \triangle x^2)$가 된다. 그럼 해당 지점 Ek의 벡터 rk는 (ok, mk)로 정의할 수 있다. 이렇게 정의된 rk는 앞에서 계산한 해당 엣지의 방향 $\phi_k$를 키(key)로 하여 테이블에 저장한다. 여기서 포인트는 방향 $\phi_k$를 '키로 하여' rk를 저장한다는 점이다. 이 말을 다시 풀어보자면, 만일 다른 지점에 위치한 엣지 Ej가 Ek와 동일한 엣지 방향($\phi_k$ )을 갖는다면, Ej를 기준으로 중심점 C와 계산한 벡터 rj도 같은 키($\phi_k$)를 기준으로 저장되어야 한다는 의미다.
해당 작업을 edge map을 기준으로 반복하게 되면, 서로 다른 위치에 있으나 같은 엣지 방향을 갖는 지점들에 대하여 벡터 r이 같은 군집으로 저장될 것이고, 이 과정을 통해 만들어진 r-table은 다음과 같은 형태가 될 것이다. 다시 말하지만, 이는 결국 주어진 (A)를 엣지에 따라 그 방향과 해당 방향에 따라 존재할 수 있는 다양한 벡터의 집합으로 표현한 형태가 된다. 일종의 해당 도형을 정의하는 '사양서'가 준비된 셈이다.
'사양서'가 준비되었으니, 해당 사양에 맞는 지점을 (B) target image에서 찾기만 하면 된다. 이를 위해 물체가 위치한 지점을 찾기 위한 누적 테이블을 이용하여 '투표(voting)'을 진행하면 된다. (B)에서 (A)의 위치를 찾는 것이 목표이므로, 투표를 수행할 누적 테이블은 영상 (B)와 동일한 크기로 설정한다. 이후 (A)와 마찬가지로 영상 (B)에 대하여 (1) edge map, (2) orienation map을 생성한다. 여기서도 orientation map은 양자화한다.
앞에서와 마찬가지로 영상 (B)의 edge map은 엣지의 위치를 찾기 위해 만든 영상으로, 중요한 것은 해당 엣지의 위치에 해당하는 엣지의 방향이다. 만일 (B)의 엣지 맵에서 Ea(ya, xa)가 주어질 경우, orientation map에서 Oa(ya, xa)에 위치한 엣지의 방향을 알 수 있다. 만일 Oa의 방향이 $\phi_a$라고 하면, r-table에서 해당 $\phi_a$에 속한 벡터값들이 존재할 것이다. 이를 rq, rw, re라고 하면, 이 세 벡터를 이용하여 (ya, xa)를 기준으로 3개의 새로운 지점을 다음과 같이 계산할 수 있다.
$y_{rq} = ya + m_q * sin(\alpha_q)$
$x_{rq} = xa + m_q * cos(\alpha_q)$
$y_{rw} = ya + m_w * sin(\alpha_w)$
$x_{rw} = xa + m_w * cos(\alpha_w)$
$y_{re} = ya + m_e * sin(\alpha_e)$
$x_{re} = xa + m_e * cos(\alpha_e)$
이렇게 산출한 세 지점, $(y_{rq}, x_{rq})$, $(y_{rw}, x_{rw})$, $(y_{re}, x_{re})$은 결과적으로 영상(B) 상에서 찾고자 하는 물체 (A)의 후보 위치가 되며, 누적 테이블에서 해당 후보 위치에 +1을 수행하면 된다.
(= voting) 이러한 작업을 영상 (B)의 모든 엣지에 대하여 반복 수행하고, 최종적으로 가장 많이 투표받은(+1이 된) 지점이 해당 물체의 위치(중심의 위치)가 된다.
영상(B)에서 추출된 엣지와 그 엣지의 방향을 추출했을 때, 해당 방향을 기준으로 r-table에서 존재하는 벡터들을 기준으로 후보 위치들을 계산한 후 누적하는 방식을 취했을 때, 영상(B)에서 찾고자 하는 물체의 위치가 가장 많은 표를 받는 이유를 직관적으로 생각해보면 아래의 그림과 같다. 위의 그림을 다시 가져와 보면, $\phi_k$를 키로 갖는 엣지들의 경우, 이에 해당하는 rj, rk를 기준으로 각각 두 개의 후보 위치를 계산하게 된다. 이때, Ek를 기준으로 rj로 계산된 후보 위치라든지, Ej를 기준으로 rk로 계산된 후보 위치의 경우, 중중복하여 투표를 받을 가능성은 낮다(Ea 기준 rb, Eb 기준 rc로 계산된 후보 위치도 마찬가지). 반면 실제 중심이었던 위치는 다른 지점 Ea, Eb로 부터도 중복하여 후보 위치로 투표받을 가능성이 높다는 것을 알 수 있다.
다만, 위 예시의 경우, 영상(A)와 영상(B)에 존재하는 우산의 크기와 회전 각도는 동일하다는 전제를 갖는다. 즉, 스케일괴 회전 변화는 고려하지 않은 GHT 방식으로, 스케일과 회전이 있는 경우 이 방식은 유효하지 않다. 이 경우 동일한 방식으로 GHT를 수행하되, y, x로 표현되는 위치와 함께 스케일, 회전을 함께 고려하는 4차원의 누적 테이블을 생성해야 한다. 다만, 이 경우 연산량이 많아 속도가 느릴 수 있다는 점에 유의해야 한다.
아래는 이러한 GHT를 파이썬으로 간단하게 구현해본 예시다.
(단 해당 코드에서는 canny, direction 함수 부분을 직접 구현하여 사용하였으나 해당 소스 부분은 제공하지 않는다. 따라서 간단히 테스트해볼 목적이라면 해당 부분만 cv2 등을 활용하거나 다시 직접 구현하면 된다.)
def generalizedHough(img1, img2):
"""
:param img1: 찾고자 하는 대상 영상(reference image)
ㄴ dtype : cv2.imgObj(grayscale)
:param img2: 찾고자 하는 물체를 갖고 있는 참조 이미지(target image)
ㄴ dtype : cv2.imgObj(grayscale)
단, 해당 모델은 스케일과 로테이션을 고려하지 않는다.
"""
# set an arbitrary center point of ref. image
# 설정한 점을 img1에 존재하는 물체의 중심이라고 가정한다.
cY, cX = img1.shape[0]/2, img1.shape[1]/2 # center of y, x
# STEP 1 : reference 영상 전처리
# get the orientation map of ref. image / get edge
refOri = direction(imgObj=img1, quantize=True) # orientation of ref. image
refEdge = canny(imgObj = img1, tHigh = 120, tLow = 60, sigma = 0.5)
# STEP 2 : Generate r-table(=orientation table)
r = {} # r-table
edges = np.where(refEdge==1)
for ey, ex in zip(edges[0], edges[1]):
# key index
p = refOri[ey, ex] # pi
# description
dY, dX = cY-ey, cX-ex
t = np.arctan2(dY, dX) # theta
d = np.sqrt(dY**2+dX**2) # distance
if p in r.keys():
r[p].append((t, d))
else:
r[p] = [(t, d)]
# STEP 3 : find location in img2 using "Accumulator"
acc = np.zeros_like(img2) # accumulator
targetEdge = canny(imgObj = img2, tHigh = 120, tLow = 60, sigma = 0.5)
targetOri = direction(imgObj = img2, quantize=True)
tEdges = np.where(targetEdge == 1) # target edges
for ey, ex in zip(tEdges[0], tEdges[1]):
tP = targetOri[ey, ex]
if tP in r.keys():
for t, d in r[tP]:
canY = ey + d*np.sin(t)
canX = ex + d*np.cos(t)
if canY < acc.shape[0] and canX < acc.shape[1]:
acc[int(canY), int(canX)] += 1
location = np.where(acc == np.max(acc))
plt.imshow(acc, cmap="gray")
plt.scatter(location[1], location[0], marker="+", c="red")
plt.show()
if __name__ == "__main__":
img1Dir = "./../data/umbrella.jpg"
img2Dir = "./../data/test_image4.jpg"
img1 = cv2.imread(img1Dir, cv2.IMREAD_GRAYSCALE)
img2 = cv2.imread(img2Dir, cv2.IMREAD_GRAYSCALE)
generalizedHough(img1=img1, img2=img2)