文档章节

维纳滤波图像处理

KyJason
 KyJason
发布于 2016/12/12 20:06
字数 1336
阅读 225
收藏 0

维纳滤波处理模糊图像

代码链接: https://github.com/hoijui/DIP/tree/master/ex04/src/main/native 感谢作者。

滤波结果

wiener.h

//wiener.h
#ifndef WIENER_H
#define WIENER_H
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

// function headers of not yet implemented functions

/**
 * Function applies inverse filter to restore a degraded image.
 * @param degraded  degraded input image
 * @param filter  filter which caused degradation
 * @return  restored output image
 */
Mat inverseFilter(Mat& degraded, Mat& filter);

/**
 * Function applies wiener filter to restore a degraded image.
 * @param degraded  degraded input image
 * @param filter  filter which caused degradation
 * @param snr  signal to noise ratio of the input image
 * @return  restored output image
 */
Mat wienerFilter(Mat& degraded, Mat& filter, double snr);

/**
 * Creates a filter kernel matrix of a certain type and size.
 * @param kernel  where to store the generated kernel to
 * @param kSize  size of the kernel (kSize * kSize)
 * @param name  name of the type of kernel to generate; possible values:
 *   "gaussian", "uniform"
 */
void createKernel(Mat& kernel, int kSize, string name = "gaussian");

/**
 * Performs a circular shift in (dx,dy) direction.
 * @param in  input matrix
 * @param out  circular shifted matrix
 * @param dx  shift in x-direction
 * @param dy  shift in y-direction
 */
void circShift(Mat& in, Mat& out, int dx, int dy);

// function headers of given functions

/**
 * Function degrades a given image with gaussian blur and additive gaussian noise.
 * @param img  input image
 * @param degradedImg  degraded output image
 * @param filterDev  standard deviation of kernel for gaussian blur
 * @param snr  signal to noise ratio for additive gaussian noise
 * @return the used gaussian kernel
 */
Mat degradeImage(Mat& img, Mat& degradedImg, double filterDev, double snr);

/**
 * Function displays image (after proper normalization).
 * @param win  Window name
 * @param img  Image that shall be displayed
 * @param cut  whether to cut or scale values outside of [0,255] range
 */
void showImage(const char* win, Mat img, bool cut = true);




#endif // WIENER_H

wiener.cpp

//wiener.cpp
#include <iostream>
#include <opencv2/opencv.hpp>

// function headers of not yet implemented functions
#include "wiener.h"

static void circShiftXXX(Mat& in, Mat& out, int dx, int dy){

    int dstI, dstJ;

    for(int i=0; i<in.rows;i++)
    {
        dstI = i+dx;
        if(dstI<0)
        {
            dstI += out.rows;
        } else if(dstI>=out.rows){
            dstI -= out.rows;
        }

        for(int j=0; j<in.cols;j++)
        {
            dstJ = j+dy;
            if(dstJ<0)
            {
                dstJ += out.cols;
            } else if(dstJ>=out.cols){
                dstJ -= out.cols;
            }

            out.at<float>(dstI,dstJ) = in.at<float>(i,j);
        }
    }
}

