DIP 실습 - Band pass filter / Frequency domain

---

과제

img1.jpg

img2.jpg

img3.jpg

---

풀이 코드

#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
using namespace cv;
using namespace std;

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;
}

Mat doDft(Mat src_img) { // freqeuncy domain으로 transform
	Mat float_img;
	src_img.convertTo(float_img, CV_32F);

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

	return complex_img;
}

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);
	// manitude 취득
	//log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)

	return magImg;
}

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

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

	return phaImg;
}

Mat myNormalize(Mat src) { // 정규화과정 (0~255)
	Mat dst;
	src.copyTo(dst);
	normalize(dst, dst, 0, 255, NORM_MINMAX);
	dst.convertTo(dst, CV_8UC1);

	return dst;
}

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;
}

Mat setComplex(Mat mag_img, Mat pha_img) {
	exp(mag_img, mag_img); // 앞에서 로그 취하고 1 더해줘서, 반대 과정을 해줌
	mag_img -= Scalar::all(1);

	Mat planes[2];
	polarToCart(mag_img, pha_img, planes[0], planes[1]);

	Mat complex_img;
	merge(planes, 2, complex_img);

	return complex_img;
}

Mat doIdft(Mat complex_img) { // 주파수 영역을 공간 영역으로 바꿈
	Mat idftcvt;
	idft(complex_img, idftcvt);

	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);

	return dst_img;
}

Mat doBPF(Mat srcImg) {
	//DFT
	Mat padImg = padding(srcImg);
	Mat complexImg = doDft(padImg);
	Mat centerComplexImg = centralize(complexImg);
	Mat magImg = getMagnitude(centerComplexImg);
	Mat phaImg = getPhase(centerComplexImg);
	Mat mag_img2, bandmask;

	double minVal, maxVal;
	Point minLoc, maxLoc;

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

	// BPF 마스크
	Mat maskImg_out = Mat::zeros(magImg.size(), CV_32F); // zero matrix 생성
	circle(maskImg_out, Point(maskImg_out.cols / 2, maskImg_out.rows / 2), 60, Scalar::all(1), -1, -1, 0);

	Mat maskImg_in = Mat::zeros(magImg.size(), CV_32F); // one matrix 생성
	circle(maskImg_in, Point(maskImg_in.cols / 2, maskImg_in.rows / 2), 20, Scalar::all(1), -1, -1, 0);

	// 내부와 외부 마스크를 빼서 밴드 패스 마스크 생성
	subtract(maskImg_out, maskImg_in, bandmask);
	multiply(bandmask, magImg, mag_img2);

	// bandmask와 mag_img2를 수평으로 결합하여 하나의 이미지로 출력
	Mat combinedImg;
	std::vector<cv::Mat> imgs = { bandmask, mag_img2 };
	hconcat(imgs, combinedImg);
	// 병합된 이미지를 보여주는 창 생성
	imshow("Band-Pass Filter Visualization", combinedImg);
	waitKey(0);

	// 마스크 적용 후 주파수 이미지를 정규화
	normalize(mag_img2, mag_img2, (float)minVal, (float)maxVal, NORM_MINMAX);
	Mat complexImg2 = setComplex(mag_img2, phaImg);
	Mat dstImg = doIdft(complexImg2);

	return myNormalize(dstImg);
}

