【opencv】教程代码 —ImgProc (8) 通过傅里叶变换对图像去除周期性噪声

8. periodic_noise_removing_filter.cpp通过傅里叶变换对图像去除周期性噪声

/**
* @brief 你将学会如何在傅立叶域中去掉周期性噪声
* @author Karpushin Vladislav, karpushin@ngs.ru, https://github.com/VladKarpushin
*/
#include <iostream>   // 引入C++的标准输入输出库
#include "opencv2/highgui.hpp"   // 引入opencv库，用于图像的高级GUI操作
#include <opencv2/imgcodecs.hpp>   // 引入opencv库，用于读写图片
#include <opencv2/imgproc.hpp>  // 引入opencv库，用于图像处理


using namespace cv;  // 使用OpenCV的命名空间
using namespace std;   // 使用标准的命名空间


//声明函数原型
void fftshift(const Mat& inputImg, Mat& outputImg);
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H);
void synthesizeFilterH(Mat& inputOutput_H, Point center, int radius);
void calcPSD(const Mat& inputImg, Mat& outputImg, int flag = 0);


//定义一个字符串变量，用于存储命令行参数
const String keys =
"{help h usage ? |             | print this message   }"
"{@image          |period_input.jpg | input image name     }"
;


int main(int argc, char* argv[])   // 定义主函数，输入参数为命令行的参数数量及详细内容
{
    CommandLineParser parser(argc, argv, keys);   // 创建一个命令行解析对象
    string strInFileName = parser.get<String>("@image");   // 获取输入的图像文件名
    samples::addSamplesDataSearchSubDirectory("doc/tutorials/imgproc/periodic_noise_removing_filter/images");   //添加搜索子目录


    Mat imgIn = imread(samples::findFile(strInFileName), IMREAD_GRAYSCALE);  // 以灰度模式读取输入的图像文件
    if (imgIn.empty()) // 检查图像是否加载成功
    {
        cout << "ERROR : Image cannot be loaded..!!" << endl;   // 如果加载不成功，输出错误信息
        return -1;  // 并返回错误代码
    }


    imshow("Image corrupted", imgIn);   // 显示原始图像
    imgIn.convertTo(imgIn, CV_32F);   // 将图像格式转换为浮点型


    // 只处理偶数像素的图像
    Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);   // 创建一个矩形，该矩形的大小为偶数像素的图像大小
    imgIn = imgIn(roi);   // 将矩形应用于图像，只保留在矩形内的像素


    // 计算图像的功率谱密度(PSD)
    Mat imgPSD;   // 定义一个Mat对象，用于存储PSD
    calcPSD(imgIn, imgPSD);   // 计算PSD
    fftshift(imgPSD, imgPSD);   // 进行快速傅立叶变换
    normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);   // 对PSD进行归一化处理


    // 计算傅立叶变换的H函数
    Mat H = Mat(roi.size(), CV_32F, Scalar(1));   // 定义一个Mat对象，用于存储H函数
    const int r = 21;   // 定义一个变量，表示半径
    synthesizeFilterH(H, Point(705, 458), r);   // 生成H函数
    synthesizeFilterH(H, Point(850, 391), r);
    synthesizeFilterH(H, Point(993, 325), r);


    // 进行滤波操作
    Mat imgOut;   // 定义一个Mat对象，用于存储滤波后的图像
    fftshift(H, H);    // 对H函数进行快速傅立叶变换
    filter2DFreq(imgIn, imgOut, H);   // 对图像进行滤波操作
    imgOut.convertTo(imgOut, CV_8U);   // 将滤波后的图像类型转换为8位无符号整型
    normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);   // 对滤波后的图像进行归一化处理
    imwrite("result.jpg", imgOut);   // 将滤波后的图像写入到result.jpg文件中
    imwrite("PSD.jpg", imgPSD);   // 写入PSD到PSD.jpg文件中
    fftshift(H, H);   // 对H函数进行快速傅立叶变换
    normalize(H, H, 0, 255, NORM_MINMAX);   //对H函数进行归一化处理
    imshow("Debluring", imgOut);   // 显示去噪后的图像
    imwrite("filter.jpg", H);   // 将H函数写入到filter.jpg文件中
    waitKey(0);   // 等待用户按键，如果用户不按键，程序会一直等待
    return 0;   // 主函数结束，返回0
}


//! [fftshift]
void fftshift(const Mat& inputImg, Mat& outputImg) // 定义一个fftshift函数，用于将图像的四个象限进行位置交换，参数为输入的图像inputImg和输出的图像outputImg
{
    outputImg = inputImg.clone();  // 克隆输入图像到输出图像
    int cx = outputImg.cols / 2;   // 计算图像列数的一半，即x轴上的中心点
    int cy = outputImg.rows / 2;   // 计算图像行数的一半，即y轴上的中心点


    //定义四个矩形块，均以左上角为起始点，大小为图像的四分之一，分别代表四个象限
    Mat q0(outputImg, Rect(0, 0, cx, cy));    // q0是左上角的区域
    Mat q1(outputImg, Rect(cx, 0, cx, cy));   // q1是右上角的区域
    Mat q2(outputImg, Rect(0, cy, cx, cy));   // q2是左下角的区域
    Mat q3(outputImg, Rect(cx, cy, cx, cy));  // q3是右下角的区域


    Mat tmp;   // 定义一个临时的Mat对象，用于在交换象限时临时存储数据
    q0.copyTo(tmp);  // 将q0复制到临时对象tmp
    q3.copyTo(q0);   // 将q3复制到q0
    tmp.copyTo(q3);  // 将tmp复制到q3, 到此q0与q3象限的图像数据已交换完毕


    q1.copyTo(tmp);  // 将q1复制到临时对象tmp
    q2.copyTo(q1);   // 将q2复制到q1
    tmp.copyTo(q2);  // 将tmp复制到q2,  此步操作在于交换q1与q2象限的图像数据，此时fftshift函数执行完毕
}
//! [fftshift]


