[C++/OpenCV] 공간 주파수, 이산 푸리에 변환, DFT

 

 

영상을 데이터로 표현하는 데에는 크게 두 가지 영역으로 나누어 설명한다. 화소값이 직접 표현된 공간 영역과 우주 공간과 같은 변환 영역이다.

 

변환 영역은 직교 변환에 의해 얻어진 영상 데이터의 다른 표현이다. 여기서는 화소값이 직접 표현되는 것이 아니고 변환 계수(coefficient)로 표현된다. 대표적인 변환은 DCT(Discrete Cosine Transform)와 DFT(Discrete Fourier Transform)이 있다. 

 

공간 주파수란?

영상 신호에서의 주파수는 공간상에서 화소 밝기의 변화율이라 할 수 있다. 이런 의미에서 공간 주파수라는 표현을 사용한다.

 

공간 주파수는 밝기가 얼마나 빨리 변화하는가에 다라 고주파 영역과 저주파 영역으로 분류한다.

- 고주파 영역: 화소 밝기가 급변하는 영역 / 배경 or 객체 내부
- 저주파 영역: 화소 밝기의 변화가 거의 없거나 점진적으로 변화하는 영역 / 엣지(Edge)

 

주파수 영역 변환의 필요성?

변환을 통하여 영상을 주파수 영역별로 분리할 수 있으면, 각 주파수 영역별 처리가 가능하다. 엣지 부분을 구성하는 고주파 성분을 제거하면 경계가 흐려진 영상을 만들 수 있고, 저주파 성분을 제거한다면 에지 영상을 만들 수 있을 것이다.

 

그러나 현실 영상에서는 저주파와 고주파 성분이 혼합되어 존재하기 때문에 주파수 영역에서의 영상 처리가 필요하다. 주파수 변환으로 얻어진 계수의 특정 주파수 영역에 원하는 영상 처리를 적용한 뒤 공간영역 영상으로 변환하면 원하는 주파수 성분만 존재하는 영상을 만들 수 있다.

 

이산 푸리에 변환(Discrete Fourier Transform, DFT)

푸리에 변환은 "주기를 가진 신호는 정현파/여현파의 합으로 표현할 수 있다"는 전제를 기본으로 한다. 따라서 사인(sin) 또는 코사인(cos) 함수의 선형 조합으로 특정 주기의 신호를 구성할 수 있다. 여기서 조합되는 각각의 사인(sin) 또는 코사인(cos) 함수를 기저 함수(basis function)이라 부르며, 기저 함수에 곱해지는 값들은 주파수의 계수(coefficient)가 된다. 이 계수가 각 주파수 성분의 크기에 해당하며, 신호를 주파수로 변환하는 것은 각 주파수의 기저 함수들에 대한 계수를 찾는 것이다.

 

다음은 이산 푸리에 변환(Discrete Fourier Transform, DFT)(위)과 그 역변환(아래)의 수식이다. 

 

DFT(위)과 그 역변환(아래)

 

$g[n]$은 디지털 신호, $G(k)$는 주파수 $k$에 대한 푸리에 변환 계수이다. 또한 이산 신호이기 때문에 $k$와 $n$은 신호의 원소 개수($N$)만큼 정수로 주어진다. 

 

2차원 공간상의 영역에서 DFT를 수행하려면 1차원 DFT를 가로 방향과 세로 방향으로 연속해서 두 번 적용해야 한다. 다음 수식은 2차원 DFT에 대한 수식이다. 괄호부분은 가로 방향에 대한 1차원 푸리에 변환이다.

 

2차원 공간상의 영역에서 DFT

 

아래는 2차원 DFT 역변환에 대한 수식이다.

 

2차원 공간상의 영역에서 DFT 역변환

 

그러나 이게 끝이 아니다. 푸리에 변환의 기저 함수에 허수를 이용했기 때문에 실수부와 함께 허수부에 대한 고려도 해야 한다. 다음은 기저함수를 사인과 코사인으로 변경하고 실수부와 허수부를 구분하여 1차원 DFT을 나타낸 것이다. 푸리에 변환과 그 역변환이 사인과 코사인 함수에서 각도의 부호만 반대이다.

 

기저함수를 사인과 코사인으로 변경하고 실수부와 허수부를 구분한 1차원 DFT

 

푸리에 변환을 수행하면 복소수의 행렬이 결과로 생성된다. 이것을 영상으로 확인하기 위해서는 복소수의 실수부와 허수부를 벡터로 간주하여 벡터의 크기를 구하면 된다. 이것을 주파수 스펙트럼이라 한다. 또한 실수부와 허수부의 각도를 이용해서 주파수 위상을 계산할 수도 있다.

 

