# 频率域滤波【1】：低通滤波

//快速傅里叶变换
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;
//为了进行快速的傅里叶变换，我们经行和列的扩充,找到最合适扩充值
//进行边缘扩充,扩充值为零
//快速的傅里叶变换（双通道：用于存储实部 和 虚部）
}
//求傅里叶变换的幅度谱 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 LPFILTER_TYPE { ILP_FILTER = 0, BLP_FILTER = 1, GLP_FILTER = 2 };

//构造低通滤波器
Mat createLPFilter(Size size, Point center, float radius, int type, int n = 2)
{
Mat lpFilter = Mat::zeros(size, CV_32FC1);
int rows = size.height;
int cols = size.width;
return lpFilter;
//构造理想低通滤波器
if (type == ILP_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);
lpFilter.at<float>(r, c) = 1;
else
lpFilter.at<float>(r, c) = 0;
}
}
}
//构造巴特沃斯低通滤波器
if (type == BLP_FILTER)
{
for (int r = 0; r < rows; r++)
{
for (int c = 0; c<cols; c++)
{
lpFilter.at<float>(r, c) = float(1.0 / (1.0 + pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / radius, 2.0*n)));
}
}
}
//构造高斯低通滤波
if (type == GLP_FILTER)
{
for (int r = 0; r< rows; r++)
{
for (int c = 0; c < cols; c++)
{
lpFilter.at<float>(r, c) = float(exp(-(pow(c - center.x, 2.0) + pow(r - center.y, 2.0)) / (2 * pow(radius, 2.0))));
}
}
}
return lpFilter;
}

Mat LPFilter(Mat &src,Mat &F, Point maxLoc, int type,int radius)
{
Mat lpFilter;//低通滤波器
Mat F_lpFilter;//低通傅里叶变换
Mat FlpSpectrum;//低通傅里叶变换的傅里叶谱灰度级
Mat result;

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

//低通傅里叶变换的傅里叶谱
amplitudeSpectrum(F_lpFilter, FlpSpectrum);
//低通傅里叶谱的灰度级的显示
FlpSpectrum = graySpectrum(FlpSpectrum);
/* -- 第七、八步：对低通傅里叶变换执行傅里叶逆变换，并只取实部 -- */
dft(F_lpFilter, 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 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  = LPFilter(image, F, maxLoc, 0, 25);

cout << "";
}