static Mat inverseAndWiener(Mat& s, Mat& p, double snr, bool inverse) {

    const bool wiener = !inverse;

    // Pad input image to avoid ringing artifacts along image borders.
    int bH = p.cols;
    int bV = p.rows;
    Mat sBorder;
    copyMakeBorder(s, sBorder, bV, bV, bH, bH, BORDER_REPLICATE);

    // Allocate some memory like it is going out of style.
    Mat pBigShifted = Mat::zeros(sBorder.size(), CV_32F);
    Mat P = Mat::zeros(sBorder.size(), CV_32F);
    Mat S = Mat::zeros(sBorder.size(), CV_32F);
    Mat OApprox = Mat::zeros(sBorder.size(), CV_32F);
    Mat oApprox = Mat::zeros(sBorder.size(), CV_32F);

    // Shift kernel.
    const int pHalf = p.rows / 2;
    circShiftXXX(p, pBigShifted, -pHalf, -pHalf);

    // Transform shifted kernel and degrated input image into frequency domain.
    // Note: DFT_COMPLEX_OUTPUT means that we want the complex result to be stored
    //       in a two-channel matrix as opposed to the default compressed output.
    dft(pBigShifted, P, DFT_COMPLEX_OUTPUT);
    dft(sBorder, S, DFT_COMPLEX_OUTPUT);

    if (inverse) {
        const double epsilon = 0.05f;

        // Remove frequencies whose magnitude is below epsilon * max(freqKernel magnitude).
        double maxMagnitude;
        minMaxLoc(abs(P), 0, &maxMagnitude);

        const double threshold =  maxMagnitude * epsilon;
        for (int ri = 0; ri < P.rows; ri++) {
            for (int ci = 0; ci < P.cols; ci++) {
                if (norm(P.at<Vec2f>(ri, ci)) < threshold) {
                    P.at<Vec2f>(ri, ci) = threshold;
                }
            }
        }
    }

    // OpenCV only provides a multiplication operation for complex matrices, so we need
    // to calculate the inverse (1/H) of our filter spectrum first. Since it is complex
    // we need to compute 1/H = H*/(HH*) = H*/(Re(H)^2+Im(H)^2), where H* -> complex conjugate of H.

    // Multiply spectrum of the degrated image with the complex conjugate of the frequency spectrum
    // of the filter.
    const bool conjFreqKernel = true;
    mulSpectrums(S, P, OApprox, DFT_COMPLEX_OUTPUT, conjFreqKernel); // I * H*

    // Split kernel spectrum into real and imaginary parts.
    Mat PChannels[] = {Mat::zeros(sBorder.size(), CV_32F), Mat::zeros(sBorder.size(), CV_32F)};
    split(P, PChannels); // 0:real, 1:imaginary

    // Calculate squared magnitude (Re(H)^2 + Im(H)^2) of filter spectrum.
    Mat freqKernelSqMagnitude = Mat::zeros(sBorder.rows, sBorder.cols, CV_32F);
    magnitude(PChannels[0], PChannels[1], freqKernelSqMagnitude); // freqKernelSqMagnitude = magnitude
    pow(PChannels[0], 2, freqKernelSqMagnitude); // freqKernelSqMagnitude = magnitude^2 = Re(H)^2 + Im(H)^2

    if (wiener) {
        // Add 1 / SNR^2 to the squared filter kernel magnitude.
        freqKernelSqMagnitude += 1 / pow(snr, 2.0);
    }

    // Split frequency spectrum of degradedPadded image into real and imaginary parts.
    Mat OApproxChannels[] = {Mat::zeros(sBorder.size(), CV_32FC1), Mat::zeros(sBorder.size(), CV_32F)};
    split(OApprox, OApproxChannels);

    // Divide each plane by the squared magnitude of the kernel frequency spectrum.
    // What we have done up to this point: (I * H*) / (Re(H)^2 + Im(H)^2) = I/H
    divide(OApproxChannels[0], freqKernelSqMagnitude, OApproxChannels[0]); // Re(I) / (Re(H)^2 + Im(H)^2)
    divide(OApproxChannels[1], freqKernelSqMagnitude, OApproxChannels[1]); // Im(I) / (Re(H)^2 + Im(H)^2)

    // Merge real and imaginary parts of the image frequency spectrum.
    merge(OApproxChannels, 2, OApprox);

    // Inverse DFT.
    // Note: DFT_REAL_OUTPUT means that we want the output to be a one-channel matrix again.
    dft(OApprox, oApprox, DFT_INVERSE | DFT_SCALE | DFT_REAL_OUTPUT);

    // Crop output image to original size.
    oApprox = oApprox(Rect(bH, bV, oApprox.cols - (bH * 2), oApprox.rows - (bV * 2)));

    return oApprox;
}

Mat inverseFilter(Mat& degraded, Mat& filter) {
    return inverseAndWiener(degraded, filter, -1.0, true);
}

Mat wienerFilter(Mat& degraded, Mat& filter, double snr) {
    return inverseAndWiener(degraded, filter, snr, false);
}

