Feature Matching using Variable Circular Template for Multi-resolution Image Registration

다중 해상도 영상 등록을 위한 가변 원형 템플릿을 이용한 특징 정합

Abstract

Image registration is an essential process for image fusion, change detection and time series analysis using multi-sensor images. For this purpose, we need to detect accurately the difference of scale and rotation between the multi-sensor images with difference spatial resolution. In this paper, we propose a new feature matching method using variable circular template for image registration between multi-resolution images. The proposed method creates a circular template at the center of a feature point in a coarse scale image and also a variable circular template in a fine scale image, respectively. After changing the scale of the variable circular template, we rotate the variable circular template by each predefined angle and compute the mutual information between the two circular templates and then find the scale, the angle of rotation and the center location of the variable circular template, respectively, in fine scale image when the mutual information between the two circular templates is maximum. The proposed method was tested using Kompsat-2, Kompsat-3 and Kompsat-3A images with different spatial resolution. The experimental results showed that the error of scale factor, the error of rotation angle and the localization error of the control point were less than 0.004, $0.3^{\circ}$ and one pixel, respectively.

다중 센서 영상을 영상 융합, 변화 탐지, 시계열 분석에 활용하기 위해서는 두 영상 간의 영상 등록 과정이 필수적이다. 영상 등록을 위해서는 서로 다른 공간 해상도를 가지는 다중 센서 영상 사이의 스케일과 회전각도 차이를 정확히 검출해야 한다. 본 논문에서는 다중 해상도 영상 간의 영상 등록을 위하여 가변 원형 템플릿을 이용한 새로운 특징 정합 기법을 제안한다. 제안하는 정합 기법은 스케일이 작은 영상의 특징점을 중심으로 원형 템플릿을 설정하고 스케일이 큰 영상에서는 가변 원형 템플릿을 생성한다. 가변 원형 템플릿의 스케일을 일정한 스케일 단위로 변경한 후에 가변 원형 템플릿을 일정 각도 단위로 회전시키면서 두 원형 템플릿 사이의 상호 정보량이 최대가 될 때의 가변 원형 템플릿의 스케일, 회전 각도 그리고 중심 위치를 각각 검출한다. 제안한 방법을 서로 다른 공간 해상도를 가지는 Kompsat(Korea Multi-Purpose Satellite) 2호, 3호, 3A호 영상 조합에 적용한 결과, 스케일 팩터 오차는 0.004 이하, 회전 각도 오차는 $0.3^{\circ}$ 이하, 제어점의 위치 오차는 1 화소 이하의 정합 성능을 보였다.

1. 서론

원격 탐사 영상을 이용한 변화 탐지, 시계열 분석과 같이 동일 대상을 다른 촬영 조건에서 획득한 영상을 비교하기 위해서는 두 영상 간의 영상 좌표를 일치시키는 영상 등록 과정이 필수적이다. 영상 등록 과정은 두 영상 좌표의 변환 기준점이 되는 특징점 검출과 두 영상에 공통으로 존재하는 특징점 쌍을 구하는 특징점 정합이 핵심적인 과정이다. 특징점 검출은 템플릿의 중심 화소를 기준으로 밝기값의 변화량을 측정하여 이웃 화소보다 변화량이 큰 화소를 특징점으로 검출한다(Moravec,1977). 이러한 특징점의 기하학적인 형태는 코너 형태가 일반적이며 그레디언트 행렬의 고유치를 이용하여 코너를 검출하는 다양한 방법들이 제안되었다(Harris and Stephens, 1988; Shi and Tomasi, 1994; Carneiro and Jepson, 2002; Ye, 2014a). 한편, 두 영상에서 검출된 특징점 간의 정합을 위해서 특징점 주변의 영상 정보를 담고 있는 영상 기술자(descriptor)를 이용한다. 영상 기술자는 특징점 주변의 영상 정보를 그 성분으로 하는 벡터이며 대표적으로 Scale Invariant Feature Transform(SIFT)(Lowe, 2004)가 널리 알려진 영상 기술자이다. SIFT 기반의 방법들은 영상 기술자가 특징점 주변의 영상 특성을 반영하는 히스토그램을 기본적으로 이용하기 때문에 스케일 및 회전 변화 등과 같은 기하학적 변형에 대해 강건한 특징을 가진다(Song et al., 2014; Gong et al.,2013; Sun et al., 2014). 그러나 SIFT 방법은 blob-like 형태의 특징점 검출에는 유리하지만 건물이나 도로의 교차점 등과 같은 코너점의 검출에는 적합하지 않은 단점이있다(Wu et al., 2012). SIFT 방법은 검출한 특징점의 위치가 실제 코너 위치가 아닌 코너 주변에서 검출이 되기 때문에 영상 등록을 위한 코너 검출기로서는 적합하지 않다.

Shechtman and Iran(2007)는 특징점을 중심으로하는 영상 템플릿과 특징점 주변의 영상과의 유사도를 계산하여 correlation surface를 생성한 후에 이를 극좌표 형식으로 표현하는 Local Self-Similarity(LSS) 기술자를 제안하였다. 영상의 밝기값 변화에 대해 강건한 특성을 가지는 LSS 기술자의 성능을 일부 개선한 방법도 제안이 되었다(Liu et al., 2012; Liu and Zeng, 2012). LSS 기술자의 단점은 다양한 특징에 대한 기술자의 식별 능력(discriminability)이 상대적으로 낮다는 점이다(Sedaghat and Ebadi, 2015). Correlation surface를 극좌표 형식으로 표현한 후에 극좌표 중심에서 일정 길이의 반지름 이내에 속하는 원형의 영역을, 중심에서의 거리와 방향에 따라 각각 일정한 개수의 작은 영역으로 분할한 후에 각영역의 대표 값으로 correlation surface 최대을 할당하여 LSS 기술자 벡터를 생성한다. 이러한 기술자 벡터생성 방식은 SIFT 기반의 기술자와 같이 특징점 주변의 정보를 벡터 형태로 표현하는데 유용한 방식이나, 특징점 주변 영역을 일정한 크기의 작은 영역으로 분할하고 대표값을 지정하는 과정을 통해서 밝기값 또는 에지 정보 등의 작은 변화에 대한 기술자의 민감도가 저하되는 요인으로 작용한다. 즉, 특징점을 중심으로 생성된 LSS벡터와 특징점에 가까운 화소의 LSS 벡터는 거의 유사한 벡터 성분을 가지게 되기 때문에 특징점과 특징점 근처의 화소에 대한 기술자의 식별 능력은 감소하게 된다.

