//快速傅里叶变换
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;
}
//自定义滤波器
bool drawing_box = false;
Point downPoint;
Rect rectFilter;
bool gotRectFilter = false;
void mouseRectHandler(int event, int x, int y, int, void*)
{
switch (event)
{
//按下鼠标左键
case CV_EVENT_LBUTTONDOWN:
drawing_box = true;
//记录起点
downPoint = Point(x, y);
break;
//移动鼠标
case CV_EVENT_MOUSEMOVE:
if (drawing_box)
{
//鼠标移动到 downPoint 的右下角
if (x >= downPoint.x && y >= downPoint.y)
{
rectFilter.x = downPoint.x;
rectFilter.y = downPoint.y;
rectFilter.width = x - downPoint.x;
rectFilter.height = y - downPoint.y;
}
//鼠标移到到 downPoint 的右上角
if (x >= downPoint.x && y <= downPoint.y)
{
rectFilter.x = downPoint.x;
rectFilter.y = y;
rectFilter.width = x - downPoint.x;
rectFilter.height = downPoint.y - y;
}
//鼠标移动到 downPoint 的左上角
if (x <= downPoint.x && y <= downPoint.y)
{
rectFilter.x = x;
rectFilter.y = y;
rectFilter.width = downPoint.x - x;
rectFilter.height = downPoint.y - y;
}
//鼠标移动到 downPoint 的左下角
if (x <= downPoint.x && y >= downPoint.y)
{
rectFilter.x = x;
rectFilter.y = downPoint.y;
rectFilter.width = downPoint.x - x;
rectFilter.height = y - downPoint.y;
}
}
break;
//松开鼠标左键
case CV_EVENT_LBUTTONUP:
drawing_box = false;
gotRectFilter = true;
break;
default:
break;
}
}
//自定义滤波器
Mat self_Filter(Mat &src, Mat &F, Mat &spectrum, Point maxLoc)
{
Mat selfFilter;自定义滤波器
Mat F_selfFilter;//自定义傅里叶变换
Mat FselfSpectrum;//自定义傅里叶变换的傅里叶谱灰度级
Mat result;
string windowName = "幅度谱的灰度级";
F_selfFilter = F.clone();
namedWindow(windowName, CV_WINDOW_AUTOSIZE);
setMouseCallback(windowName, mouseRectHandler, NULL);
for (;;)
{
spectrum(rectFilter).setTo(0);
/* -- 第六步:自定义滤波器与傅里叶变换点乘 -- */
F_selfFilter(rectFilter).setTo(Scalar::all(0));
imshow(windowName, spectrum);
//按下 "Esc" 键退出编辑
if (waitKey(10) == 27)
break;
}
//自定义傅里叶变换的傅里叶谱
amplitudeSpectrum(F_selfFilter, FselfSpectrum);
//自定义傅里叶谱的灰度级的显示
FselfSpectrum = graySpectrum(FselfSpectrum);
/* -- 第七、八步:对自定义傅里叶变换执行傅里叶逆变换,并只取实部 -- */
dft(F_selfFilter, 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);
//带通
//result = bpFilter(image, F, maxLoc, 0, 50, 40);
//自定义滤波器
result = self_Filter(image, F, spectrum, maxLoc);
cout << "";
}
原图:
傅里叶:
自定义:
反变换: