OpenCV-一维频域滤波器

作者:翟天保Steven
版权声明:著作权归作者所有,商业转载请联系作者获得授权,非商业转载请注明出处

之前分享了许多图像频域滤波器的功能,如理想滤波器、高斯滤波器、巴特沃斯滤波器等等,都是二维形式的写法;有粉丝提出,不知道一维数据该怎么进行类似的处理。其实搞懂原理,二维降一维也是非常方便的,为了方便大家使用,本文将用C++和OpenCV实现一维频域滤波器功能。

示例为巴特沃斯,如果想换成别的,只要把滤波核更换即可。

二维与一维的不同之处:

// 图像边界处理(一维)
cv::Mat image_make_border_onedim(cv::Mat &src)
{
    int h = cv::getOptimalDFTSize(src.rows); // 获取DFT变换的最佳高度

    cv::Mat padded;
    // 常量法扩充图像边界,常量 = 0
    cv::copyMakeBorder(src, padded, 0, h - src.rows, 0, 0, cv::BORDER_CONSTANT, cv::Scalar::all(0));
    padded.convertTo(padded, CV_32FC1);

    return padded;
}

// fft变换后进行频谱搬移(一维)
void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1)
{
    // 以下的操作是移动图像  (零频移到中心)
    int cy = plane0.rows / 2;
    cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy));  // 元素坐标表示为(cx, cy)
    cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy));

    cv::Mat temp;
    part1_r.copyTo(temp);  //上与下交换位置(实部)
    part2_r.copyTo(part1_r);
    temp.copyTo(part2_r);

    cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy));  //元素坐标(cx,cy)
    cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy));

    part1_i.copyTo(temp);  //上与下交换位置(虚部)
    part2_i.copyTo(part1_i);
    temp.copyTo(part2_i);
}

// pow操作
Mat powZ(cv::InputArray src, double power) {
    cv::Mat dst;
    cv::pow(src, power, dst);
    return dst;
}

// sqrt操作
Mat sqrtZ(cv::InputArray src) {
    cv::Mat dst;
    cv::sqrt(src, dst);
    return dst;
}

// 频率域滤波
cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur)
{
    cv::Mat mask = scr == scr;
    scr.setTo(0.0f, ~mask);

    //创建通道,存储dft后的实部与虚部(CV_32F,必须为单通道数)
    cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) };

    cv::Mat complexIm;
    cv::merge(plane, 2, complexIm); // 合并通道 (把两个矩阵合并为一个2通道的Mat类容器)
    cv::dft(complexIm, complexIm); // 进行傅立叶变换,结果保存在自身

    // 分离通道(数组分离)
    cv::split(complexIm, plane);

    // 以下的操作是频域迁移
    fftshift_onedim(plane[0], plane[1]);

    // *****************滤波器函数与DFT结果的乘积****************
    cv::Mat blur_r, blur_i, BLUR;
    cv::multiply(plane[0], blur, blur_r);  // 滤波(实部与滤波器模板对应元素相乘)
    cv::multiply(plane[1], blur, blur_i);  // 滤波(虚部与滤波器模板对应元素相乘)
    cv::Mat plane1[] = { blur_r, blur_i };

    // 再次搬移回来进行逆变换
    fftshift_onedim(plane1[0], plane1[1]);
    cv::merge(plane1, 2, BLUR); // 实部与虚部合并

    cv::idft(BLUR, BLUR);       // idft结果也为复数
    BLUR = BLUR / BLUR.rows / BLUR.cols;

    cv::split(BLUR, plane);//分离通道,主要获取通道

    return plane[0];
}

// 巴特沃斯低通滤波核函数(一维)
cv::Mat butterworth_low_kernel_onedim(cv::Mat &scr, float sigma, int n)
{
    cv::Mat butterworth_low_pass(scr.size(), CV_32FC1); //,CV_32FC1
    float D0 = sigma;//半径D0越小,模糊越大;半径D0越大,模糊越小
    for (int i = 0; i < scr.rows; i++) {
        float d = abs(float(i - scr.rows / 2));
        butterworth_low_pass.at(i, 0) = 1.0f / (1.0f + pow(d / D0, 2 * n));
    }
    return butterworth_low_pass;
}

