关键公式,《数字图像处理》第三版,P147:
1、频谱图具有周期性,归一化为图像大小;
2、频谱图具有对称性,以复数形式表示图像二维,幅值为实部与虚部平方和;
https://blog.csdn.net/daduzimama/article/details/80109139
https://www.zhihu.com/question/67234546/answer/403292397
https://www.cnblogs.com/Imageshop/archive/2019/10/08/11625895.html
示例一:参考HDevelop中的remove_texture_fft.hdev,采用OpenCV实现去除图像中固定纹理,实际上是构建了一个带阻滤波器;这里只是演示,只去除了部分纹理频率,没有完全去除干净,效果肯定没有HDevelop中的好:
Mat ConvertToFreq(Mat imgInput)
{
int w = getOptimalDFTSize(imgInput.cols);
int h = getOptimalDFTSize(imgInput.rows);
Mat padded;
copyMakeBorder(imgInput, padded, 0, h - imgInput.rows, 0, w - imgInput.cols, BORDER_CONSTANT, Scalar::all(0));
padded.convertTo(padded, CV_32FC1);
for (int i = 0; i < padded.rows; i++)//中心化操作,其余操作和上一篇博客的介绍一样
{
float* ptr = padded.ptr<float>(i);
for (int j = 0; j < padded.cols; j++) ptr[j] *= pow(-1, i + j);
}
Mat plane[] = { padded,Mat::zeros(padded.size(),CV_32F) };
Mat imgFreqComplex;
merge(plane, 2, imgFreqComplex);
dft(imgFreqComplex, imgFreqComplex);
return imgFreqComplex;
}
Mat ConvertFromFreq(Size size, Mat imgInput)
{
Mat imgTemp;
idft(imgInput, imgTemp);
vector<Mat> plane;
split(imgTemp, plane);
magnitude(plane[0], plane[1], plane[0]);
Mat imgResult;
normalize(plane[0], plane[0], 0, 255, NORM_MINMAX, CV_8UC1);
Rect rect(0, 0, size.width, size.height);
try
{
imgResult = plane[0](rect);
}
catch (...)
{
imgResult = plane[0].clone();
}
return imgResult;
}
Mat ConvertVisualization(Mat imgInput)
{
vector<Mat> plane;
split(imgInput, plane);
magnitude(plane[0], plane[1], plane[0]);
plane[0] += Scalar::all(1);
log(plane[0], plane[0]);
normalize(plane[0], plane[0], 1, 0, NORM_MINMAX);
return plane[0];
}
void testFrequency1()
{
Mat imgInput = imread("plan_01.png", 0);
Mat imgFreqComplex = ConvertToFreq(imgInput);//中心化后复数形式的频谱图
Mat imgFreqVisual = ConvertVisualization(imgFreqComplex);//取对数后的幅值频谱图
Mat imgFreq8U;//归一化8位1通道的幅值频谱图
normalize(imgFreqVisual, imgFreq8U, 0, 255, NORM_MINMAX, CV_8UC1);
Mat imgFreqFilter;
blur(imgFreq8U, imgFreqFilter, Size(5, 5));
Mat imgFreqThresh;
threshold(imgFreqFilter, imgFreqThresh, 140, 255, THRESH_BINARY); //经验值
Mat imgMorph;
morphologyEx(imgFreqThresh, imgMorph, MORPH_DILATE, getStructuringElement(MORPH_ELLIPSE, Size(5, 5)));
vector<contour> vecContFind;
findContours(imgMorph, vecContFind, RETR_LIST, CHAIN_APPROX_NONE);
vector<contour> vecContSelect;//筛选排除掉中心区域大量的低频分量
for (size_t j = 0; j < vecContFind.size(); j++)
{
contour contTemp = vecContFind[j];
float fContArea = contourArea(contTemp);
if (fContArea < 300)//经验值
vecContSelect.push_back(contTemp);
}
Mat imgFreqMask = 255 - Contour2Mask(imgMorph.size(), vecContSelect);
Mat imgFreqBitwise;//去除纹理相关频率的频谱图
bitwise_and(imgFreqComplex, imgFreqComplex, imgFreqBitwise, imgFreqMask);
Mat imgOutput = ConvertFromFreq(imgInput.size(), imgFreqBitwise);//从频域转换回时域,已经归一化为CV_8UC1
}
效果图如下,分别对应原图、去除纹理后、原图对应的复数形式频谱图、去除纹理后对应的复数形式频谱图、取对数后的幅值频谱图、二值化后的幅值频谱图、膨胀后的幅值频谱图、排除掉中心区域的低频分量后的幅值频谱图。
示例二:在频域中实现高斯滤波,实现了高通低通高斯滤波,对应图像的锐化和模糊;
Mat GeneGaussFilter(Size size, float fSigma, FilterType type)
{
Mat imgGauss(size, CV_32FC2);
if (HIGHT_PASS == type)
{
for (int i = 0; i < imgGauss.rows; i++)
{
float* q = imgGauss.ptr<float>(i);
for (int j = 0; j < imgGauss.cols; j++)
{
float d = pow(i - imgGauss.rows / 2, 2) + pow(j - imgGauss.cols / 2, 2);
q[2 * j] = 1 - expf(-d / fSigma);
q[2 * j + 1] = 1 - expf(-d / fSigma);
}
}
}
else
{
for (int i = 0; i < imgGauss.rows; i++)
{
float* p = imgGauss.ptr<float>(i);
for (int j = 0; j < imgGauss.cols; j++)
{
float d = pow(i - imgGauss.rows / 2, 2) + pow(j - imgGauss.cols / 2, 2);
p[2 * j] = expf(-d / fSigma);
p[2 * j + 1] = expf(-d / fSigma);
}
}
}
return imgGauss;
}
void testFrequency2()
{
Mat imgInput = imread("lena.jpg", 0);
Mat imgFreqComplex = ConvertToFreq(imgInput);
Mat imgFreqGauss = GeneGaussFilter(imgFreqComplex.size(), 2 * 10 * 10, HIGHT_PASS);
Mat imgFreqBlur;
multiply(imgFreqComplex, imgFreqGauss, imgFreqBlur);//矩阵元素对应相乘法,注意,和矩阵相乘区分
Mat imgOutput = ConvertFromFreq(imgInput.size(), imgFreqBlur);
Mat imgVisual0 = ConvertVisualization(imgFreqComplex);
Mat imgVisual1 = ConvertVisualization(imgFreqGauss);
Mat imgVisual2 = ConvertVisualization(imgFreqBlur);
}
效果图如下,分别对应原图、滤波后、原图的频谱、高斯核的频谱。