opencv 低通滤波器

#include "stdafx.h"
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>


using namespace cv;
using namespace std;

void shiftDFT(cv::Mat& fImage);
void create_lowpass_filter(cv::Mat &dft_Filter, float fc, int type);
void lowpass_dft(const cv::Mat& image, cv::Mat& dft_filter_img, float fc);

int main()
{
	float scale = 1.0;
	cv::Mat image = cv::imread("FileName.tif", cv::IMREAD_GRAYSCALE);

	if (image.empty())
	{
		return 0;
	}

	cv::Mat src;
	cv::resize(image, src, cv::Size(), scale, scale);



	imwrite("src.tif", src);
	float fc = 0.3;
	cv::Mat dft_filter_img;
	lowpass_dft(src, dft_filter_img, fc);
	imwrite("dft_filter_img.bmp", dft_filter_img);

	return 0;
}

void create_lowpass_filter(cv::Mat &dft_Filter, float fc, int type)
{
	cv::Point centre = cv::Point(dft_Filter.rows / 2, dft_Filter.cols / 2);
	double radius;

	float x, y;
	int state = -1;
	x = dft_Filter.cols * 0.5;
	y = dft_Filter.rows * 0.5;
	double tempD;
	double D0 = fc * min(dft_Filter.rows, dft_Filter.cols) / 2.0;
	for (int i = 0; i < dft_Filter.rows; i++)
	{
		float* ptr = dft_Filter.ptr<float>(i);
		for (int j = 0; j < dft_Filter.cols; j++)
		{
			if (i > y && j > x)
			{
				state = 3;
			}
			else if (i > y)
			{
				state = 1;
			}
			else if (j > x)
			{
				state = 2;
			}
			else
			{
				state = 0;
			}

			switch (state)
			{


			case 0:
				tempD = pow(i, (double)2.0) + pow(j, (double)2.0); tempD = sqrt(tempD); break;
			case 1:
				tempD = pow((dft_Filter.rows - i), (double)2.0) + pow(j, (double)2.0); tempD = sqrt(tempD); break;
			case 2:
				tempD = pow(i, (double)2.0) + pow((dft_Filter.cols - j), (double)2.0); tempD = sqrt(tempD); break;
			case 3:
				tempD = pow((dft_Filter.rows - i), (double)2.0) + pow((dft_Filter.cols - j), (double)2.0); tempD = sqrt(tempD); break;
			default:
				break;
			}

			switch (type)
			{
			case 0://理想滤波器
				 if(tempD <= D0)  
				 {  
					 ptr[j * 2] = 1.0;
					 ptr[j * 2 + 1] = 1.0;
				 }  
				 else  
				 {  
					 ptr[j * 2] = 0.0;
					 ptr[j * 2 + 1] = 0.0;
				 }  
				break;
			case 1://2阶巴特沃思低通滤波器传递函数  
				tempD = 1 / (1 + pow(tempD / D0, 2 * 2));
				ptr[j * 2] = tempD;
				ptr[j * 2 + 1] = tempD;
				break;
			case 2://二维高斯低通滤波器传递函数  
				tempD = exp(-0.5 * pow(tempD / D0, 2));
				ptr[j * 2] = tempD;
				ptr[j * 2 + 1] = tempD;
				break;
			case 3://衰减系数为2的二维指数低通滤波器传递函数 
				tempD = exp(-pow(tempD / D0, 2));
				ptr[j * 2] = tempD;
				ptr[j * 2 + 1] = tempD;
				break;

			}
		}
	}


}


void lowpass_dft(const cv::Mat& image, cv::Mat& dft_filter_img, float fc)
{
	int M = getOptimalDFTSize(image.rows);
	int N = getOptimalDFTSize(image.cols);
	cv::Mat padded;
	copyMakeBorder(image, padded, 0, M - image.rows, 0, N - image.cols, BORDER_CONSTANT, Scalar::all(0));

	cv::Mat planes[] = { cv::Mat_<float>(padded), cv::Mat::zeros(padded.size(), CV_32F) };
	cv::Mat complexImg(planes[0].size(), planes[0].type(), CV_32FC2);
	merge(planes, 2, complexImg);
	///傅里叶正变换
	dft(complexImg, complexImg);

	//显示频域图
	split(complexImg, planes);
	magnitude(planes[0], planes[1], planes[0]);
	Mat mag = planes[0].clone();
	mag += Scalar::all(1);
	log(mag, mag);
	shiftDFT(mag);
	normalize(mag, mag, 0, 1, NORM_MINMAX);

	///创建滤波器
	cv::Mat filterImg = cv::Mat(complexImg.size(), complexImg.type());
	create_lowpass_filter(filterImg, fc, 2);

	显示滤波器图
	split(filterImg, planes);
	magnitude(planes[0], planes[1], planes[0]);
	Mat magFilter = planes[0].clone();
	magFilter += Scalar::all(1);
	log(magFilter, magFilter);
	shiftDFT(magFilter);
	normalize(magFilter, magFilter, 0, 1, NORM_MINMAX);
	
	//同时显示频域图和滤波器图
	cv::Mat dftDebug = mag + magFilter;

	//频域滤波
	mulSpectrums(complexImg, filterImg, complexImg, 0);
	///傅里叶逆变换
	idft(complexImg, complexImg, cv::DFT_SCALE);

	split(complexImg, planes);
	dft_filter_img.create(complexImg.size(), CV_8UC1);
	planes[0].convertTo(dft_filter_img, CV_8UC1);
	dft_filter_img = dft_filter_img(cv::Rect(0, 0, image.cols, image.rows));
}

void shiftDFT(cv::Mat& fImage)
{
	Mat tmp, q0, q1, q2, q3;
	//first crop the image, if it has an odd number of rows or columns  
	int cx = fImage.cols / 2;
	int cy = fImage.rows / 2;


	// rearrange the quadrants of Fourier image  
	// so that the origin is at the image center  


	q0 = fImage(Rect(0, 0, cx, cy));
	q1 = fImage(Rect(cx, 0, cx, cy));
	q2 = fImage(Rect(0, cy, cx, cy));
	q3 = fImage(Rect(cx, cy, cx, cy));


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


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