频率域滤波【4】：带阻滤波

//快速傅里叶变换
void fft2Image(InputArray _src, OutputArray _dst)
{
	//得到Mat类型
	Mat src = _src.getMat();
	//判断位深
	CV_Assert(src.type() == CV_32FC1 || src.type() == CV_64FC1);
	CV_Assert(src.channels() == 1 || src.channels() == 2);
	int rows = src.rows;
	int cols = src.cols;
	//为了进行快速的傅里叶变换，我们经行和列的扩充,找到最合适扩充值
	Mat padded;
	int rPadded = getOptimalDFTSize(rows);
	int cPadded = getOptimalDFTSize(cols);
	//进行边缘扩充,扩充值为零
	copyMakeBorder(src, padded, 0, rPadded - rows, 0, cPadded - cols, BORDER_CONSTANT, Scalar::all(0));
	//快速的傅里叶变换（双通道：用于存储实部 和 虚部）
	dft(padded, _dst, DFT_COMPLEX_OUTPUT);
}
//求傅里叶变换的幅度谱 amplitude spectrum
void amplitudeSpectrum(InputArray _srcFFT, OutputArray _dstSpectrum)
{
	//判断傅里叶变换是两个通道
	CV_Assert(_srcFFT.channels() == 2);
	//分离通道
	vector<Mat> FFT2Channel;
	split(_srcFFT, FFT2Channel);
	//计算傅里叶变换的幅度谱 sqrt(pow(R,2)+pow(I,2))
	magnitude(FFT2Channel[0], FFT2Channel[1], _dstSpectrum);
}
//傅里叶谱的灰度级显示
Mat graySpectrum(Mat spectrum)
{
	Mat dst;
	log(spectrum + 1, dst);
	//归一化
	normalize(dst, dst, 0, 1, NORM_MINMAX);
	//为了进行灰度级显示，做类型转换
	dst.convertTo(dst, CV_8UC1, 255, 0);
	return dst;
}
//求相位谱 phase spectrum
Mat phaseSpectrum(Mat _srcFFT)
{
	//相位谱
	Mat phase;
	phase.create(_srcFFT.size(), CV_64FC1);
	//分离通道
	vector<Mat> FFT2Channel;
	split(_srcFFT, FFT2Channel);
	//计算相位谱
	for (int r = 0; r<phase.rows; r++)
	{
		for (int c = 0; c < phase.cols; c++)
		{
			//实部 虚部
			double real = FFT2Channel[0].at<double>(r, c);
			double imaginary = FFT2Channel[1].at<double>(r, c);
			// atan2 的返回值范围 [0,180] [-180,0]
			phase.at<double>(r, c) = atan2(imaginary, real);
		}
	}
	return phase;
}


//带阻
//枚举类型：理想带阻滤波器,，巴特沃斯带阻滤波器，高斯带阻滤波器
enum BRFILTER_TYPE { IBR_FILTER = 0, BBR_FILTER = 1, GBR_FILTER = 2 };
//构建带阻滤波器
Mat createBRFilter(Size _size, Point center, float _radius, float _bandWidth, int type, int n = 2)
{
	Mat _bpFilter(_size, CV_32FC1);
	int rows = _size.height;
	int cols = _size.width;
	CV_Assert(_bandWidth / 2 <= _radius);
	//构造理想带阻滤波器
	if (type == IBR_FILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
				float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
				if (sqrt(norm2) <= _radius + _bandWidth / 2 && sqrt(norm2) >= _radius - _bandWidth / 2)
					_bpFilter.at<float>(r, c) = 0;
				else
					_bpFilter.at<float>(r, c) = 1;
			}
		}
	}
	//构造巴特沃斯带阻滤波器
	if (type == BBR_FILTER)
	{
		for (int r = 0; r < rows; r++)
		{
			for (int c = 0; c<cols; c++)
			{
				float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
				if (norm2 != pow(_radius, 2))
					_bpFilter.at<float>(r, c) = 1.0 / (1.0 + pow(sqrt(norm2)*_bandWidth / (norm2 - pow(_radius, 2)), 2 * n));
				else
					_bpFilter.at<float>(r, c) = 0;
			}
		}
	}
	//构造高斯带阻滤波
	if (type == GBR_FILTER)
	{
		for (int r = 0; r< rows; r++)
		{
			for (int c = 0; c < cols; c++)
			{
				float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
				if (norm2 != 0)
					_bpFilter.at<float>(r, c) = 1.0 - exp(-pow((norm2 - pow(_radius, 2)) / (_bandWidth*sqrt(norm2)), 2));
				else
					_bpFilter.at<float>(r, c) = 1;
			}
		}
	}
	return _bpFilter;
}