void circShift(Mat& in, Mat& out, int dx, int dy) {

    const int h = in.rows;
    const int w = in.cols;

//	out = Mat::zeros(h, w, in.type());

    for (int y = 0; y < h; ++y) {
        int yNew = y + dy;
        if (yNew < 0) {
            yNew = yNew + h;
        } else if (yNew >= h) {
            yNew = yNew - h;
        }

        for (int x = 0; x < w; ++x) {
            int xNew = x + dx;
            if (xNew < 0) {
                xNew = xNew + w;
            } else if (xNew >= w) {
                xNew = xNew - w;
            }

            out.at<float>(yNew, xNew) = in.at<float>(y, x);
        }
    }
}

mian.cpp

//mian.cpp

#include<iostream>
#include<vector>
#include <opencv2/opencv.hpp>


#include "wiener.h"

using namespace cv;
using namespace std;

Mat degradeImage(Mat& img, Mat& degradedImg, double filterDev, double snr);
void showImage(const char* win, Mat img, bool cut);
int main()
{

    const char* win_1 = "Original Image";
    const char* win_2 = "Degraded Image";
    const char* win_3 = "Restored Image: Inverse filter";
    const char* win_4 = "Restored Image: Wiener filter";
    namedWindow(win_1);
    namedWindow(win_2);
    namedWindow(win_3);
    namedWindow(win_4);

    // load image, path in argv[1]
    cout << "load image" << endl;
    Mat img = imread("lena.jpg", 0);
    // convert U8 to 32F
    img.convertTo(img, CV_32FC1);
    cout << " > done" << endl;

    // show and safe gray-scale version of original image
    showImage(win_1, img);
    imwrite("original.png", img);

    // degrade image
    cout << "degrade image" << endl;
    double filterDev = 9;
    double snr = 10; //10000;
    Mat degradedImg;
    Mat gaussKernel = degradeImage(img, degradedImg, filterDev, snr);
    cout << " > done" << endl;

    // show and safe degraded image
    showImage(win_2, degradedImg);
    imwrite("degraded.png", degradedImg);

    // inverse filter
    cout << "inverse filter" << endl;
    Mat restoredImgInverseFilter = inverseFilter(degradedImg, gaussKernel);
    cout << " > done" << endl;

    // show and safe restored image
    showImage(win_3, restoredImgInverseFilter);
    imwrite("restored_inverse.png", restoredImgInverseFilter);

    // wiener filter
    cout << "wiener filter" << endl;
    Mat restoredImgWienerFilter = wienerFilter(degradedImg, gaussKernel, snr);
    cout << " > done" << endl;

    // show and safe restored image
    showImage(win_4, restoredImgWienerFilter, false);
    imwrite("restored_wiener.png", restoredImgWienerFilter);

    // wait
    waitKey(0);

    return 0;
}

/*
      *************************
      ***   GIVEN FUNCTIONS ***
      *************************
      */

Mat degradeImage(Mat& img, Mat& degradedImg, double filterDev, double snr) {

    int kSize = round(filterDev * 3)*2 - 1;

    Mat gaussKernel = getGaussianKernel(kSize, filterDev, CV_32FC1);
    gaussKernel = gaussKernel * gaussKernel.t();
    filter2D(img, degradedImg, -1, gaussKernel);

    Mat mean, stddev;
    meanStdDev(img, mean, stddev);

    Mat noise = Mat::zeros(img.rows, img.cols, CV_32FC1);
    randn(noise, 0, stddev.at<double>(0) / snr);
    degradedImg = degradedImg + noise;
    threshold(degradedImg, degradedImg, 255, 255, CV_THRESH_TRUNC);
    threshold(degradedImg, degradedImg, 0, 0, CV_THRESH_TOZERO);

    return gaussKernel;
}

void showImage(const char* win, Mat img, bool cut) {

    Mat tmp = img.clone();

    if (tmp.channels() == 1) {
        if (cut) {
            threshold(tmp, tmp, 255, 255, CV_THRESH_TRUNC);
            threshold(tmp, tmp, 0, 0, CV_THRESH_TOZERO);
        } else {
            normalize(tmp, tmp, 0, 255, CV_MINMAX);
        }

        tmp.convertTo(tmp, CV_8UC1);
    } else {
        tmp.convertTo(tmp, CV_8UC3);
    }
    imshow(win, tmp);
}