기존의 영상 기술자를 이용한 정합 기법은 특징점을 중심으로 좁은 영역의 영상 정보를 특징 벡터로 기술하고 이 특징 벡터와의 유사도가 높은 특징 벡터를 가지는 대응점을 정합 과정에서 찾는다. 이 과정에서 특징 벡터가 추출되는 영역 자체가 작으면 특징점 주변의 영상 정보가 충분히 반영되지 못하는 문제가 발생될 수 있으며, 또한 특징 벡터를 추출하는 과정 자체가 특징점 주변의 영상 전체 밝기값 정보 대신에 벡터로 추출된 일부영상 정보만을 이용하는 근본적인 한계를 가진다. 반면, 특징점 기술에 사용되는 주변 영역을 충분히 넓게 하면 특징점 정합을 위한 유사도 계산 과정에서 잡음이나 다른 요인에 의해 유사도를 저하시키는 일부 요인이 작용하더라도 그 영향을 적게 받을 수 있다. Ye(2014b)의 연구에 따르면 정합에 사용되는 원형 템플릿의 크기가 점진적으로 증가하면 정합오차도 함께 감소하고, 원형 템플릿의 반지름이 100~500 화소의 비교적 넓은 범위로 증가하면 정합오차도 1화소 이내의 낮은 수준으로 수렴하였다. 이와 같은 결과는 특징점 중심의 좁은 영역의 유사도 비교에 비해 넓은 영역의 유사도 비교가 가지는 장점을 잘 보여준다.

앞서 소개한 기존의 영상 기술자 중심의 특징점 정합기법들은 스케일 및 회전 변화를 포함하는 영상 간의 특징점 유사도를 비교하기 위해서 스케일 변화, 회전 변화의 영향을 적게 받는 특징점 추출 방법, 특징점 유사도 계산 방법을 제시한다. 그러나 영상 등록을 위한 이러한 특징점 중심의 유사도 계산 접근 방식 대신에 특징점을 포함하는 넓은 영역 간의 영상 변환(image transformation)관점에 본다면 영상 등록은 두 비교 영상 간의 어파인 파라미터를 결정하는 과정으로 이해할 수 있다. 두 영상 간의 일반적인 어파인 파라미터를 결정하기 위해서는 최소 3개 이상의 대응점이 필요하지만 두 영상 내부의 정합 특징점을 중심으로 하는 두 영상 간 변환 관점에서 보면 두 영상 간의 어파인 변환은 두 영상 간의 정합 특징점을 변환 중심으로 하는 영상 간의 스케일 변환(scale transformation), 회전 변환(rotation transformation) 그리고 이동 변환(translation transformation)의 각 파라미터를 결정하는 과정과 등가의 의미를 가지게 된다. 스케일 변환, 회전 변환, 이동 변환의 각 파라미터를 구하게 되면한 쪽 영상의특징점을 중심으로 각각의 영상변환을 수행하여 다른 한쪽의 영상 좌표와 일치시킬 수 있다. 이 과정은 그 자체로서 특징점을 중심으로 하는 일종의 글로벌 영상 등록(global image registration)으로 볼 수 있으며 삼각망기반의 국부 변환(triangulation-based local transformation)기법(Ye, 2014c) 등과 같이 보다 정밀한 영상 등록을 위한 초기 영상 등록 방법으로 사용할 수 있다. 두 영상 사이의 스케일 및 회전 차이가 보정되면, 두 영상 간의 기하학적인 차이가 보정이 되어 추가적인 정밀 영상 등록 과정에서 유사도 계산의 정확도를 높일 수 있다. 이와 같은 접근 방법에서 해결해야 하는 문제는 두 영상 간의 정합특징점은 일반적으로 미리 알려져 있지 않기 때문에 각 어파인 파라미터를 구하는 과정과 정합 특징점을 구하는 과정이 동시에 진행되어야 한다.

본 논문에서는 가변 원형 템플릿 정합 기법(Ye, 2018)을 기초로 하여, 기존의 특징점 검출 및 정합 방식을 대체할 수 있는 어파인 파라미터를 이용한 다중 해상도영상의 특징 정합 기법을 소개한다. 먼저 영상 간의 기하학적 보정이 특징점 정합 유사도 향상에 미치는 영향을 분석하고 스케일이 서로 다른 영상 간의 원형 템플릿을 이용한 영상 변환 방법을 제시한다. 이어서 특징점의 유사도 계산을 위한 원형 템플릿의 설계 방법과 원형 템플릿 간의 상호 정보량에 기반한 유사도 계산 방법을 제시한다. 제안한 정합 방법의 성능 검증을 위해서 Kompsat 3호 영상을 이용하여 다양한 조건에서 생성된 스케일 영상들과 실제 서로 다른 스케일을 가지는 Kompsat 2호, 3호, 3A호 영상을 이용하여 성능 평가를 각각 수행하였다.

