维纳滤波处理模糊图像
代码链接: 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);
}
4132
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