// 巴特沃斯低通滤波(一维)
cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n)
{
    // H = 1 / (1+(D/D0)^2n)   n表示巴特沃斯滤波器的次数
    // 阶数n=1 无振铃和负值    阶数n=2 轻微振铃和负值  阶数n=5 明显振铃和负值   阶数n=20 与ILPF相似
    cv::Mat padded = image_make_border_onedim(src);
    cv::Mat butterworth_kernel = butterworth_low_kernel_onedim(padded, d0, n);
    cv::Mat result = frequency_filter(padded, butterworth_kernel);
    return result;
}
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include<opencv2 opencv.hpp>
#include<ctime>

using namespace cv;
using namespace std;

// &#x8BFB;&#x53D6;&#x6570;&#x636E;
cv::Mat ReadOneDimData(string &file)
{
    cv::Mat result;
    std::ifstream in(file);
    std::string str;

    getline(in, str);
    while (str != "")
    {
        int a = str.find(',');
        string front = str.substr(0, a);
        string back = str.substr(a + 1);
        float wave = stof(front);
        float h = stof(back);
        cv::Mat temp = (cv::Mat_<float>(1, 2) << wave, h);
        result.push_back(temp);
        getline(in, str);
    }

    return result;
}

// &#x663E;&#x793A;&#x6570;&#x636E;&#x66F2;&#x7EBF;
cv::Mat ShowDataCurve(cv::Mat input)
{
    // &#x5F52;&#x4E00;&#x5316;
    cv::Mat curve = input.clone();
    cv::normalize(curve, curve, 0, 255, cv::NORM_MINMAX);
    curve.convertTo(curve, CV_8UC1);

    // &#x7ED8;&#x5236;&#x7B80;&#x6613;&#x66F2;&#x7EBF;&#x56FE;
    int size = input.rows;
    cv::Mat pic = cv::Mat::zeros(256, size, CV_8UC1);
    for (int i = 0; i < size; ++i)
    {
        uchar h = curve.at<uchar>(i, 0);
        pic.at<uchar>(255 - h, i) = 255;
    }
    return pic;
}

// &#x56FE;&#x50CF;&#x8FB9;&#x754C;&#x5904;&#x7406;&#xFF08;&#x4E00;&#x7EF4;&#xFF09;
cv::Mat image_make_border_onedim(cv::Mat &src)
{
    int h = cv::getOptimalDFTSize(src.rows); // &#x83B7;&#x53D6;DFT&#x53D8;&#x6362;&#x7684;&#x6700;&#x4F73;&#x9AD8;&#x5EA6;

    cv::Mat padded;
    // &#x5E38;&#x91CF;&#x6CD5;&#x6269;&#x5145;&#x56FE;&#x50CF;&#x8FB9;&#x754C;&#xFF0C;&#x5E38;&#x91CF; = 0
    cv::copyMakeBorder(src, padded, 0, h - src.rows, 0, 0, cv::BORDER_CONSTANT, cv::Scalar::all(0));
    padded.convertTo(padded, CV_32FC1);

    return padded;
}

// fft&#x53D8;&#x6362;&#x540E;&#x8FDB;&#x884C;&#x9891;&#x8C31;&#x642C;&#x79FB;&#xFF08;&#x4E00;&#x7EF4;&#xFF09;
void fftshift_onedim(cv::Mat &plane0, cv::Mat &plane1)
{
    // &#x4EE5;&#x4E0B;&#x7684;&#x64CD;&#x4F5C;&#x662F;&#x79FB;&#x52A8;&#x56FE;&#x50CF;  (&#x96F6;&#x9891;&#x79FB;&#x5230;&#x4E2D;&#x5FC3;)
    int cy = plane0.rows / 2;
    cv::Mat part1_r(plane0, cv::Rect(0, 0, 1, cy));  // &#x5143;&#x7D20;&#x5750;&#x6807;&#x8868;&#x793A;&#x4E3A;(cx, cy)
    cv::Mat part2_r(plane0, cv::Rect(0, cy, 1, cy));

    cv::Mat temp;
    part1_r.copyTo(temp);  //&#x4E0A;&#x4E0E;&#x4E0B;&#x4EA4;&#x6362;&#x4F4D;&#x7F6E;(&#x5B9E;&#x90E8;)
    part2_r.copyTo(part1_r);
    temp.copyTo(part2_r);

    cv::Mat part1_i(plane1, cv::Rect(0, 0, 1, cy));  //&#x5143;&#x7D20;&#x5750;&#x6807;(cx,cy)
    cv::Mat part2_i(plane1, cv::Rect(0, cy, 1, cy));

    part1_i.copyTo(temp);  //&#x4E0A;&#x4E0E;&#x4E0B;&#x4EA4;&#x6362;&#x4F4D;&#x7F6E;(&#x865A;&#x90E8;)
    part2_i.copyTo(part1_i);
    temp.copyTo(part2_i);
}

