SIFT 원리, 구현

해당 포스트는 "열혈강의 영상처리 프로그래밍" 책의 내용을 요약한 것이다.

사정상 이전 포스트 "NCC, 교차 검토, 특징 기술자" 포스트에서 SIFT 특징 기술자 생성기를 설명하다가 중단했다. 이번 포스트에서는 SIFT 특징 기술자 생성기의 원리에 대해서 알아본다. SIFT 특징 기술자를 이용한 영상 정합을 해보겠다. 

※ SIFT 특징 검출기

SIFT 특징 검출기의 각 피라미드 단계에서 특징점 검출에 사용하는 정보는 해리스 코너 검출기와 유사한데, 기본은 영상에 대한  2차 미분 정보이다. SIFT는 영상 피라미드를 만들기 때문에 서로 다른 크기의 가우시안 마스크를 통해 계산을 수행한다. 그 후 영상의 경계선 정보에 해당하는 2차 미분값을 구하고 나면 해리스 코너 검출기와 비슷하게 비최댓값 억제 등의 기법을 사용하여 최종 결과에 해당하는 특징점들을 걸러낸다.

※ SIFT 특징 기술자 생성기

SIFT 특징 기술자에서는 SIFT 특징점 검출기를 사용하여 얻은 점들에 대해 검출된 스케일에 맞게 특징 기술자를 생성한다. 예를 들어, 특징점이 영상 피라미드의 1단계에서 검출되었다면 비교적 좁은 영역을 표현하는 특징 기술자를 얻고, 반대로 특징점이 영상 피라미드의 4단계에서 검출되었다면 더 큰 영여을 표현하는 특징 기술자를 얻는다. SIFT 특징 기술자도 기본은 영상의 미분 정보이고 미분 방향에 대한 히스토그램 정보로 구성된다. SIFT 특징 검출기로 얻은 특징점들은 해당 특징점이 검출된 피라키드 단계에 따른 스케이로가 미분으로 구한 방향 정보를 함께 가지고 있다. 이러한 정보를 다음 그림에서 사각형으로 표시하였다.

이 정보를 이용하여 입력 영상에서 특징점을 중심으로 16*16 픽셀 크기의 패치를 만든다. 이렇게 얻은 패치는 아래그림의 왼쪽과 같다. 특징 기술자는 이 패치를 총 16개의 4*4 픽셀 크기의 영역으로 분할하고 각각의 작은 영역에서 미분 방향 히스토그램을 구한다. 단, 히스토그램을 구할 때는 특징점의 위치를 중심으로 패치에 그린 내접원 안에 있는 픽셀만 다루도록 한다. 미분 방향 히스토그램을 구하려면 먼저 아래의 그림 가운데와 같이 각 픽셀의 미분방향을 구하고 이를 총 여덟 종류의 방향으로 분류하여 같은 방향끼리 미분값의 크기를 더한다. 그러면 아래 그림 오른쪽과 같이 4*4 픽셀 영역에 대한 미분값의 히스토그램으로 이루어진 특징 기술자를 얻는다.

 

이 특징 기술자는 여덟 개의 원소를 가지는 벡터로 표현할 수 이고 16*16픽셀에서 16개를 만들 수 있기 때문에 8*16=128개가 실제 SIFT 특정 기술자의 원소 갯수가 된다. 만약 두 특징점이 유사하다면 128차원의 특징 기술자는 비슷한 값ㅇ르 가지며, 특징점의 모습이 다르면 특징 기술자 값들이 큰 차이를 보인다.

SIFT 특징점 검출기는 최근에 개발된 알고리즘으로 구현하기도 어렵고 상업적으로 사용하기 어렵다. 따라서 직접 구현하기 보다는 원 개발자의 코드나 OpenCV 등에 구현된 라이브러리를 사용하는 것이 좋다. 여기서는 원 개발자가 개발한 프로그램(siftWin32.exe)을 사용한다. 해당 프로그램에 입력영상을 넘겨주면 다음과 같은 SIFT 특징 기술자 파일이 생성된다.

