주파수 영역 필터링(Filtering in Frequency Domain)

영상 신호의 주파수 영역에서 분석하여 공간 영역에서 잘 보이지 않던 정보를 쉽게 발견할 수 있다. 

Fourier Transform은 time domain에서의 영상을 frequency domain으로의 성분으로 변환하는 것이다. 

즉 임의의 입력 신호를 다양한 주파수를 갖는 주기함수들의 합으로 분해하여 표현하는 것이다. 

이 푸리에 변환으로 F(w) 스펙트럼을 얻을 수 있고 w는 실수부와 허수부로 구성된다. 

 

이를 Specturm 또는 magnitude라고 부르고, F(w)의 각도 Φ를 phase angle(또는 phase spectrum)라고 부른다 .

 

 

Discrete Fourier Transform

2D DFT

Mat doDft(Mat src_img) {
	Mat float_img;
	src_img.convertTo(float_img, CV_32F);

	Mat complex_img;
	dft(float_img, complex_img, DFT_COMPLEX_OUTPUT);

	return complex_img;
}

CV_32F(32-bit floating-point number) 를 사용해서 Mat의 고정 픽셀 유형을 지정한다. 

들어온 이미지를 주파수 도메인으로 변환한다. 

 

주파수 도메인으로 변환된 이미지를 기반으로

Magnitude 영상을 취득할 수 있는데, 

Mat getMagnitude(Mat complexImg) {
	Mat planes[2];
	split(complexImg, planes);
	// 실수부, 허수부 분리
	
	Mat magImg; 
	magnitude(planes[0], planes[1], magImg);
	magImg += Scalar::all(1);
	log(magImg, magImg);// magitude 취득

	return magImg;
}

opencv의 내장함수 magnitude를 이용해서 complexImg의 크기를 구한다. 

magnitude 함수는 2D벡터의 x좌표 행렬, y좌표 행렬과 2D벡터의 크기 행렬을 구해

결과를 담을 행렬을 인자로 받아 벡터의 크기를 계산한다. 

Mat myNormalize(Mat src) {
	Mat dst;
	src.copyTo(dst);
	normalize(dst, dst, 0, 255, NORM_MINMAX);
	dst.convertTo(dst, CV_8UC1);

	return dst;
}

정규화 과정을 거쳐서, 영상이 0에서 255의 범위를 갖도록 한다. 

Mat getPhase(Mat complexImg) {
	Mat planes[2];
	split(complexImg, planes);

	Mat phaImg;
	phase(planes[0], planes[1], phaImg);
	return phaImg;
}

 

벡터의 방향을 계산하는 phase함수를 사용해 벡터의 방향 행렬을 담는다. 

 

 

좌표계 중앙 이동 (centralize) 

Mat centralize(Mat complex) {
	Mat planes[2];
	split(complex, planes);
	int cx = planes[0].cols / 2;
	int cy = planes[1].rows / 2;

	Mat q0Re(planes[0], Rect(0, 0, cx, cy));
	Mat q1Re(planes[0], Rect(cx, 0, cx, cy));
	Mat q2Re(planes[0], Rect(0, cy, cx, cy));
	Mat q3Re(planes[0], Rect(cx, cy, cx, cy));

	Mat tmp;
	q0Re.copyTo(tmp);
	q3Re.copyTo(q0Re);
	tmp.copyTo(q3Re);
	q1Re.copyTo(tmp);
	q2Re.copyTo(q1Re);
	tmp.copyTo(q2Re);

	Mat q0Im(planes[1], Rect(0, 0, cx, cy));
	Mat q1Im(planes[1], Rect(cx, 0, cx, cy));
	Mat q2Im(planes[1], Rect(0, cy, cx, cy));
	Mat q3Im(planes[1], Rect(cx, cy, cx, cy));

	q0Im.copyTo(tmp);
	q3Im.copyTo(q0Im);
	tmp.copyTo(q3Im);
	q1Im.copyTo(tmp);
	q2Im.copyTo(q1Im);
	tmp.copyTo(q2Im);

	Mat centerComplex;
	merge(planes, 2, centerComplex);
	return centerComplex;
}

왼쪽 상단으로 원점이 맞춰져 있는 영상의 초점을 중앙으로 이동한다. 

영상을 4등분하여 조합해서 재배치한 방식이다. 

 

Discrete Fourier Transform

Mat setComplex(Mat mag_img, Mat pha_img) {
	exp(mag_img, mag_img);
	mag_img -= Scalar::all(1);
	// magnitude 계산을 반대로 수행

	Mat planes[2];
	polarToCart(mag_img, pha_img, planes[0], planes[1]);
	// 극 좌표계 -> 직교 좌표계 (각도와 크기로부터 2차원 좌표)

	Mat complex_img;
	merge(planes, 2, complex_img);
	// 실수부, 허수부 합체

	return complex_img;
}

변환된 magnitude 이미지와 phase 이미지를 토대로 원본 ComplexImg를 가져오는 setComplex 함수이다. 