//带阻滤波器
Mat brFilter(Mat &src, Mat &F, Point maxLoc, int type, int radius,int bandWid)
{
	Mat brFilter;//带阻滤波器
	Mat F_brFilter;//带阻傅里叶变换
	Mat FbrSpectrum;//带阻傅里叶变换的傅里叶谱灰度级
	Mat result;

	/* -- 第四步：构造带阻滤波器 -- */
	brFilter = createBRFilter(F.size(), maxLoc, radius, bandWid, type);


	/*-- 第六步：带阻滤波器和图像快速傅里叶变换点乘 --*/
	F_brFilter.create(F.size(), F.type());
	for (int r = 0; r < F_brFilter.rows; r++)
	{
		for (int c = 0; c < F_brFilter.cols; c++)
		{
			//分别取出当前位置的快速傅里叶变换和理想带阻滤波器的值
			Vec2f F_rc = F.at<Vec2f>(r, c);
			float lpFilter_rc = brFilter.at<float>(r, c);
			//带阻滤波器和图像的快速傅里叶变换对应位置相乘
			F_brFilter.at<Vec2f>(r, c) = F_rc*lpFilter_rc;
		}
	}

	//带阻傅里叶变换的傅里叶谱
	amplitudeSpectrum(F_brFilter, FbrSpectrum);
	//带阻傅里叶谱的灰度级的显示
	FbrSpectrum = graySpectrum(FbrSpectrum);
	/* -- 第七、八步：对带阻傅里叶变换执行傅里叶逆变换，并只取实部 -- */
	dft(F_brFilter, result, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);
	/* -- 第九步：同乘以(-1)^(x+y) -- */
	for (int r = 0; r < result.rows; r++)
	{
		for (int c = 0; c < result.cols; c++)
		{
			if ((r + c) % 2)
				result.at<float>(r, c) *= -1;
		}
	}
	//注意将结果转换 CV_8U 类型
	result.convertTo(result, CV_8UC1, 1.0, 0);
	/* -- 第十步：截取左上部分,大小等于输入图像的大小 --*/
	result = result(Rect(0, 0, src.cols, src.rows)).clone();
	return result;
}

void main()
{	/* -- 第一步：读入图像矩阵 -- */
	Mat image = imread("lena.jpg", 0);
	Mat fI;
	image.convertTo(fI, CV_32FC1);
	/* -- 第二步：每一个数乘以(-1)^(r+c) -- */
	for (int r = 0; r < fI.rows; r++)
	{
		for (int c = 0; c < fI.cols; c++)
		{
			if ((r + c) % 2)
				fI.at<float>(r, c) *= -1;
		}
	}
	/* -- 第三、四步：补零和快速傅里叶变换 -- */
	Mat F;//图像的快速傅里叶变换
	fft2Image(fI, F);
	//傅里叶谱
	Mat amplSpec;
	amplitudeSpectrum(F, amplSpec);
	//傅里叶谱的灰度级显示
	Mat spectrum = graySpectrum(amplSpec);
	//找到傅里叶谱的最大值的坐标
	Point maxLoc;
	minMaxLoc(spectrum, NULL, NULL, NULL, &maxLoc);

	Mat result;
	//低通
	//result  = LPFilter(image, F, maxLoc, 0, 25);

	//高通
	//result = ihpFilter(image, F, maxLoc, 0, 25);

	//带阻
	result = brFilter(image, F, maxLoc, 0, 50, 40);

	cout << "";
}

原图：

傅里叶：

带阻：

反变换：