2. 연구 방법

1) 기하학적 보정과 유사도 향상

이종 센서 영상의 영상 취득 조건의 차이로 인해서 발생하는 기하학적 차이는 영상 등록 과정에서 특징 사이의 유사도를 저하시키는 주요 원인이 된다. Fig. 1(a)와 (b)는 동일 지역을 2002년 2월 7일에 촬영한 11비트 방사해상도와 1m 공간해상도를 가지는 전정색 Ikonos위성 영상과 2005년 1월 15일에 촬영한 11비트 방사해상도와 60 cm 공간해상도를 가지는 전정색 Quickbird위성 영상의 예를 보여준다. Fig. 1(c)는 Ikonos 위성 영상을 기준으로 Quickbird 위성 영상을 어파인 변환한 영상으로, Fig. 1 하단의 확대된 영상에 보듯이 어파인 변환 후의 Quickbird 위성 영상과 Ikonos 위성 영상과의 기하학적 차이가 감소한 것을 볼 수 있다. 어파인 변환을 통해 Fig. 1(c)의 건물 코너 중심에 연결된 두 개의 흰색 화살표 사이의 각도가 Fig. 1(a)의 Ikonos 영상의 코너 화살표 사이의 각도와 유사하게 변환된 것을 통해 기하학적 차이의 감소를 쉽게 확인 할 수 있다.

Fig. 1. Ikonos image and Quickbrd images beforeand after affine transformation. (a) Ikonos image. (b) Quickbird image before affine transformation. (c) Quickbird image after affine transformation.

어파인 변환을 통한 기하학적 보정을 통해 영상 등록을 위한 특징점 정합 과정에서 정합 성능의 향상을 기대할 수 있다. Fig. 2는 Fig. 1(a)의 흰색 원의 중심 위치에 있는 Ikonos 위성 영상의 코너 화소에 대한 대응점을 어파인 변환 전후의 Quickbird 위성 영상에서 찾기 위하여 상호 정보량(mutual information)을 이용하여 유사도를 계산한 결과를 보여준다. 정합 화소의 상호 정보량 최대값은 어파인 변환 적용 전후의 값이 각각 0.514, 0.613이 얻어져서 변환 적용 후의 유사도가 보다 높은 결과를 보인다. 특히, 정합점과 주변 화소들 간의 유사도의 차이는 Fig. 2(b)와 같이 어파인 변환 후에 그 차이가 크게 증가하여 오정합 가능성이 크게 낮아 질 수 있음을 볼 수 있다. 따라서 두 영상 간의 기하학적 차이를 보정하면 정합점 간의 유사도를 높이고 정합점과 주변 화소들 간의 유사도 차이를 증가시켜서 정합 성능을 향상시킬 수 있음을 알 수 있다. 두 영상 간의 기하학적 보정을 위해서는 일반적으로 두 영상 간의 특징점들을 정합한 후에 정합된 특징점들을 이용하여 영상 변환을 수행하는 단계를 거치나 본 연구에서는 원형 템플릿을 이용한 영상 변환과 특징점 정합을 동시에 수행하는 방법을 제안한다.

Fig. 2. Comparison of mutual information-based similarities between Ikonos image and before and after Quickbird images of affine transformation. (a)mutual information before affine transformation (b)mutual information after affine transformation.

2) 원형 템플릿을 이용한 영상 변환

스케일이 서로 다른 영상 간의 특징점 정합은 두 영상 간의 대응점(corresponding points)을 연결하는 어파인 변환 성분을 찾는 과정으로 해석할 수 있다. 상대적으로 스케일이 작은 영상의 대응점을(x_l, y_l), 스케일이 큰 영상의 대응점을(x_h, y_h), 대응점(x_h, y_h)의 스케일이 큰 영상 좌표 상의 어파인 변환 후의 좌표를(x_t, y_t)라할 때, 대응점(x_h, y_h)과 변환 후 좌표(x_t, y_t) 이의 관계는 일반적으로 식 (1)과 같은 어인 변환 관계로 주어지며 이를 스케일 성분 s, 회전각 θ, 이동 변위(d_x, d_y)로 나타내면 식 (2)와 같이 표현된다. 스케일이 작은 영상의 대응점(x_l, y_l)과의 유사도가 가장 높은 어파인 변환 후의 좌표(x_t, y_t)를 구하면 어파인 변환 계수를 이용하여 스케일이 큰 영상에서의 대응점(x_h, y_h)를 구할 수 있다.

\(\begin{array}{l} x_{t}=a_{1} \cdot x_{h}+a_{2} \cdot y_{h}+a_{3} \\ y_{t}=a_{4} \cdot x_{h}+a_{5} \cdot y_{h}+a_{6} \end{array}\)       (1)

\(\left[\begin{array}{l} x_{t} \\ y_{t} \end{array}\right]=\left[\begin{array}{ll} s & 0 \\ 0 & s \end{array}\right]\left[\begin{array}{cc} \cos \theta-\sin \theta \\ \sin \theta & \cos \theta \end{array}\right]\left[\begin{array}{l} x_{h} \\ y_{h} \end{array}\right]+\left[\begin{array}{l} d_{x} \\ d_{y} \end{array}\right]\)       (2)