int myKernelConv3x3(uchar* arr, int kernel[][3], int x, int y, int width, int height) {
	int sum = 0;
	int sumKernel = 0;

	for (int j = -1; j <= 1; j++) {
		for (int i = -1; i <= 1; i++) {
			if ((y + j) >= 0 && (y + j) < height && (x + i) >= 0 && (x + i) < width) {
				sum += arr[(y + j) * width + (x + i)] * kernel[i + 1][j + 1];
				sumKernel += kernel[i + 1][j + 1];
			}
		}
	}

	if (sumKernel != 0) { return sum / sumKernel; }
	else return sum;
}
Mat mysobelFilter(Mat srcImg) {
	int kernelY[3][3] = { -1, -2, -1, 
							0, 0, 0, 
							1, 2, 1 };  // 수직 에지 감지 커널

	Mat dstImg(srcImg.size(), CV_8UC1);
	uchar* srcData = srcImg.data;
	uchar* dstData = dstImg.data;
	int width = srcImg.cols;
	int height = srcImg.rows;

	for (int y = 0; y < height; y++) {
		for (int x = 0; x < width; x++) {
			int edgeStrength = abs(myKernelConv3x3(srcData, kernelY, x, y, width, height));
			dstData[y * width + x] = edgeStrength;
		}
	}

	return dstImg;
}
Mat freqSobel(Mat srcImg) {
	// DFT
	Mat padImg = padding(srcImg);
	Mat complexImg = doDft(padImg);
	Mat centerComplexImg = centralize(complexImg);
	Mat magImg = getMagnitude(centerComplexImg);
	Mat phaImg = getPhase(centerComplexImg);

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

	// 삼각형 필터 생성, 저주파 영역 제외
	Mat maskImg = Mat::zeros(magImg.size(), CV_32F);
	int offset = 70;  // 중심에서 떨어진 거리
	int thickness = 50;  // 선의 두께
	// 상단 삼각형
	line(maskImg, Point(maskImg.cols / 2 + offset, 0), Point(maskImg.cols / 2 + 50, maskImg.rows / 2 - offset), Scalar::all(1), thickness);
	line(maskImg, Point(maskImg.cols / 2 - offset, 0), Point(maskImg.cols / 2 - 50, maskImg.rows / 2 - offset), Scalar::all(1), thickness);
	// 하단 삼각형
	line(maskImg, Point(maskImg.cols / 2 + offset, maskImg.rows), Point(maskImg.cols / 2 + 50, maskImg.rows / 2 + offset), Scalar::all(1), thickness);
	line(maskImg, Point(maskImg.cols / 2 - offset, maskImg.rows), Point(maskImg.cols / 2 - 50, maskImg.rows / 2 + offset), Scalar::all(1), thickness);
	line(maskImg, Point(maskImg.cols / 2, 0), Point(maskImg.cols / 2, maskImg.rows / 2 - offset), Scalar::all(1), 70);
	line(maskImg, Point(maskImg.cols / 2, maskImg.rows / 2 + offset), Point(maskImg.cols / 2, maskImg.rows), Scalar::all(1), 70);


	Mat mag_img2;
	imshow("freqeuncy domain", magImg);
	multiply(maskImg, magImg, mag_img2);

	imshow("band", maskImg); // bandmask 출력
	imshow("band_img", mag_img2); // bandmask가 곱해진 이미지 출력
	waitKey(0);

	// IDFT
	normalize(mag_img2, mag_img2, (float)minVal, (float)maxVal, NORM_MINMAX);
	Mat complexImg2 = setComplex(mag_img2, phaImg);
	Mat dstImg = doIdft(complexImg2);
	return myNormalize(dstImg);
}


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

	double minVal, maxVal;
	Point minLoc, maxLoc;

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

	// delete Filckering
	Mat maskImg = Mat::ones(magImg.size(), CV_32F); // 이미지 크기만큼 1행렬 생성

	line(maskImg, Point(maskImg.cols / 2, 23 * maskImg.rows / 30), Point(maskImg.cols / 2, 16 * maskImg.rows / 30), (0, 0, 0), 35);
	line(maskImg, Point(maskImg.cols / 2, 7 * maskImg.rows / 30), Point(maskImg.cols / 2, 14 * maskImg.rows / 30), (0, 0, 0), 35);

	multiply(maskImg, magImg, mag_img2);

	Mat combined;
	hconcat(maskImg, mag_img2, combined);
	imshow("Filter & Filtered Image", combined);
	waitKey(0);

	//IDFT
	normalize(mag_img2, mag_img2, (float)minVal, (float)maxVal, NORM_MINMAX);
	Mat complexImg2 = setComplex(mag_img2, phaImg);
	Mat dstImg = doIdft(complexImg2);
	return myNormalize(dstImg);
}

void taskBPF() {
	Mat src_img = imread("img1.jpg", 0);
	Mat dst_img = doBPF(src_img);

	Mat combinedImg;
	hconcat(src_img, dst_img, combinedImg);

	imshow("BPF - Original vs Filtered", combinedImg);
	waitKey(0);
	destroyAllWindows();
}

void taskSobel() {
	Mat src_img = imread("img2.jpg", 0);
	Mat dst_img = mysobelFilter(src_img);
	Mat dst_img_2 = freqSobel(src_img);

	imshow("src_img", src_img);
	imshow("Spatial Domain Sobel", dst_img);
	imshow("Frequency Domain Sobel", dst_img_2);
	waitKey(0);
	destroyAllWindows();
}

void taskFlickering() {
	Mat src_img = imread("img3.jpg", 0);
	Mat dst_img = solveFlickering(src_img);

	imshow("Before", src_img);
	imshow("After", dst_img);
	waitKey(0);
	destroyAllWindows();
}

int main() {
	taskBPF();
	taskSobel();
	taskFlickering();

	return 0;
}

doBPF함수에서 Band Pass Filter는 특정 주파수 범위의 신호를 통과시키고 다른 신호는 차단하는 주파수 도메인 마스크를 생성하도록 하였다. `maskImg_out`과 `maskImg_in` 두 마스크를 초기화하여 각각 외부 반지름과 내부 반지름을 각각 20, 60으로 정의하고, 이 두 마스크를 서로 빼서 도넛 모양의 밴드 패스 마스크를 만들었다. 이 마스크는 원본 이미지의 크기에 맞춰서 정해진 비율로 조절되며, 주파수 영역에서 이미지의 특정 세부사항을 강조하거나 제거하는 데 사용될 수 있다.

 