// pow&#x64CD;&#x4F5C;
Mat powZ(cv::InputArray src, double power) {
    cv::Mat dst;
    cv::pow(src, power, dst);
    return dst;
}

// sqrt&#x64CD;&#x4F5C;
Mat sqrtZ(cv::InputArray src) {
    cv::Mat dst;
    cv::sqrt(src, dst);
    return dst;
}

// &#x9891;&#x7387;&#x57DF;&#x6EE4;&#x6CE2;
cv::Mat frequency_filter(cv::Mat &scr, cv::Mat &blur)
{
    cv::Mat mask = scr == scr;
    scr.setTo(0.0f, ~mask);

    //&#x521B;&#x5EFA;&#x901A;&#x9053;&#xFF0C;&#x5B58;&#x50A8;dft&#x540E;&#x7684;&#x5B9E;&#x90E8;&#x4E0E;&#x865A;&#x90E8;&#xFF08;CV_32F&#xFF0C;&#x5FC5;&#x987B;&#x4E3A;&#x5355;&#x901A;&#x9053;&#x6570;&#xFF09;
    cv::Mat plane[] = { scr.clone(), cv::Mat::zeros(scr.size() , CV_32FC1) };

    cv::Mat complexIm;
    cv::merge(plane, 2, complexIm); // &#x5408;&#x5E76;&#x901A;&#x9053; &#xFF08;&#x628A;&#x4E24;&#x4E2A;&#x77E9;&#x9635;&#x5408;&#x5E76;&#x4E3A;&#x4E00;&#x4E2A;2&#x901A;&#x9053;&#x7684;Mat&#x7C7B;&#x5BB9;&#x5668;&#xFF09;
    cv::dft(complexIm, complexIm); // &#x8FDB;&#x884C;&#x5085;&#x7ACB;&#x53F6;&#x53D8;&#x6362;&#xFF0C;&#x7ED3;&#x679C;&#x4FDD;&#x5B58;&#x5728;&#x81EA;&#x8EAB;

    // &#x5206;&#x79BB;&#x901A;&#x9053;&#xFF08;&#x6570;&#x7EC4;&#x5206;&#x79BB;&#xFF09;
    cv::split(complexIm, plane);

    // &#x4EE5;&#x4E0B;&#x7684;&#x64CD;&#x4F5C;&#x662F;&#x9891;&#x57DF;&#x8FC1;&#x79FB;
    fftshift_onedim(plane[0], plane[1]);

    // *****************&#x6EE4;&#x6CE2;&#x5668;&#x51FD;&#x6570;&#x4E0E;DFT&#x7ED3;&#x679C;&#x7684;&#x4E58;&#x79EF;****************
    cv::Mat blur_r, blur_i, BLUR;
    cv::multiply(plane[0], blur, blur_r);  // &#x6EE4;&#x6CE2;&#xFF08;&#x5B9E;&#x90E8;&#x4E0E;&#x6EE4;&#x6CE2;&#x5668;&#x6A21;&#x677F;&#x5BF9;&#x5E94;&#x5143;&#x7D20;&#x76F8;&#x4E58;&#xFF09;
    cv::multiply(plane[1], blur, blur_i);  // &#x6EE4;&#x6CE2;&#xFF08;&#x865A;&#x90E8;&#x4E0E;&#x6EE4;&#x6CE2;&#x5668;&#x6A21;&#x677F;&#x5BF9;&#x5E94;&#x5143;&#x7D20;&#x76F8;&#x4E58;&#xFF09;
    cv::Mat plane1[] = { blur_r, blur_i };

    // &#x518D;&#x6B21;&#x642C;&#x79FB;&#x56DE;&#x6765;&#x8FDB;&#x884C;&#x9006;&#x53D8;&#x6362;
    fftshift_onedim(plane1[0], plane1[1]);
    cv::merge(plane1, 2, BLUR); // &#x5B9E;&#x90E8;&#x4E0E;&#x865A;&#x90E8;&#x5408;&#x5E76;

    cv::idft(BLUR, BLUR);       // idft&#x7ED3;&#x679C;&#x4E5F;&#x4E3A;&#x590D;&#x6570;
    BLUR = BLUR / BLUR.rows / BLUR.cols;

    cv::split(BLUR, plane);//&#x5206;&#x79BB;&#x901A;&#x9053;&#xFF0C;&#x4E3B;&#x8981;&#x83B7;&#x53D6;&#x901A;&#x9053;

    return plane[0];
}