Mat doIdft(Mat complex_img) {
	Mat idftcvt;
	idft(complex_img, idftcvt);
	// IDFT를 이용한 원본 영상 취득

	Mat planes[2];
	split(idftcvt, planes);

	Mat dst_img;
	magnitude(planes[0], planes[1], dst_img);
	normalize(dst_img, dst_img, 255, 0, NORM_MINMAX);
	dst_img.convertTo(dst_img, CV_8UC1);
	// 일반 영상의 type과 표현범위로 변환

	return dst_img;
}

idft는 dft, discrete fourier transform(이산 푸리에 변환)의 역변환이다.

즉 영방향 이산 푸리에 변환을 수행하는 것인데, 

샘플링한 주파수 형태를 역변환하여 

 

 

Filtering

 

Low pass filtering(LPF)는 주파수 성분에서 고주파 성분을 배제한다. 

영상의 에지 부분이 흐려지게 되어 영상이 전반적으로 스무딩하게 된다. 

// Low Pass Filtering(주파수 성분에서 고주파 성분을 배제)
Mat doLPF(Mat src_img) {
	Mat pad_img = padding(src_img);
	Mat complex_img = doDft(pad_img);
	Mat center_complex_img = centralize(complex_img);
	Mat mag_img = getMagnitude(center_complex_img);
	Mat pha_img = getPhase(center_complex_img);

	double minVal, maxVal;
	Point minLoc, maxLoc;
	minMaxLoc(mag_img, &minVal, &maxVal, &minLoc, &maxLoc);
	normalize(mag_img, mag_img, 0, 1, NORM_MINMAX);

	Mat mask_img = Mat::zeros(mag_img.size(), CV_32F);
	circle(mask_img, Point(mask_img.cols / 2, mask_img.rows / 2), 20, Scalar::all(1), -1, -1, 0);

	Mat mag_img2;
	multiply(mag_img, mask_img, mag_img2);

	normalize(mag_img2, mag_img2, (float)minVal, (float)maxVal, NORM_MINMAX);
	Mat complex_img2 = setComplex(mag_img2, pha_img);
	Mat dst_img = doIdft(complex_img2);

	return myNormalize(dst_img);
}

 

High pass filtering(HPF) 

 

LPF와는 반대로 주파수 성분에서 저주파를 배제하고 고주파 영역만 남기므로

에지 부분이 선명하게 남게 된다.  (에지 검출)

 

Mat doHPF(Mat srcImg) {

	//DFT
	Mat padImg = padding(srcImg);
	Mat complexImg = doDft(padImg);
	Mat centerComplexImg = centralize(complexImg);
	Mat magImg = getMagnitude(centerComplexImg);
	Mat phaImg = getPhase(centerComplexImg);

	//HPF
	double minVal, maxVal;
	Point minLoc, maxLoc;

	minMaxLoc(magImg, &minVal, &maxVal, &minLoc, &maxLoc);
	normalize(magImg, magImg, 0, 1, NORM_MINMAX);

	Mat maskImg = Mat::ones(magImg.size(), CV_32F);
	circle(maskImg, Point(maskImg.cols / 2, maskImg.rows / 2), 50, Scalar::all(0), -1, -1, 0);

	Mat magImg2;
	multiply(magImg, maskImg, magImg2);

	//imshow("HW1", magImg2);
	//waitKey(0);
	//destroyWindow("HW1");

	//IDFT
	normalize(magImg2, magImg2, (float)minVal, (float)maxVal, NORM_MINMAX);
	Mat complexImg2 = setComplex(magImg2, phaImg);
	Mat dstImg = doIdft(complexImg2);

	return myNormalize(dstImg);
}

 

위 예제들의 main함수는 다음과 같다!

void ex1() {
	Mat srcImg = imread("gear.jpg", 0);
	Mat dstImg = doDft(srcImg);

	dstImg = getMagnitude(dstImg);
	dstImg = myNormalize(dstImg);

	imshow("srcImg", srcImg);
	imshow("dstImg", dstImg);
	waitKey(0);
	destroyAllWindows();
}

void ex2() {
	Mat srcImg = imread("gear.jpg", 0);
	Mat dstImg = doDft(srcImg);
	dstImg = centralize(dstImg);
	dstImg = getPhase(dstImg);
	dstImg = myNormalize(dstImg);

	imshow("srcImg", srcImg);
	imshow("dstImg", dstImg);
	waitKey(0);
	destroyAllWindows();
}

void ex3() {
	Mat srcImg = imread("gear.jpg", 0);
	Mat dstImg = doLPF(srcImg);

	imshow("srcImg", srcImg);
	imshow("dstImg", dstImg);
	waitKey(0);
	destroyAllWindows();
}

void ex4() {
	Mat srcImg = imread("gear.jpg", 0);
	Mat dstImg = doHPF(srcImg);

	imshow("srcImg", srcImg);
	imshow("dstImg", dstImg);
	waitKey(0);
	destroyAllWindows();
}

Mat padding(Mat img) {
	int dftRows = getOptimalDFTSize(img.rows);
	int dftCols = getOptimalDFTSize(img.cols);

	Mat padded;
	copyMakeBorder(img, padded, 0, dftRows - img.rows, 0, dftCols - img.cols, BORDER_CONSTANT, Scalar::all(0));

	return padded;
}