이전 과제들에서 만들어서 사용했던 myKernelConv3x3과 수평방향 Sobel을 구현한 mysobelFilter를 구현해뒀다. 이는 spatial domain에서 filtering을 하여 frequency domain에서의 filtering과 비교하기 위한 용도이다.

`freqSobel` 함수에서는 주파수 영역에서 소벨 필터와 유사한 기능을 수행하여 이미지의 경계 감지를 위해 설계된 필터를 구현했다. 이 함수는 원본 이미지에 대해 이산 푸리에 변환(DFT)을 수행하고, 변환된 결과를 정규화한 후 직접 설계한 삼각형 모양의 마스크를 적용했다. 이 마스크는 위아래로 꼭짓점 하나가 중심을 향한 삼각형을 하나씩 배치한 필터 모양에서 저주파부를 뜻하는 중심부를 제거한 모양의 마스크이다. 주파수 영역에서 수직방향의 고주파 성분을 강조하여 이미지의 수평 edge를 감지할 수 있다. 최종적으로, 마스크가 적용된 주파수 영역 이미지를 역 푸리에 변환하여 시각적으로 확인할 수 있는 공간 영역의 이미지로 변환하여 수평 edge를 강조한 image를 시각적으로 확인할 수 있도록 하였다.

 

`solveFlickering` 함수는 비디오나 이미지에서 발생하는 flickering 현상을 제거하기 위해 설계되었습니다. 이 함수는 주파수 영역에서 flickering에 해당하는 수직의 중,고 주파 대역을 차단하는 마스크를 적용했습니다. 정확히 얼만큼 어느 부분을 제거해야할지 몰라서 직접 테스트를 해가며 구체적인 수치를 찾아야했습니다. 그 결과 마스크는 이미지의 중앙에서 수직으로 두 개의 검은 선을 그려서, 각각의 선이 이미지 높이의 23/30부분에서 16/30부분, 14/30부분에서 7/30부분 까지의 범위를 차단하도록 설계하는 것이 flickering 제거 성능이 준수한 것을 확인하였습니다. 이 선들의 두께는 35 픽셀로 설정하였습니다. 이렇게 생성된 마스크는 원본 이미지의 주파수 대표 영역에서 flickering을 유발할 수 있는 주파수를 효과적으로 제거함으로써 보다 안정적인 이미지를 생성하였습니다. 최종적으로 마스크가 적용된 주파수 영역을 역 푸리에 변환하여 공간 영역으로 복구하고, 결과 이미지를 정규화하여 출력하였습니다.

 

---

결과

1번 필터

BPS역할을 하는 도넛모양의 필터이고 해당 필터를 주파수영역으로 변환된 이미지에 씌운 모습이다.

 

1번 결과

원본 영상에 Band Pass Filter를 적용한 영상이다. Band Pass Filter이므로 설정한 bandmask에 해당하는 주파수 성분만 남아 출력된다.

 

2번 필터

정확하게 삼각형모양을 만들지는 못했지만 주파수 영역에서 수직방향의 밝게 강조된 부분들은 모두 포함하고 저주파 영역은 걸러내는 필터가 구현된 것을 확인 할 수 있다. 이를 주파수영역으로 변환된 영상에 씌운 결과도 확인할 수 있다.

 

2번 결과 (원본 vs Spatial Domain Sobel vs Frequency Domain Sobel)

Spatial Domain의 Sobel filter와 보이는 것은 조금 다르지만 원하던 바인 수평방향의 edge를 살리는 것은 Frequency Domain Sobel도 잘 수행한 것을 볼 수 있다. 옷의 전체적인 테두리는 spatial domain Sobel이 더 명확하게 포착했으나 오히려 몸통 부분에서는 Frequency Domain Sobel이 좀더 나은 결과를 만들어냈다.

 

3번 필터

Flickering은 수평방향으로 생긴 검은 줄무늬 이므로 수평방향의 주파수 내용을 제거하는 필터를 만들었다. 직접 다양한 테스트를 해본 결과 저주파부는 자연스러운 내용이므로 살려야했고 고주파부를 아예 제거하니 원본이 사라지는 경향이 보여서 최소한의 저주파부를 남기고 고주파부도 어느정도 남긴 세로 방향의 중고주파를 제거하니 좋은 결과가 나왔다.

 

3번 결과 (원본 vs flickering 제거)

결과를 확인하면 완벽하게 제거하지는 못했지만 심각한 flickering은 제거된 것을 확인할 수 있다.

---