文档章节

维纳滤波图像处理

KyJason
 KyJason
发布于 2016/12/12 20:06
字数 1336
阅读 205
收藏 0
点赞 1
评论 3

维纳滤波处理模糊图像

代码链接: 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
粉丝 8
博文 50
码字总数 36257
作品 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

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

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

AllenMoore ⋅ 01/27 ⋅ 0

基于深层神经网络的语音 增强方法研究

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

chenxiao60 ⋅ 2016/06/07 ⋅ 0

OpenCV学习之旅-简介

1、什么是图像,对图像进行处理是神马操作 一副图像可以定义为二维的函数z = f(x,y),其中x、y是其空间坐标,而其值z的大小就是函数在该点的灰度值。 例图.png 比如我用Matlab打开了一张256...

C6C ⋅ 2017/08/28 ⋅ 0

android 学习之图像处理系统(一)

android 学习之图像处理系统(一) 源代码:android图像处理系统1.0 下图是软件运行的截图: 本文只做了图像的打开和简单处理算法(图像的变亮、中值滤波、平滑滤波)的功能。其中,当软件运行...

长平狐 ⋅ 2012/10/08 ⋅ 0

RCP版的数字图像处理软件

这个项目是大三上学期学习数字图像处理时,我用Eclipse RCP开发的一个软件,实现了各种图像处理,包括反色,灰度化,直方图均衡化,空间域的低通滤波和高通滤波等等,频域率由于时间和能力限...

宅男潇涧 ⋅ 2013/05/19 ⋅ 2

开源图形库--CxImage

CxImage是一个可以用于MFC 的C++图像处理类库类,它可以打开,保存,显示,转换各种常见格式的图像文件,比如BMP, JPEG, GIF, PNG, TIFF, MNG, ICO, PCX, TGA, WMF, WBMP, JBG, J2K 等格式的...

匿名 ⋅ 2009/08/21 ⋅ 0

【图像处理】空间滤波、中值滤波(Spatial Filtering and Median Filtering)

实验要求   编写一个能够完成两幅图像之间加、减、乘、除四种算术运算的函数。另外,对于两幅图像的乘法,所编写的乘法程序还要能够完成一幅图像乘以一个常数的功能。使用图Fig1.10(4)和F...

u013165921 ⋅ 01/15 ⋅ 0

【图像处理】基于OpenCV底层实现的直方图匹配

image processing 系列: 【图像处理】图片旋转 【图像处理】高斯滤波、中值滤波、均值滤波 直方图匹配算法。又称直方图规定化。简单说。就是依据某函数、或者另外一张图片的引导,使得原图改...

技术mix呢 ⋅ 2017/11/09 ⋅ 0

PIL (Python图像处理库) 1.1.7 发布

PIL 1.1.7 发布了,该版本改进了对 PNG 压缩的处理,支持隔行扫描的PNG文件,改进了对各种 TGA 的支持,修复了一些错误等等。 图像处理工具包PIL(Python Image Library),该软件包提供了基本的...

红薯 ⋅ 2011/10/08 ⋅ 1

没有更多内容

加载失败,请刷新页面

加载更多

下一页

vim编辑模式、命令模式

编辑模式 vim要从一般模式进入编辑模式只要按字母 i 、I、a、A、o、O键就可以了 要从编辑模式回到一般模式按键盘上的Esc键即可。 按键 作用 i 在当前字符前插入 I 在光标所在行的行首插入 o ...

黄昏残影 ⋅ 26分钟前 ⋅ 0

OSChina 周五乱弹 —— 如果有一天不当程序员了

Osc乱弹歌单(2018)请戳(这里) 【今日歌曲】 @guanglun :分享off的单曲《我唱情歌给你听》 《我唱情歌给你听》- off 手机党少年们想听歌,请使劲儿戳(这里) @小小编辑 :#如果不做程序...

小小编辑 ⋅ 32分钟前 ⋅ 4

从 Confluence 5.3 及其早期版本中恢复空间

如果你需要从 Confluence 5.3 及其早期版本中的导出文件恢复到晚于 Confluence 5.3 的 Confluence 中的话。你可以使用临时的 Confluence 空间安装,然后将这个 Confluence 安装实例升级到你现...

honeymose ⋅ 今天 ⋅ 0

Java8新增的DateTimeFormatter与SimpleDateFormat的区别

两者最大的区别是,Java8的DateTimeFormatter也是线程安全的,而SimpleDateFormat并不是线程安全。 在并发环境下使用SimpleDateFormat 为了能够在多线程环境下使用SimpleDateFormat,有这三种...

人觉非常君 ⋅ 今天 ⋅ 0

多线程如何控制执行顺序

线程的生命周期说明: 当线程被创建并启动以后,它既不是一启动就进入了执行状态,也不是一直处于执行状态,在线程的生命周期中,它要经过新建(New)、就绪(Runnable)、运行(Running)、...

MarinJ_Shao ⋅ 今天 ⋅ 0

用ZBLOG2.3博客写读书笔记网站能创造今日头条的辉煌吗?

最近两年,著名的自媒体网站今日头条可以说是火得一塌糊涂,虽然从目前来看也遇到了一点瓶颈,毕竟发展到了一定的规模,继续增长就更加难了,但如今的今日头条规模和流量已经非常大了。 我们...

原创小博客 ⋅ 今天 ⋅ 0

MyBatis四大核心概念

本文讲解 MyBatis 四大核心概念(SqlSessionFactoryBuilder、SqlSessionFactory、SqlSession、Mapper)。 MyBatis 作为互联网数据库映射工具界的“上古神器”,训有四大“神兽”,谓之:Sql...

waylau ⋅ 今天 ⋅ 0

以太坊java开发包web3j简介

web3j(org.web3j)是Java版本的以太坊JSON RPC接口协议封装实现,如果需要将你的Java应用或安卓应用接入以太坊,或者希望用java开发一个钱包应用,那么用web3j就对了。 web3j的功能相当完整...

汇智网教程 ⋅ 今天 ⋅ 0

2个线程交替打印100以内的数字

重点提示: 线程的本质上只是一个壳子,真正的逻辑其实在“竞态条件”中。 举个例子,比如本题中的打印,那么在竞态条件中,我只需要一个方法即可; 假如我的需求是2个线程,一个+1,一个-1,...

Germmy ⋅ 今天 ⋅ 0

Django第一期

安装Django 去https://www.djangoproject.com/download/ 下载最新版的Django,然后解压放到Anaconda\Lib\site-packages目录下,然后cmd进入此目录,输入安装命令: python setup.py install ...

大不了敲一辈子代码 ⋅ 今天 ⋅ 0

没有更多内容

加载失败,请刷新页面

加载更多

下一页

返回顶部
顶部