© 著作权归作者所有

共有 人打赏支持
KyJason
粉丝 9
博文 62
码字总数 38933
作品 0
杭州
程序员
加载中

评论(3)

KyJason
KyJason

引用来自“黑驴啸西风”的评论

大神,能解释一下circShiftXXX函数的作用么,我模拟了一下,直接懵逼了。。。
circshift 就是把像素向dx dy的方向平移过去,多出来的部分补在缺失的地方
黑驴啸西风
大神,能解释一下circShiftXXX函数的作用么,我模拟了一下,直接懵逼了。。。
chl1990
chl1990
hjjj
【图像处理】参数维纳滤波(Parametric Wiener Filter)

实验要求   (a) 编写一个给图像中添加高斯噪声的程序,程序的输入参数为噪声的均值与方差。   (b) 编写程序实现公式(5.6-11)所示的污损滤波;   (c) 如图5.26(b)所示,对图像5.26(a) 进...

u013165921
01/15
0
0
如何还原模糊的车牌----维纳反卷积滤波算法

Wiener deconv opencv学习笔记---维纳反卷积滤波算法演示 ,算法参考[维基百科:https://en.wikipedia.org/wiki/Wiener_deconvolution 本例将使用算法实现模糊车牌的清晰化处理,使用Wiener反卷...

gaoshine
07/10
0
0
【工具使用系列】关于 MATLAB 图像复原 & 图像重建,你需要知道的事

如何进行图像复原 & 图像重建 添加噪声 用imnoise函数为图像添加噪声 用给定分布产生空间随机噪声 周期噪声 估计噪声参数 消除噪声 仅有噪声的复原 空间噪声滤波器 自适应空间滤波器 周期噪声...

AllenMoore
01/27
5
0
基于深层神经网络的语音 增强方法研究

近年来,随着深层神经网络(在语音识别领域的成功应用,给了语音增强任务的研宄人员很多启发。的深层非线性结构可以被设计成一个精细的降噪滤波器。同时基于大数据训练,可以充分学 习带噪语...

chenxiao60
2016/06/07
187
0
[Python图像处理] 四.图像平滑之均值滤波、方框滤波、高斯滤波及中值滤波

版权声明:本文为博主原创文章,转载请注明CSDN博客源地址!共同学习,一起进步~ https://blog.csdn.net/Eastmount/article/details/82216380 该系列文章是讲解Python OpenCV图像处理知识,前...

Eastmount
09/02
0
0

没有更多内容

加载失败,请刷新页面

加载更多

Vue学习资料

一直以为Vue是依赖nodejs的。 作为前端也可以耦合性就很低了。 //npm包管理器 进行管理npm install vue//初始化一个项目vue init//本地调试npm run dev//编译完成 ...

大灰狼wow
29分钟前
1
0
fullcalendar重新渲染

uiCalendarConfig.calendars.lesson_calendar.fullCalendar('removeEvents');var ym = uiCalendarConfig.calendars.lesson_calendar.fullCalendar('getView').title;$scope.get_lesson(y......

人来疯啊
33分钟前
1
0
多渠道打包总结

https://www.jianshu.com/p/2130db7584c8 https://blog.csdn.net/u011153817/article/details/50772496...

塔塔米
42分钟前
1
0
android -------- Data Binding的使用 ( 六) 自定义属性

今天来说说DataBinding在自定义属性的使用 默认的android命名空间下,我们会发现并不是所有的属性都能直接通过data binding进行设置,比如margin,padding,还有自定义View的各种属性。 默认...

切切歆语
48分钟前
1
0
收邮件 下载附件

uses IdMessage, IdMessageParts, IdAttachment, IdGlobalProtocols, ...;procedure SaveAttachmentsFromFile(FileName: String)var IdMessage: TIdMessage; MsgPart: Ti......

vga
54分钟前
1
0

没有更多内容

加载失败,请刷新页面

加载更多

返回顶部
顶部