【opencv 450 Image Processing】Periodic Noise Removing Filter周期性去噪滤波器

最新推荐文章于 2024-03-26 14:21:43 发布

十年一梦实验室

最新推荐文章于 2024-03-26 14:21:43 发布

阅读量2.9k

点赞数

分类专栏： opencv 文章标签： opencv c++ 计算机视觉

本文链接：https://blog.csdn.net/cxyhjl/article/details/123895739

版权

opencv 专栏收录该内容

87 篇文章 12 订阅

订阅专栏

Periodic Noise Removing Filter

周期性去噪滤波器 OpenCV: Periodic Noise Removing Filter

Goal

在本教程中，您将学习：如何去除傅里叶域中的周期性噪声

Theory

笔记

解释是基于书[95]。此页面上的图像是真实世界的图像。

周期性噪声在傅立叶域中产生尖峰，通常可以通过视觉分析检测到。

How to remove periodic noise in the Fourier domain?

如何去除傅里叶域中的周期性噪声？

通过频域滤波（frequency domain filtering）可以显着降低周期性噪声（Periodic noise）。在此页面上，我们使用具有适当半径的陷波抑制滤波器（notch reject filter）来完全封闭傅里叶域中的噪声尖峰。 陷波滤波器（notch filter）拒绝中心频率附近的（around a center frequency）预定义邻域中的频率。陷波滤波器的数量是任意的。凹口区域（notch areas ）的形状也可以是任意的（例如矩形或圆形）。在此页面上，我们使用三个圆形陷波滤光片。图像的功率谱致密化（Power spectrum densify ）用于噪声尖峰的视觉检测。

Source code

您可以在 OpenCV 源代码库的 samples/cpp/tutorial_code/ImgProc/periodic_noise_removing_filter/periodic_noise_removing_filter.cpp 中找到源代码。

/**
* @brief You will learn how to remove periodic noise in the Fourier domain
* 您将学习如何 去除 傅里叶域中的周期性噪声
* @author Karpushin Vladislav, karpushin@ngs.ru, https://github.com/VladKarpushin
*/
#include <iostream>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>

using namespace cv;
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);

int main()
{
    Mat imgIn = imread("input.jpg", IMREAD_GRAYSCALE);
    if (imgIn.empty()) //检查图像是否加载
    {
        cout << "ERROR : Image cannot be loaded..!!" << endl;
        return -1;
    }

    imgIn.convertTo(imgIn, CV_32F);

//! [main]
    // 它只需要处理偶数图像it needs to process even image only
    Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
    imgIn = imgIn(roi);

    // PSD calculation (start)
    Mat imgPSD;
    calcPSD(imgIn, imgPSD);//计算图像的功率谱密度
    fftshift(imgPSD, imgPSD);//fft变换后进行频谱搬移 
    normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);
    // PSD calculation (stop)

    //H calculation (start)
    Mat H = Mat(roi.size(), CV_32F, Scalar(1));
    const int r = 21;
    synthesizeFilterH(H, Point(705, 458), r);//根据中心频率和半径形成理想圆形陷波抑制滤波器的传递函数
    synthesizeFilterH(H, Point(850, 391), r);
    synthesizeFilterH(H, Point(993, 325), r);
    //H calculation (stop)
    // filtering (start)
    Mat imgOut;
    fftshift(H, H);//fft变换后进行频谱搬移 
    filter2DFreq(imgIn, imgOut, H);//imgIn在频域进行H滤波后输出时域imgOut
    // filtering (stop)
//! [main]

    imgOut.convertTo(imgOut, CV_8U);//转换格式用于显示
    normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);//归一化
    imwrite("result.jpg", imgOut);//显示结果
    imwrite("PSD.jpg", imgPSD);//
    fftshift(H, H);
    normalize(H, H, 0, 255, NORM_MINMAX);
    imwrite("filter.jpg", H);
    return 0;
}

//! [fftshift] fft变换后进行频谱搬移 
void fftshift(const Mat& inputImg, Mat& outputImg)
{
    outputImg = inputImg.clone();
    int cx = outputImg.cols / 2;
    int cy = outputImg.rows / 2;
    Mat q0(outputImg, Rect(0, 0, cx, cy));
    Mat q1(outputImg, Rect(cx, 0, cx, cy));
    Mat q2(outputImg, Rect(0, cy, cx, cy));
    Mat q3(outputImg, Rect(cx, cy, cx, cy));
    Mat tmp;
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);
    q2.copyTo(q1);
    tmp.copyTo(q2);
}
//! [fftshift]

//! [filter2DFreq] inputImg在频域进行H滤波后输出时域outputImg
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
{
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);
    dft(complexI, complexI, DFT_SCALE);//傅里叶变换

    Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
    Mat complexH;
    merge(planesH, 2, complexH);
    Mat complexIH;//计算 PSD（功率谱密度）
    mulSpectrums(complexI, complexH, complexIH, 0);//在频域  滤波

    idft(complexIH, complexIH);//逆向傅里叶变换
    split(complexIH, planes);
    outputImg = planes[0];
}
//! [filter2DFreq]

//! [synthesizeFilterH] 合成过滤器     inputOutput_H功率谱密度
void synthesizeFilterH(Mat& inputOutput_H, Point center, int radius)
{
    Point c2 = center, c3 = center, c4 = center;
    c2.y = inputOutput_H.rows - center.y;
    c3.x = inputOutput_H.cols - center.x;
    c4 = Point(c3.x,c2.y);
    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]