위 그림에서 첫 줄의 두 숫자는 특징점의 개수와 특징 기술자 벡터의 길이(차원)을 나타낸다. 총 1021개의 특징점을 얻었고 각 특징점은 총 128개의 정수 값으로 이루어진 특징 기술자 벡터를 가진다. 두 번째 줄부터는 각 특징점에 대한 정보들이 이어진다. 여기서 첫 줄은 특징점의 가로와 세로 좌표, 특징점의 스케일, 특징점의 방향 정보가 담겨 있고 이어서 각 특징점 기술자의 벡터정보가 기록되어 있다. 특징점 스케일은 첫 번째 그림에서 본 영역의 반지름을 나타내고, 특징점의 방향은 영역이 기울어진 방향을 나타낸다. 

이렇게 특징 기술자를 각 입력 영상에 대해 얻으면 원본 영상을 사용하지 않고 특징 기술자 값만을 사용하여 두 영상을 정합할 수 있다. 특징 기술자의 유사도는 두 특징 기술자 벡터 사이의 거리로 정의할 수 있고 이 거리가 짧을 수록 두 특징점이 유사하다는 의미가 된다. 

우리는 앞 포스트에서 패치 기반 정합 기법을 사용할 때 교차 검토를 통해 대응 관계의 유효성을 검토했다. 여기서는 또 다른 유효성 검사 기법인 두번 째로 가까운 이웃을 이용한 기법을 사용한다. 두 번째로 가까운 이웃이란 특징 기술자의 거리가 가장 가깝게 나온 특징점 바로 다음으로 가까운 거리를 가진 특징점을 가리킨다. 통계적으로 보았을 때 실제로 참인 대응 관계이면 특징 기술자의 거리가 두 번째로 가까운 거리부터는 큰 차이를 보인다는 점을 이용하여 찾아낸 대응 관계의 유효성을 판별하는 것이다. 

이를 다음 그림으로 표현했다. 참인 대응 관계를 가진 특징점이라면 다른 특징점과는 특징 기술자의 거리가 두드러지게 차이가 날 것이다. SIFT 특징 기술자를 이용한 정합 기법에서는 첫 번째와 두 번째로 작은 특징 기술자 거리의 비율이 0.6보다 작으면 유효한 대응인 것으로 판단한다.

※ SIFT 특징 기술자를 이용한 영상 정합

일단 특징 기술자가 기록되는 파일을 읽어오는 코드를 작성해야 한다. 따라서 특징점과 특징 기술자 정보를 저장할 수 있는 구조체를 다음과 같이 정의한다.

struct Keypoint{
	float			x, y;					// SIFT 특징점의 x, y 좌표
	float			scale, ori;				// SIFT 특징점의스케일과 방향
	unsigned char	descrip[DIM_SIFT+1];		// 특징 기술자 배열
};

#define DIM_SIFT 128   

특징 기술자의 각 원소는 0에서 255사이의 정수값이기 때문에 unsigned char형으로 특징 기술자의 배열은 선언한다. 또한 SIFT 특징 기술자는 128개의 원소를 가지고 있으므로 #define으로 128을 선언했다. 특징 기술자 배열에 1을 더해주는 이유는 변수를 문자열로 인식할 수 있어 NULL 문자를 저장할 수 있게 하기 위해서이다. 파일에 담긴 한 영상의 특징 기술자들은 Keypoint구조체의 배열에 저장되는데 다음 메서드를 통해 배열을 채울 수 있다.

bool ReadDescriptors(FILE *pFile, Keypoint* pKeyPt, int nDesc)
{
	for (int i=0 ; i<nDesc ; i++)
	{
		Keypoint& pt = pKeyPt[i];
		// 특징점의 좌표, 스케일, 방향 정보를 읽응ㅁ
		if (fscanf_s(pFile, "%f %f %f %f", &pt.y, &pt.x, &pt.scale, &pt.ori) != 4)
			return false;

		// 특징 기술자 배열을 읽음
		for (int j=0 ; j<DIM_SIFT ; j++)
		{
			//fscanf_s는 읽은 원소의 갯수 반환
			if (fscanf_s(pFile, "%d", &pt.descrip[j]) != 1 || 
				pt.descrip[j] < 0 || pt.descrip[j] > 255)
				return false;
		}
	}
    return true;
}

다음은 특징 기술자를 비교하는 함수이다.