소스코드 (C++)

다음은 DFT를 C++로 짠 코드이다.

#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;

void log_mag(Mat complex, Mat& dst) { //푸리에 변환을 수행하면 복소수의 행렬이 결과로 생성
	Mat planes[2]; //복소수의 실수부와 허수부를 벡터로 간주
	split(complex, planes); //2채널 행렬 분리
	magnitude(planes[0], planes[1], dst); //벡터의 크기
	log(dst + 1, dst);
	normalize(dst, dst, 0, 255, NORM_MINMAX); //정규화(저주파 영역과 고주파 영역의 계수값을 정규화)
	dst.convertTo(dst, CV_8U);
}

void shuffling(Mat mag_img, Mat& dst) {
	int cx = mag_img.cols / 2;
	int cy = mag_img.rows / 2;
	Rect q1(cx, 0, cx, cy); //1사분면
	Rect q2(0, 0, cx, cy); //2
	Rect q3(0, cy, cx, cy); //3
	Rect q4(cx, cy, cx, cy); //4

	dst = Mat(mag_img.size(), mag_img.type());
	mag_img(q1).copyTo(dst(q3)); //사분면 맞바꿈
	mag_img(q3).copyTo(dst(q1));
	mag_img(q2).copyTo(dst(q4));
	mag_img(q4).copyTo(dst(q2));
}

Mat DFT_1D(Mat row, int dir) { //1차원 신호의 이산 푸리에 변환
	int n = row.cols;
	Mat dst(row.size(), CV_32FC2);

	for (int i = 0; i < n; i++) {
		Vec2f complex(0, 0);
		for (int j = 0; j < n; j++) {
			float theta = dir * (-2) * CV_PI * i * j / n; //기저함수 각도 계산
			Vec2f value = row.at<Vec2f>(j);
			complex[0] += value[0] * cos(theta) - value[1] * sin(theta);
			complex[1] += value[1] * cos(theta) + value[0] * sin(theta);
		}
		dst.at<Vec2f>(i) = complex;  //한 원소의 DFT 계산 결과 저장
	}

	if (dir == -1) { //-1이면 역변환
		dst /= n;
	}

	return dst;
}

void DFT_2D(Mat complex, Mat& dst, int dir) { //2차원 신호 푸리에 변환
	complex.convertTo(complex, CV_32F);
	Mat tmp(complex.size(), CV_32FC2, Vec2f(0, 0));
	tmp.copyTo(dst);

	for (int i = 0; i < complex.rows; i++){  //가로 방향 푸리에 변환
		Mat one_row = complex.row(i);
		Mat dft_row = DFT_1D(one_row, dir);  //1개 행 변환
		dft_row.copyTo(tmp.row(i)); //저장
	}

	transpose(tmp, tmp); //전치
	for (int i = 0; i < tmp.rows; i++) {  //세로 방향 푸리에 변환
		Mat one_row = tmp.row(i);
		Mat dft_row = DFT_1D(tmp.row(i), dir); 
		dft_row.copyTo(dst.row(i));
	}
	transpose(dst, dst);
}

int main() {
	Mat image = imread("Lenna.png", 0);
	Mat complex, dft_coef, dft_img, idft_coef, shuffle, idft_img[2];  
	Mat tmp[] = { image, Mat::zeros(image.size(), CV_8U) };
	merge(tmp, 2, complex); //복소수 행렬 구성

	DFT_2D(complex, dft_coef, 1); //2차원 DFT
	log_mag(dft_coef, dft_img);
	shuffling(dft_img, shuffle);

	DFT_2D(dft_coef, idft_coef, -1); // 푸리에 역변환으로 원본 영상 복원
	split(idft_coef, idft_img);
	idft_img[0].convertTo(idft_img[0], CV_8U);

	imshow("image", image);
	imshow("dft_img", dft_img);
	imshow("shuffling", shuffle);
	imshow("idft_img", idft_img[0]);
	waitKey();
	return 0;
}

 

결과

- 입력 영상

 

 

- DFT 수행 후 스펙트럼 영상

  : 영상의 중심이 저주파 영역, 바깥쪽이 고주파 영역

 

 

- DFT 역변환으로 다시 영상 복원

 

 

 

 

 

 

 

* 이 글은 <OpenCV로 배우는 영상 처리 및 응용> (정성환, 생능출판) 책으로 공부하며 작성한 것이다.