대응점(x_l, y_l)과 어파인 변환 후 좌표(x_t, y_t)와의 유사도를 측정하기 위하여 스케일이 작은 영상의 대응점(x_l, y_l)과 스케일이 큰 영상의 후보 대응점(x_h, y_h)을 중심으로 원형 템플릿을 각각 생성한다. 원형 템플릿 내에 위치하는 화소들은 두 대응점 간의 유사도 계산에 사용되며 대응점을 중심으로 영상이 회전하여도 원형 템플릿 내부에 포함되는 화소들은 변하지 않는 회전 불변(rotation invariance) 특성을 가진다. 따라서 대응점(x_l, y_l) 중심의원형 템플릿은 고정하고 스케일이 큰 영상의 후보 대응점(x_h, y_h) 중심의 원형 템플릿은 일정 각도 단위로 회전시키면서 두 템플릿의 유사도가 최고가 되는 회전 각도를 측정하여 두 템플릿의 상대적인 회전 각 차이, 즉 두 영상의 상대적인 회전 각도 차이를 구한다.

한편, 스케일이 작은 영상과 스케일이 큰 영상의 경우 상대적인 스케일 차이가 존재하기 때문에 두 원형 템플릿을 동일한 크기로 설정을 하면 스케일이 큰 영상의원형 템플릿은 스케일이 작은 영상의 템플릿보다 실제로는 더 작은 영역을 포함하게 된다. 따라서 두 원형 템플릿의 유사도를 비교하기 전에 먼저 두 원형 템플릿의 스케일을 일치시키는 과정이 필요하다. 이를 위해서 Fig. 3에 보인 바와 같이 스케일이 작은 영상 플릿 내의 위치 \(r_{i}^{l}\)에 대응되는 스케일이 큰 영상 템플릿 내의 위치 \(r_{i}^{h}\)를 스케일 파라미터 s_k를 이용하여 \(r_{i}^{h}=s_{k} \cdot r_{i}^{l} \)로 계산한 후에 두 원형 템플릿 반지름 상에 동일한 개수의 화소를 배치한다. 원형 템플릿의 반지름 상에 동일 개수의 화소를 배치하면 스케일 파라미터 s_k의 변화에도 불구하고 두 원형 템플릿의 유사도 계산에 사용되는 화소들의 개수를 일치시킬 수 있다.

Fig. 3. Affine transformation of the fine-scale image and similarity measurement between the affine transformed fine-scale image and coarse-scale image. The similarity is computed using the circular template of the coarse-scale image and the circular template obtained after each parameter combination of the affine transformation in the fine-scale image.

3) 원형 템플릿의 설계와 유사도 계산

본 논문에서는 원형 템플릿에 속하는 화소의 위치(x_ij, y_ij)를 결정하기 위하여 Fig. 4와 같이 극좌표(polar coordinates) 방식을 이용하여 스케일이 작은 영상에서 중심 좌표가 (x_l, y_l)이고 원점에서의 거리가 \(r_{i}^{l}\), 각도가 θ_j인 화소의 위치를 식 (3)과 같이 표현한다. 이와 유사한 방식으로 스케일이 큰 영상에서 대응되는 화소의 위치(x_ij, y_ij)는 식 (4)와 같이 주어진다.