//! [filter2DFreq]
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H) // 定义一个filter2DFreq函数，用于在傅立叶域进行滤波，参数包括输入图像inputImg，输出图像outputImg，和频率域滤波函数H
{
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };  // 创建两个Mat对象数组，即平面，第一个平面复制输入图像的像素数据，第二个平面为同等大小的零矩阵，用于存储傅立叶变换的虚部
    Mat complexI;  // 定义一个复数Mat，用于存放复数形式的图像像素数据
    merge(planes, 2, complexI);  // 将前面创建的两个平面合并为一个Mat，即复数形式的图像像素数据
    
    dft(complexI, complexI, DFT_SCALE);  // 对复数平面的图像像素进行傅立叶变换


    Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };  // 对滤波器也进行相同的平面创建操作
    Mat complexH;   // 定义一个复数Mat，用于存放复数形式的滤波器像素数据
    merge(planesH, 2, complexH);  // 将滤波器的两个平面合并为一个Mat，即复数形式的滤波器像素数据
    
    Mat complexIH;  // 定义一个结果复数Mat，用于存放傅立叶变换后的图像像素与滤波器的乘积
    mulSpectrums(complexI, complexH, complexIH, 0);  // 对两个进行傅立叶变换后的Mat对象进行乘积操作，结果赋值为complexIH


    idft(complexIH, complexIH);  //对乘积操作后的结果执行逆傅立叶变换，以便恢复原图像
    
    split(complexIH, planes);  // 将逆变换后的复数Mat拆分为两个平面
    outputImg = planes[0];  // 结果中实数部分的平面即为滤波后的图像
}
//! [filter2DFreq]


//! [synthesizeFilterH]
void synthesizeFilterH(Mat& inputOutput_H, Point center, int radius)  // 定义synthesizeFilterH函数，用于在滤波器中添加半径为radius的平面波，该平面波位于center, 参数包括滤波器的引用inputOutput_H，平面波的中心点center和半径radius。
{
    Point c2 = center, c3 = center, c4 = center;  // 定义中心点的映射点c2,c3,c4，这里以中心点center为基准
    c2.y = inputOutput_H.rows - center.y;  // c2是center相对于y轴的对称点
    c3.x = inputOutput_H.cols - center.x;  // c3是center相对于x轴的对称点
    c4 = Point(c3.x,c2.y);  // c4是center相对于原点的对称点


    //在滤波器上画出4个圆，并将其值设为0，其中圆心为center, 半径为radius，这四个圆对应频域上应该抑制的频率分量，对应抑制图像的噪声。
    circle(inputOutput_H, center, radius, 0, -1, 8);
    circle(inputOutput_H, c2, radius, 0, -1, 8);
    circle(inputOutput_H, c3, radius, 0, -1, 8);
    circle(inputOutput_H, c4, radius, 0, -1, 8);
}
//! [synthesizeFilterH]


// 函数通过FFT计算功率谱密度(PSD)
// flag = 0 表示返回PSD
// flag = 1 表示返回log(PSD) 
//! [calcPSD]
void calcPSD(const Mat& inputImg, Mat& outputImg, int flag)
{
   // 定义一个包含两个Mat的数组planes，第一个平面为输入图像的复制物，第二个平面为输入图像尺寸相同的全零矩阵
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    Mat complexI;  // 定义一个复数型的Mat对象complexI
    merge(planes, 2, complexI); // 调用merge函数将平面数组合并成一个多通道数组complexI


    dft(complexI, complexI); // 利用dft对合并后的数组进行傅立叶变换


    split(complexI, planes); // 利用split函数将傅立叶变换后的复数型数组拆分到两个通道，planes[0]为复数的实部，planes[1]为复数的虚部


    // 将两个平面的(0,0)像素值设为0，消除直流分量
    planes[0].at<float>(0) = 0;  
    planes[1].at<float>(0) = 0;


    // 计算PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
    Mat imgPSD;  // 定义一个Mat对象imgPSD，用来存放功率谱密度图像
    magnitude(planes[0], planes[1], imgPSD);   // 计算每个像素的幅度值，其值为复数的模，即sqrt(Re^2 + Im^2)
    pow(imgPSD, 2, imgPSD); // 对幅度矩阵全体元素再求平方，得到PSD
    outputImg = imgPSD; // 将计算得到的功率谱密度图赋值给输出图像outputImg


    // flag=1时，PSD经过log处理，即logPSD = log(1 + PSD)，log对数变换可以扩展低值区域的对比度，使图像更易观察。
    if (flag)
    {
        Mat imglogPSD;
        imglogPSD = imgPSD + Scalar::all(1); // calc log(1 + PSD) = logPSD
        log(imglogPSD, imglogPSD); // 取PSD的对数，得到log(PSD+1)
        outputImg = imglogPSD; // 将计算得到的logPSD赋值给输出图像outputImg
    }
}
//!

代码的主要功能是在傅立叶域中消除图像的周期性噪声。这是通过动态生成和应用滤波函数实现的，通过对影响图像质量的特定频率成分进行消除，达到了消除周期性噪声的效果。代码中同时实现了傅立叶逆变换，以便将处理结果重新转换回图像空间，生成最终的去噪结果。