維納濾波圖像處理

維納濾波處理模糊圖像

代碼連接: 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);
}