// &#x5DF4;&#x7279;&#x6C83;&#x65AF;&#x4F4E;&#x901A;&#x6EE4;&#x6CE2;&#x6838;&#x51FD;&#x6570;&#xFF08;&#x4E00;&#x7EF4;&#xFF09;
cv::Mat butterworth_low_kernel_onedim(cv::Mat &scr, float sigma, int n)
{
    cv::Mat butterworth_low_pass(scr.size(), CV_32FC1); //&#xFF0C;CV_32FC1
    float D0 = sigma;//&#x534A;&#x5F84;D0&#x8D8A;&#x5C0F;&#xFF0C;&#x6A21;&#x7CCA;&#x8D8A;&#x5927;&#xFF1B;&#x534A;&#x5F84;D0&#x8D8A;&#x5927;&#xFF0C;&#x6A21;&#x7CCA;&#x8D8A;&#x5C0F;
    for (int i = 0; i < scr.rows; i++) {
        float d = abs(float(i - scr.rows / 2));
        butterworth_low_pass.at<float>(i, 0) = 1.0f / (1.0f + pow(d / D0, 2 * n));
    }
    return butterworth_low_pass;
}

// &#x5DF4;&#x7279;&#x6C83;&#x65AF;&#x4F4E;&#x901A;&#x6EE4;&#x6CE2;&#xFF08;&#x4E00;&#x7EF4;&#xFF09;
cv::Mat butterworth_low_pass_filter_onedim(cv::Mat &src, float d0, int n)
{
    // H = 1 / (1+(D/D0)^2n)   n&#x8868;&#x793A;&#x5DF4;&#x7279;&#x6C83;&#x65AF;&#x6EE4;&#x6CE2;&#x5668;&#x7684;&#x6B21;&#x6570;
    // &#x9636;&#x6570;n=1 &#x65E0;&#x632F;&#x94C3;&#x548C;&#x8D1F;&#x503C;    &#x9636;&#x6570;n=2 &#x8F7B;&#x5FAE;&#x632F;&#x94C3;&#x548C;&#x8D1F;&#x503C;  &#x9636;&#x6570;n=5 &#x660E;&#x663E;&#x632F;&#x94C3;&#x548C;&#x8D1F;&#x503C;   &#x9636;&#x6570;n=20 &#x4E0E;ILPF&#x76F8;&#x4F3C;
    cv::Mat padded = image_make_border_onedim(src);
    cv::Mat butterworth_kernel = butterworth_low_kernel_onedim(padded, d0, n);
    cv::Mat result = frequency_filter(padded, butterworth_kernel);
    return result;
}

int main(void)
{
    // &#x8BFB;&#x53D6;&#x6570;&#x636E;
    string datafile = "data4.txt";
    cv::Mat data = ReadOneDimData(datafile);

    // &#x83B7;&#x53D6;&#x539F;&#x59CB;&#x66F2;&#x7EBF;
    cv::Mat test = data.col(1).clone();
    cv::Mat show1 = ShowDataCurve(test);

    // &#x83B7;&#x53D6;&#x4F4E;&#x901A;&#x6570;&#x636E;&#x66F2;&#x7EBF;
    float D0 = 50.0f;
    Mat lowpass = butterworth_low_pass_filter_onedim(test, D0, 2);
    cv::Mat show2 = ShowDataCurve(lowpass);

    imshow("curve1", show1);
    imshow("curve2", show2);
    waitKey(0);

    system("pause");
    return 0;
}

</float></uchar></uchar></float></ctime></opencv2></vector></sstream></fstream></iostream>

图1 原始数据

图2 滤波处理后数据

为了方便大家看出差异,我写了一个简单的数据曲线显示函数,将就着看哈~上述示例数据为光谱仪的一组数据,原始数据是高频变化的信息,经过低通滤波器处理后,保留低频成分,过滤掉高频成分,因此曲线变得平滑了。

另外,如果我的代码有什么问题,欢迎大家提出异议批评指正,一同进步~

如果文章帮助到你了,可以点个赞让我知道,我会很快乐~加油!

Original: https://blog.csdn.net/zhaitianbao/article/details/125934465
Author: 翟天保Steven
Title: OpenCV-一维频域滤波器

原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/703359/

转载文章受原作者版权保护。转载请注明原作者出处!

(0)

大家都在看

亲爱的 Coder【最近整理,可免费获取】👉 最新必读书单  | 👏 面试题下载  | 🌎 免费的AI知识星球