int _CalcSIFTSqDist(const Keypoint& k1, const Keypoint& k2)
{
	const unsigned char *pk1, *pk2;

	pk1 = k1.descrip;
	pk2 = k2.descrip;

	int dif, distsq = 0;
	for (int i=0 ; i<DIM_SIFT ; i++) 
	{
		dif = (int)pk1[i] - (int)pk2[i];
		distsq += dif * dif;
	}
	return distsq;
}

마지막으로 위에서 구현한 코드를 바탕으로 SIFT 기반 영상 정합을 수행하는 함수이다. 

void CImageMatchDlg::OnBnClickedButtonMatchSift()
{
	// TODO: 여기에 컨트롤 알림 처리기 코드를 추가합니다.
	system("siftWin32.exe <scene.pgm >scene.txt");

	FILE *pFile1 = NULL, *pFile2 = NULL; 

	char szFilter[] = "Descriptor File (*.TXT) | *.TXT; | All Files(*.*)|*.*||";	
	CFileDialog dlg(TRUE, NULL, NULL, OFN_HIDEREADONLY, szFilter);

	if (IDOK == dlg.DoModal()) // 파일 대화 상자 열기
	{
		CString strPathName = dlg.GetPathName();
		fopen_s(&pFile1, strPathName, "r");
	}
	if (IDOK == dlg.DoModal()) // 파일 대화 상자 열기
	{
		CString strPathName = dlg.GetPathName();
		fopen_s(&pFile2, strPathName, "r");
	}
	if (!pFile1 || !pFile2)
		return;

	//각 영상의 특징 기술자의 수를 읽음
	int nDesc1, nDesc2, nDim1, nDim2;
	if (fscanf_s(pFile1, "%d %d", &nDesc1, &nDim1) != 2)
		return;
	if (fscanf_s(pFile2, "%d %d", &nDesc2, &nDim2) != 2)
		return;

	ASSERT(nDim1==DIM_SIFT && nDim2==DIM_SIFT);

	//결과 출력용 영상 할당
	m_imageOut = CByteImage(m_imageIn1.GetWidth() + m_imageIn2.GetWidth(), 
						MAX(m_imageIn1.GetHeight(), m_imageIn2.GetHeight()), 3);
	int nWidth1 = m_imageIn1.GetWidth();
	m_imageOut.Paste(Gray2RGB(m_imageIn1), 0, 0);
	m_imageOut.Paste(Gray2RGB(m_imageIn2), nWidth1, 0);

	//특징 기술자 배열을 할당하고 배열에 값 저장
	Keypoint* pKeyPt1 = new Keypoint[nDesc1];
	Keypoint* pKeyPt2 = new Keypoint[nDesc2];

	ReadDescriptors(pFile1, pKeyPt1, nDesc1);
	ReadDescriptors(pFile2, pKeyPt2, nDesc2);

	for (int i=0 ; i<nDesc1 ; i++) 
	{
		int nIdx2 = -1;
		int minDist1 = INT_MAX, minDist2 = INT_MAX;

		for (int j=0 ; j<nDesc2 ; j++) 
		{
			int distSq = _CalcSIFTSqDist(pKeyPt1[i], pKeyPt2[j]);

			if (distSq < minDist1) // 가장 가까운 특징점 갱신
			{
				minDist2 = minDist1;
				minDist1 = distSq;
				nIdx2 = j;
			}
			else if (distSq < minDist2) // 두번째로 가까운 특징점 갱신
			{
				minDist2 = distSq;
			}
		}

		if (10*10*minDist1 < 6*6*minDist2) // 유효한 대응 관계 판단
		{
			DrawLine(m_imageOut, pKeyPt1[i].x, pKeyPt1[i].y, 
				pKeyPt2[nIdx2].x+nWidth1, pKeyPt2[nIdx2].y, 255, 0, 0);
		}
	}

	ShowImage(m_imageOut, "정합 결과");

	delete [] pKeyPt1;
	delete [] pKeyPt2;
}

'영상처리 프로그래밍' 카테고리의 다른 글

주파수 스펙트럼과 위상각의 시각화  (1)	2017.06.22
고속 푸리에 변환/FFT(2)  (1)	2017.06.22
NCC, 교차 검토, 특징 기술자  (0)	2017.06.19
영상 정합  (0)	2017.06.19
소벨/캐니 경계(Edge), 허프 변환 직선, 해리스 코너  (0)	2017.06.18