\(\begin{array}{l} x_{i j}=x_{l}+r_{i}^{l} \cos \theta_{j} \quad\left(0 \leq i \leq N_{r}, 0 \leq j       (3)

\(\begin{array}{l} x_{i j}=x_{h}+s_{k} \cdot r_{i}^{h} \cos \left(\theta_{j}+\theta_{d}\right) \\ y_{i j}=y_{h}+s_{k} \cdot r_{i}^{h} \sin \left(\theta_{j}+\theta_{d}\right) \end{array}\)       (4)

Fig. 4. Determination of the location of pixels in the circular templates using the parameters of rotation angle θ_j, sale factor s_k and displacement (x_h, y_h) in (a) the coarse-scale image and in (b) the fine-scale image.

여기서 N_r은 원형 템플릿의 반지름 위에 위치하는 전체 화소수, N_θ는 원형 템플릿의 원주 위에 위치하는 전체 화소수를 의미한다. Fig. 4는 N_r=8, N_θ=8인 원형 템플릿의 생성 예를 보여준다. 0 ≤ i ≤ N_r, 0 ≤ j ≤ N_θ 사이의 모든 화소의 위치를 거리 ri와 각도 θ_j의 2차원 좌표평면에 도시하면 Fig. 5와 같이 2차원 배열 형태로 표시된다. 수평 방향의 축은 템플릿 중심에서 방사적으로 떨어진 위치 r_i를 나타내고, 수직 방향의 축은 반시계 방향의 회전 각도 θj를 나타낸다.

Fig. 5. Two-dimensional arrays for the pixels of circular templates in (a) the coarse-scale image and in (b) the fine-scale image after affine transformation.

스케일이 작은 영상과 스케일이 큰 영상의 (x_ij, y_ij) 위치에 존재하는 화소의 밝기값을 각각 I_l(x_ij, y_j), I_h(x_ij, y_ij)라 하면 스케일이 작은 영상과 스케일이 큰 영상의 원형 템플릿에 포함되는 화소들의 집합 P_l과 P_h는 식 (5)와 식 (6)과 같이 주어진다.

\(P_{l}=\left\{l\left(x_{i j}, y_{i j}\right) \mid x_{i j}=x_{l}+r_{i}^{l} \cos \theta_{j}, y_{i j}=y_{l}+r_{i}^{l} \sin \theta_{j}\right\}\)       (5)

\(\begin{aligned} P_{h}\left(k, \theta_{d}, d_{x}, d_{y}\right)=\{& I_{h}\left(x_{i j}, y_{i j}\right) \mid \\ & x_{i j}=x_{h}+s_{k} \cdot r_{i}^{h} \cos \left(\theta_{j}+\theta_{d}\right), \\ &\left.y_{i j}=y_{h}+s_{k} \cdot r_{i}^{h} \sin \left(\theta_{j}+\theta_{d}\right)\right\} \end{aligned}\)       (6)

여기서 x_h = x_l + d_x, y_h = y_l + d_y이며 k는 스케일이 큰 영상의 스케일 범위를 나타나는 인덱스로의 0 ≤ k ≤ N_s범위를 가진다. θ_d는 스케일이 큰 영상에서 생성된 템플릿의 회전 각도를 나타내며 식 (7)과 같이 주어진다. 기본값으로 T=180이 주어지며 회전 각도를 보다 미세하게 조정하기 위해서는 T를 보다 값으로 조정한다.

\(\theta_{d}=\frac{\pi}{T} \cdot d(1 \leq d<360)\)       (7)

스케일이 작은 영상과 스케일이 큰 영상의 원형 템플릿에 포함되는 화소들의 집합 P_l과 P_h의 유사도를 측정하기 위하여 식 (8)과 같이 주어지는 상호 정보량 mI를계산한다.

\(\begin{array}{c} M I\left(P_{l}, P_{h}\left(k, \theta{d}, d_{x}, d_{y}\right)\right)= \\ H\left(P_{l}\right)+H\left(P_{h}\left(k, \theta_{d}, d_{x}, d_{y}\right)\right)-H\left(P_{l}, P_{h}\left(k, \theta_{d}, d_{x}, d_{y}\right)\right) \end{array}\)    (8)

여기서 H는 엔트로피를 나타낸다. 상호 정보량 MI를 계산하기 위하여 먼저 스케일이 작은 영상의 특징 좌표(x_l, y_l)를 검출하기 위해서 일정 범위 내에서 코너 응답값이 최대가 되는 코너점의 위치를 기하 평균 기반의 특징점 검출기를 이용하여 검출한다(Ye, 2014a). 검출된 특징 좌표(x_l, y_l)를 중심으로 반지름이 R인 원형 템플릿에 속하는 화소 집합 Pl을 설정한 후에 스케일이 큰 영상의 스케일 0 ≤ k ≤ N_s 범위에 속하는 각 스케일 k에 대해 식 (6)을 이용하여 원형 템플릿 화소 집합 P_h를 설정한다. 이때 원형 템플릿을 Fig. 6과 같이 θ_d 만큼 순차적으로 증가시키면서 식 (6)을 이용하여 (x_ij, y_ij) 위치를 갱신한다. 이 과정을 통해 초기에 만들어진 원형 템플릿 화소 집합 P_h는 갱신된 (x_ij, y_ij)에 위치하는 화소 집합으로 재구성된다. 재구성된 원형 템플릿의 중심 위치(x_h, y_h)를 일정한 탐색 범위 내의 d_x와 d_y 값에 따라 변경하면서 식 (9)와 같이 상호 정보량 MI가 최대값이 될 때의 각 파라미터(\(k^{*}, \theta_{d}^{*}, d_{x}^{*}, d_{y}^{*}\))의 값을 구한다.

\(\left(k^{*}, \theta_{d}^{*}, d_{x}^{*}, d_{y}^{*}\right)=\operatorname{argmax} M I\left(P_{l}, P_{h}\left(k, \theta_{d}, d_{x}, d_{y}\right)\right)\)       (9)

Fig. 6. Updating the location of each pixel of the circular template by changing the angle θ_d. (a) Initial pixel locations. (b) Updated pixel locations

상호 정보량 계산에 소요되는 시간을 줄이기 위하여 원형 윈도우 생성시 템플릿 중심에서 방사적으로 떨어진 위치 r_i 사이의 간격과 회전 각도 θ_j 사이의 간격을 각각 일정한 간격을 두고 생성한다. 원형 템플릿의 반지름 R이 300 화소이면 총 N_r=300개의 r_i를 생성할 수 있으나, 3 화소 간격으로 화소 값을 취하면 N_r=100개의 r_i를 생성하게 된다. 이 경우 상호 정보량 계산에 이용되는 화소의 수가 1/3로 감소하여 연산 시간도 함께 감소한다. 회전 각도 θ_j의 간격도 1° 간격으로 생성하면 각 r_i위치마다 N_θ=360개의 화소를 생성할 수 있으나 5° 간격으로 생성하면 N_θ=72개의 화소가 생성되어 상호 정보량 계산에 이용되는 화소의 수가 1/5로 감소하게 된다. 또한 스케일이 큰 영상에서의 정합 위치(d_x, d_y)를 탐색하는데 소요되는 시간을 줄이기 위하여 스케일이 큰 영상의 전체 화소 그레디언트의 히스토그램에서 상위 임계 비율 T_h 이상인 화소를 정합 후보점으로 한정하면 탐색 시간을 크게 단축할 수 있다.

3. 실험 결과 및 분석

본 논문에서는 Table 1과 같이 세 가지 유형의 영상 자료를 이용하여 제안한 알고리듬의 성능을 검증하였다. 먼저 A 유형은 Kompsat-3 영상을 스케일이 큰 영상으로 삼고 스케일팩터 s_k의 범위가 1.2~3.0이 되도록 Kompsat3의 해상도를 변환하고 회전 각도 θ_d={0°, 90°, 180°,270°}인 4종류의 영상을 생성한 후에 10%의 균일 잡음(uniform noise)를 추가하여 총 40장의 작은 스케일 영상을 사용하였다. 두 번째 B 유형은 넒은 범위의 스케일팩터(2.0≤s_k≤10.0)와 회전 각도를 적용한 총 36장의 작은 스케일 영상을 사용하였다. 세 번째 C 유형은 동일한 지역을 촬영한 Kompsat-2, Kompsat-3, Kompsat-3A 영상을 조합하여 사용하였다.

Table 1. Three types of the experimental images

본 논문에서 사용된 실험 영상은 Niger의 Niamey를 촬영한 Kompsat 2호, 3호, 3A호 영상으로 공간해상도가 각각 1.0m, 0.7m, 0.55m이다(Fig. 7). 제안한 방법의 성능평가를 위하여 A 유형의 Kompsat-3 영상에 대해 Fig. 8과 같이 기하 평균 기반의 코너 검출기를 이용하여 코너 응답이 최대가 되는 제어점을 검출한 후에 다양한 스케일 팩터와 회전 성분을 가지면서 잡음 성분이 추된 Kompsat-3 시뮬레이션 영상을 생성하였다(Fig. 9).

Fig. 7. Experimental Kompsat images of Niamey. (a) Kompsat-2 image (April 18, 2009). (b) Kompsat-3 image (Jan. 21, 2013).(c) Kompsat-3A image (Nov. 29, 2015).

Fig. 8. Result of corner detection from Kompsat-3 image. (a) Kompsat-3 panchromatic image. (b) Centers of circles represent the corners detected using geometric mean-based corner detector. The center of white crosshair represents the corner with maximum corner response, which is used as the control point in the estimation of affine parameters

Fig. 9. Some sets of the tested Kompsat-3 images. (a) Original image. (b) Image rescaled by the scaling factor s_k=2 and corrupted by 10% uniform noise. (c) Image rescaled by the scaling factor s_k=3 and rotated by the angle of 90° under 10% uniform noise. (d) Image rescaled by the scaling factor s_k=3and rotated by the angle of 180° under 10% uniform noise.

스케일 팩터는 Table 2와 같이 1.2부터 3.0까지 0.2 간격으로 생성하고 회전 성분은 0°, 90°, 180°, 270°로 생성하였다. 잡음 성분은 생성된 모든 영상에 대해 10%의 균일 잡음을 추가하였다. Kompsat-3 영상과 생성된 스케일이 작은 다양한 Kompsat-3 영상 사이의 어파인 파라미터 \(\left(k^{*}, \theta_{d}^{*}, d_{x}^{*}, d_{y}^{*}\right)\)를 구하기 위하여 반지름 R=300인 원형 템플릿에 대해, N_r=100, N_θ=72인 원형 템플릿을스케일이 작은 Kompsat-3 영상의 제어점 위에 생성하였다. 스케일이 큰 영상인 Kompsat-3 영상의 제어점 후보는 그레디언트 히스토그램의 상위 임계 비율 T_h=0.05상인 화소로 한정하였다. 스케일 팩터 s_k는 0.1 간격으로, 템플릿의 회전 각도 θ_d는 0.1° 간격으로 변화시키면서 다양한 조건에서 생성된 각 영상과의 상호 정보량 유사도가 최대가 될 때의 스케일 팩터 s_k, 회전량 θ_d, 정합 위치(d_x, d_y)를 각각 구하였다. Table 2는 다양한 스케일 팩터 및 회전각 조건에서 검출한 스케일 팩터를 보여주며 총 40가지 실험 조건에서 스케일 팩터를 모두 정확하게 검출하였다. Table 2는 스케일이 작은 영상과 스케일이 큰 영상의 상대적인 회전각들을 모든 스케일에서 정확하게 검출한 결과를 보여준다.

Table 2. Results of detection of scale factors s_k for the coarse-scale images

Table 3. Results of detection of rotation angles

Table 4는 검출된 제어점의 실제 기준점과의 위치 오차를 보여준다. 검출된 제어점의 행과 열 방향 오차의 평균값을 각각 구한 것으로 각 스케일의 4가지 회전각에 대해 0.50 ~1.62 화소 범위의 평균 오차 결과를 얻었다. 스케일 팩터가 증가할수록 두 영상의 상대적인 스케일차이가 증가하게 되어 스케일 팩터 s_k=3인 경우, 스케일이 큰 영상에서 3×3 크기의 화소 영역은 스케일이 작은 영상에서 1화소 크기의 영역에 해당된다. 따라서 스케일이 큰 영상에서의 위치 오차는 스케일 팩터 증가에 따라 비례적으로 증가하기 때문에 위치 오차 분석에서 스케일 팩터의 영향을 배제하기 위하여 평균 오차와 스케일 팩터의 비(ratio)를 구하였다. 이 값은 스케일이 작은 영상 관점에서의 위치 오차로 이해할 수 있으며 1.2~3.0 사이의 스케일 팩터에 대해 0.42~0.69 화소 범위의 우수한 위치 오차 성능을 보였다.

Table 4. Error in the localization of the control point

Fig. 10(a)는 스케일이 작은 영상에 대한 스케일이 큰 영상의 회전량 &thea;_d가 0° 및 90°일 때, 스케일 팩터 s_k가 0.05 간격으로 변할 때의 상호 정보량의 변화를 보여준다. 상호 정보량의 최고 값은 인접한 스케일 팩터 s_k를 포함하여 다른 스케일 팩터의 상호 정보량과 비교해서 회전량 θ_d가 0° 및 90°인 두 경우 모두 큰 차이가 남을 볼 수 있다. Fig. 10(b)는 스케일 팩터 s_k=2.0일 때, 회전량 θ_d가 변할 때 상호 정보량의 변화를 보여준다. 회전량 θ_d가 90°일 때 상호 정보량이 최대가 되고 그 외의 다른 각도에서의 상보 정보량과 그 차이가 뚜렷함을 알 수 있다.

Fig. 10. Illustration of the detection of scale factor and rotation angle. (a) Detection of the scale factor for the coarse-scale image in the case of the scaling factor 2.0. (b)mutual information versus rotation angles in the case of the scaling factor 2.0 and the rotation angle of 90°under 10% uniform noise.

보다 넓은 범위의 스케일 팩터에 대한 실험을 위해 2.0≤s_k≤10.0 사이의 스케일 팩터와 4가지 회전 각도 θ_d={0°, 90°, 180°, 270°}를 적용한 유형 B의 저해상도 영상을 사용하였다. Fig. 11은 스케일 팩터가 2.0, 4.0, 6.0, 8.0, 10.0이고 회전 성분을 포함한 다양한 저해상도 영상의 예를 보여준다. 제안한 방법은 Table 5와 같이 각각의 스케일 팩터와 회전 각도 조건에서 모두 정확하게 스케일 팩터를 검출하였다. 두 비교 영상 간의 회전 각도도 Table 6과 같이 모든 조건에서 정확하게 검출하였다. Table 7은 검출된 제어점의 실제 기준점과의 위치 오차를 나타내며 스케일 범위 내에서 0.5~2.38 화소의 위치 오차를 보였다. 스케일을 고려한 평균 위치 오차와 스케일 팩터의 비는 0.24~0.42 범위의 값을 보였으며 스케일 팩터 s_k=10.0인 경우에도 평균 위치 오차와 스케일 팩터의 비가 0.24의 낮은 수준을 보였다. 이는 제안한 방법이 스케일이 작은 영상 관점에서 두 영상의 스케일이 2~10배 정도 차이가 나더라도 0.24~0.42 화소 수준의 높은 위치 정확도를 가짐을 보여준다.

Fig. 11. Some sets of the tested Kompsat-3 images. (a) Original image. (b) Scalng factor of 2.0 and rotation angle of 90°. (c) Scaling factor of 4.0 and rotation angle of 180°. (d) Scaling factor of 6.0and rotation angle of 270°. (e) Scaling factor of 8.0 and rotation angle of 90°. (f) Scaling factor of10.0 and rotation angle of 180°

Table 5. Results of detection of scale factors for the coarse-scale images

Table 6. Result of detection of rotation angles

Table 7. Error in the localization of control points

Fig. 12는 세 번째 C 유형의 영상 조합인 Kompsat-3 영상과 Kompsat-2 영상, Kompsat-3A 영상과 Kompsat-2 영상, Kompsat-3A 영상과 Kompsat-3 영상 간의 스케일팩터 변화에 따른 상호 정보량의 변화를 보여준다. 스케일 펙터의 변화 단위가 0.01인 경우에 각 영상 조합에서 상호 정보량이 최대가 될 때의 스케일 팩터는 각각1.43, 1.79, 1.25가 얻어졌다. Table 8부터 Table 10은 스케일 팩터의 변화 단위가 보다 정밀한 0.001인 경우의 각 영상 조합에 대한 결과를 보여준다. 스케일 팩터가 0.001 단위로 변하는 경우, 상호 정보량이 최대가 되는 스케일 팩터가 여러 값에 걸쳐서 나타났으며 최종 검출된 스케일 팩터는 상호 정보량이 최대가 되는 스케일 팩터들의 평균값을 취하고 소수 넷째 자리의 값을 반올림하여 추정하였다. Kompsat-3 영상과 Kompsat-2 영상 조합의 실험에서는 상호 정보량이 최대가 되는 스케일 팩터가 1.427~1.429 범위에서 나타났으며 Table 8과 같이 추정된 스케일 팩터는 기준 스케일 팩터와 동일한 1.428이 얻어졌다. 두 영상의 상대적인 회전 각도의 오차는 0.2°, 제어점의 위치 오차는 행과 열 방향으로 각각 1화소의 오차를 보였다. Kompsat-3A 영상과 Kompsat-2 영상 조합의 실험에서는 상호 정보량이 최대가 되는 스케일 팩터가 1.790~1.793 범위에서 나타났다. 추정된 스케일 팩터 1.792는 기준 스케일 팩터와 0.004의 오차를 보였다. 회전 각도 오차는 0.3°를 보였으며 제어점에 대응하는 기준점의 행과 열 방향의 위치는 정확하게 검출이 되었다(Table 9). Kompsat-3A 영상과 Kompsat-3 영상 조합의 실험에서는 상호 정보량이 최대가 되는 스케일 팩터는 1.254~1.257 범위에서 나타났으며 추정된 스케일 팩터 1.256은 기준 케일 팩터와 0.004의 오차를 보였다. 회전 각도 오차는 0.1°, 제어점에 대응하는 기준점의 행과 열 방향의 위치는 정확하게 검출이 되었다(Table 10).

Table 8. Results of image registration using Kompsat-2 and Kompsat-3 images

Table 9. Results of image registration using Kompsat-2 and Kompsat-3A images

Table 10. Results of image registration using Kompsat-3 and Kompsat-3A images

Fig. 12. mutual information versus scale factors. (a) Kompsat-3 image vs. Kompsat-2 image. (b)Kompsat-3A image vs. Kompsat-2 image. (c) Kompsat-3A image vs. Kompsat-3 image.

실험 결과에서 얻어진 스케일 팩터, 회전 각도, 특징점 위치 검출의 높은 정확도는 원형 템플릿의 크기 및 유사도 계산 방식과 밀접한 관계를 가진다. 실험에 사용된 원형 템플릿의 반지름 R은 300 화소이고 r_i를 3 화소단위로, 회전 각도 θ_j를 5° 간격으로 샘플링하여 N_r=100, N_θ=72로 총 7,200개의 화소 밝기값이 유사도 계산에 사용되었다. 반지름 R=300인 원형 템플릿이 실제로 포함하는 영역의 이론적인 면적은 대략 300×300×3.14=282,600 화소로 일반적인 특징 정합 방식에서 사용되는 특징의 수에 비해 매우 넓은 영역에 분포하는 화소들이 유사도 계산의 특징으로 사용이 된다. 따라서 유사도 비교시 특징점 주변의 좁은 영역이 아닌 넓은 영역에 걸쳐서 유사도 비교가 수행이 되기 때문에 일부 영역에서 정합 성능을 떨어뜨리는 요인이 작용하더라도 전체 유사도 계산 과정에서의 영향은 미미하게 작용하고 오정합 가능성을 크게 낮추는 효과를 얻을 수 있다. 넓은 영역을 포함하는 원형 템플릿의 이러한 특징은 스케일 팩터와 회전 각도를 매우 정밀하게 조절할 수 있도록 한다. 또한 상호 정보량 기반의 유사도 측정 방식은 화소의 밝기값을 직접 비교하는 상관(correlation) 관계 기반의 유사도 측정 방식과 달리 두 원형 템플릿에 포함되는 화소들을 통계적으로 분석한다. 따라서 제안하는 방법은 좁은 영역에 대한 상호 정보량 기반의 유사도 측정 대신에 넓은 영역에 대한 상호 정보량 기반의 유사도 측정 방식을 적용하기 때문에 비교 대상이 되는 두 원형 템플릿 내부 영역에서 나타나는 미세한 차이를 넓은 영역에 걸쳐 통계적으로 측정할 수 있는 장점을 가진다.

4. 결론

최근 다양한 종류의 고해상도 위성 영상이 등장하면서 이들 다중 센서 영상을 다양한 목적으로 함께 활용하고자 하는 연구가 활발하다. 다중 센서 영상을 영상융합, 변화 탐지, 시계열 분석 등에 활용하기 위해서는 두 영상 간의 영상 등록 과정이 필수적이다. 서로 다른 공간 해상도를 가지는 위성 센서 영상 사이의 영상 등록을 위해서 SIFT 기반의 특징 정합 방법을 적용할 수 있으나 검출한 특징점의 위치가 코너 주변에서 검출이 되는 문제와 함께 특징 벡터를 추출는 과정을 통해 원래의 영상 밝기값 정보 대신에 일부 특징 정보만을 이용하는 근본적인 한계를 가진다. 본 논문에서는 해상도 차이가 존재하는 두 영상 간의 영상 등록 문제를 스케일이 서로 다른 두 영상 간의 어파인 변환(affine transformation)관점에서 재해석하여 다중 해상도 영상 간의 상대적인 스케일, 회전 각도 차이, 특징점 위치를 동시에 검출할 수 있는 가변 원형 템플릿을 이용한 정합 기법을 제안하였다. 제안한 가변 원형 템플릿을 이용한 정합 기법은 스케일이 작은 영상의 특징점을 중심으로 원형 템플릿을 설정하고 스케일이 큰 영상에서 가변 원형 템플릿의 스케일을 일정한 스케일 단위로 변경하면서 스케일이 작은 영상의 원형 템플릿과의 상호 정보량이 최대가 될 때의 가변 원형 템플릿의 스케일, 회전 각도 그리고 가변 원형 템플릿의 중심 위치를 각각 검출하였다.

Kompsat-3 영상을 스케일이 큰 영상으로 삼고 스케일 팩터의 범위가 1.2~3.0이 되도록 스케일이 작은 Kompsat-3 영상을 생성한 후에 4 종류의 회전 각도와 10%의 균일 잡음을 적용한 총 40가지 실험 조건에서 스케일 팩터, 상대적인 회전 각도를 모두 정확하게 검출하였으며 스케일이 작은 영상 관점의 위치 오차인 평균 오차와 스케일 팩터의 비가 0.42~0.69인 우수한 위치 오차결과를 얻었다. 스케일 팩터 범위가 보다 넓은 2.0~10.0 사이에서 4 종류의 회전 각도 조건을 적용한 총 36가지 실험 조건에서도 스케일 팩터, 상대적인 회전 각도를 모두 정확하게 검출하였고 스케일을 고려한 평균 위치 오차와 스케일 팩터의 비는 0.24~0.42 화소 범위의 결과를 얻었다. 서로 다른 공간 해상도를 가지는 Kompsat 2호,3호, 3A호 영상 조합에 대한 실험에서는 스케일 팩터 오차는 0.004 이하, 회전 각도 오차는 0.3° 이하, 제어점의 위치 오차는 1 화소 이하의 정합 성능 결과를 얻었다.

기존의 영상 등록 방법이 스케일 및 회전에 불변한 특징점 검출 및 정합 관점에서 다양한 영상 기술자와 정합 방식을 제시하는데 반해서 제한하는 방법은 특징점 주변의 넓은 영역에 대한 영상 변환이라는 새로운 관점에서 영상 등록 문제를 해결하고자 하였다. 영상 변환관점의 새로운 접근 방법을 통해 서로 다른 해상도를 가지는 영상 간의 상대적인 스케일과 회전 그리고 특징점 위치를 정밀하게 획득할 수 있었다. 상대적인 스케일 검출과 관련하여 기존의 연구는 두 영상 내의 존재하는 특징을 중심으로 특성 스케일(characteristic scale)을 각각 검출한 후에 그 크기를 비교하여 상대적인 스케일을 추정하나 이 경우 먼저 대응하는 두 특징점을 찾아야 하는 특징점 정합 문제로 인해서 자동 영상 등록에 적용하기에는 어렵다. 본 논문은 일차적으로는 다중 해상도 영상 간의 상대적인 스케일을 자동으로 검출하는 방법을 제시하였으며, 이와 함께 두 영상 간의 상대적인 회전 각도 차이와 특징점 위치를 정확하게 검출하여 다중해상도 영상 간의 특징 정합방법을 제시한 데 그 의의가 있다.