1. 程式人生 > >基於opencv的理想低通濾波器和巴特沃斯低通濾波器

基於opencv的理想低通濾波器和巴特沃斯低通濾波器

首先看個圖瞭解下什麼是理想低通濾波器公式和圖是轉自Rolin的專欄

低通濾波器

        1.理想的低通濾波器


       其中,D0表示通帶的半徑。D(u,v)的計算方式也就是兩點間的距離,很簡單就能得到。

       使用低通濾波器所得到的結果如下所示。低通濾波器濾除了高頻成分,所以使得影象模糊。由於理想低通濾波器的過度特性過於急峻,所以會產生了振鈴現象。

        2.巴特沃斯低通濾波器

       同樣的,D0表示通帶的半徑,n表示的是巴特沃斯濾波器的次數。隨著次數的增加,振鈴現象會越來越明顯。

   

void ideal_Low_Pass_Filter(Mat src){
	Mat img;
	cvtColor(src, img, CV_BGR2GRAY);
	imshow("img",img);
	//調整影象加速傅立葉變換
	int M = getOptimalDFTSize(img.rows);
	int N = getOptimalDFTSize(img.cols);
	Mat padded;
	copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));
	//記錄傅立葉變換的實部和虛部
	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexImg;
	merge(planes, 2, complexImg);
	//進行傅立葉變換
	dft(complexImg, complexImg);
	//獲取影象
	Mat mag = complexImg;
	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));//這裡為什麼&上-2具體檢視opencv文件
	//其實是為了把行和列變成偶數 -2的二進位制是11111111.......10 最後一位是0
	//獲取中心點座標
	int cx = mag.cols / 2;
	int cy = mag.rows / 2;
	//調整頻域
	Mat tmp;
	Mat q0(mag, Rect(0, 0, cx, cy));
	Mat q1(mag, Rect(cx, 0, cx, cy));
	Mat q2(mag, Rect(0, cy, cx, cy));
	Mat q3(mag, Rect(cx, cy, cx, cy));

	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);

	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
	//Do為自己設定的閥值具體看公式
	double D0 = 60;
	//處理按公式保留中心部分
	for (int y = 0; y < mag.rows; y++){
		double* data = mag.ptr<double>(y);
		for (int x = 0; x < mag.cols; x++){
			double d = sqrt(pow((y - cy),2) + pow((x - cx),2));
			if (d <= D0){
				
			}
			else{
				data[x] = 0;
			}
		}
	}
	//再調整頻域
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
	//逆變換
	Mat invDFT, invDFTcvt;
	idft(mag, invDFT, DFT_SCALE | DFT_REAL_OUTPUT); // Applying IDFT
	invDFT.convertTo(invDFTcvt, CV_8U);
	imshow("理想低通濾波器", invDFTcvt);
}

void Butterworth_Low_Paass_Filter(Mat src){
	int n = 1;//表示巴特沃斯濾波器的次數
	//H = 1 / (1+(D/D0)^2n)
	Mat img;
	cvtColor(src, img, CV_BGR2GRAY);
	imshow("img", img);
	//調整影象加速傅立葉變換
	int M = getOptimalDFTSize(img.rows);
	int N = getOptimalDFTSize(img.cols);
	Mat padded;
	copyMakeBorder(img, padded, 0, M - img.rows, 0, N - img.cols, BORDER_CONSTANT, Scalar::all(0));

	Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
	Mat complexImg;
	merge(planes, 2, complexImg);

	dft(complexImg, complexImg);

	Mat mag = complexImg;
	mag = mag(Rect(0, 0, mag.cols & -2, mag.rows & -2));

	int cx = mag.cols / 2;
	int cy = mag.rows / 2;

	Mat tmp;
	Mat q0(mag, Rect(0, 0, cx, cy));
	Mat q1(mag, Rect(cx, 0, cx, cy));
	Mat q2(mag, Rect(0, cy, cx, cy));
	Mat q3(mag, Rect(cx, cy, cx, cy));

	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);

	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);

	double D0 = 100;

	for (int y = 0; y < mag.rows; y++){
		double* data = mag.ptr<double>(y);
		for (int x = 0; x < mag.cols; x++){
			//cout << data[x] << endl;
			double d = sqrt(pow((y - cy), 2) + pow((x - cx), 2));
			//cout << d << endl;
			double h = 1.0 / (1 + pow(d / D0, 2 * n));
			if (h <= 0.5){
				data[x] = 0;
			}
			else{
				//data[x] = data[x]*0.5;
				//cout << h << endl;
			}
			
			//cout << data[x] << endl;
		}
	}
	q0.copyTo(tmp);
	q3.copyTo(q0);
	tmp.copyTo(q3);
	q1.copyTo(tmp);
	q2.copyTo(q1);
	tmp.copyTo(q2);
	//逆變換
	Mat invDFT, invDFTcvt;
	idft(complexImg, invDFT, DFT_SCALE | DFT_REAL_OUTPUT); // Applying IDFT
	invDFT.convertTo(invDFTcvt, CV_8U);
	imshow("巴特沃斯低通濾波器", invDFTcvt);
}