// Function calculates PSD(Power spectrum density) by fft with two flags
// 函数通过带有两个标志的 fft     计算 PSD（功率谱密度）
// flag = 0 means to return PSD
// flag = 1 means to return log(PSD)
//! [calcPSD]
void calcPSD(const Mat& inputImg, Mat& outputImg, int flag)
{
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);
    dft(complexI, complexI);
    split(complexI, planes);            // 实部planes[0] = Re(DFT(I)), 虚部planes[1] = Im(DFT(I))

    planes[0].at<float>(0) = 0;
    planes[1].at<float>(0) = 0;

    // compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
    Mat imgPSD;
    magnitude(planes[0], planes[1], imgPSD);		//imgPSD = sqrt(Power spectrum density)
    pow(imgPSD, 2, imgPSD);							//it needs ^2 in order to get PSD
    outputImg = imgPSD;

    // logPSD = log(1 + PSD)
    if (flag)
    {
        Mat imglogPSD;
        imglogPSD = imgPSD + Scalar::all(1);
        log(imglogPSD, imglogPSD);
        outputImg = imglogPSD;
    }
}
//! [calcPSD]

Explanation

频域滤波周期性降噪包括功率谱密度计算（用于噪声尖峰视觉检测）、陷波抑制滤波器合成和频率滤波：

  // it needs to process even image only
    Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
    imgIn = imgIn(roi);
    // PSD calculation (start)
    Mat imgPSD;
    calcPSD(imgIn, imgPSD);
    fftshift(imgPSD, imgPSD);
    normalize(imgPSD, imgPSD, 0, 255, NORM_MINMAX);
    // PSD calculation (stop)
    //H calculation (start)
    Mat H = Mat(roi.size(), CV_32F, Scalar(1));
    const int r = 21;
    synthesizeFilterH(H, Point(705, 458), r);
    synthesizeFilterH(H, Point(850, 391), r);
    synthesizeFilterH(H, Point(993, 325), r);
    //H calculation (stop)
    // filtering (start)
    Mat imgOut;
    fftshift(H, H);
    filter2DFreq(imgIn, imgOut, H);
    // filtering (stop)

函数 calcPSD() 计算图像的功率谱密度：

void calcPSD(const Mat& inputImg, Mat& outputImg, int flag)
{
    Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);
    dft(complexI, complexI);
    split(complexI, planes);            // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))
    planes[0].at<float>(0) = 0;
    planes[1].at<float>(0) = 0;
    // compute the PSD = sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)^2
    Mat imgPSD;
    magnitude(planes[0], planes[1], imgPSD);        //imgPSD = sqrt(Power spectrum density)
    pow(imgPSD, 2, imgPSD);                         //it needs ^2 in order to get PSD
    outputImg = imgPSD;
    // logPSD = log(1 + PSD)
    if (flag)
    {
        Mat imglogPSD;
        imglogPSD = imgPSD + Scalar::all(1);
        log(imglogPSD, imglogPSD);
        outputImg = imglogPSD;
    }
}

函数 synthesizeFilterH() 根据中心频率和半径形成理想圆形陷波抑制滤波器的传递函数：

void synthesizeFilterH(Mat& inputOutput_H, Point center, int radius)
{
    Point c2 = center, c3 = center, c4 = center;
    c2.y = inputOutput_H.rows - center.y;
    c3.x = inputOutput_H.cols - center.x;
    c4 = Point(c3.x,c2.y);
    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);
}

函数 filter2DFreq() 在频域中过滤图像。函数 fftshift() 和 filter2DFreq() 复制自教程 Out-of-focus Deblur Filter。

Result

下图显示了被各种频率的周期性噪声严重破坏的图像。

Image corrupted by periodic noise被周期性噪声破坏的图像

噪声分量很容易被视为下图所示功率谱密度（Power spectrum density）中的亮点（spikes尖峰）。

Power spectrum density showing periodic noise

显示周期性噪声的功率谱密度

The figure below shows a notch reject filter with an appropriate radius to completely enclose the noise spikes.

下图显示了一个具有适当半径的陷波抑制滤波器以完全包围噪声尖峰。

Notch reject filter 陷波抑制过滤器

使用陷波抑制滤波器处理图像的结果如下所示。

改善非常明显。与原始图像相比，该图像包含明显更少的可见周期性噪声。您还可以在 YouTube https://youtu.be/Qne51TcWwAc 上找到有关此过滤想法的快速视频演示。

十年一梦实验室

关注

0
点赞
踩
5

收藏

觉得还不错? 一键收藏
2
评论
【opencv 450 Image Processing】Periodic Noise Removing Filter周期性去噪滤波器

Periodic Noise Removing Filter周期性去噪滤波器 OpenCV: Periodic Noise Removing FilterGoal在本教程中，您将学习：如何去除傅里叶域中的周期性噪声Theory笔记解释是基于书[95]。此页面上的图像是真实世界的图像。周期性噪声在傅立叶域中产生尖峰，通常可以通过视觉分析检测到。How to remove periodic noise in the Fourier domain?如何去除傅里叶域中的...
复制链接

扫一扫