Wiener decovlotuin in cpp giving back the same image
I was trying to deblur a noisy image using wiener deconvolution. I found this code which added noise to an image and removed it as well. Modifying this code only i tried to implement the exact formula given on wiki. But the output is same as input
In the code i debugged a bit and found when i performed
magI=magI/x
values in magI all became 1.
Can anyone please check if the calculation i have done are correct or not?
If so how do i prevent values in magI becoming 1
PS: I have included the full code in case anyone wants to learn and implement the code. You can jump straight to wiener2 function as the error is in there.
#include <iostream>
#include "opencv2/core/core.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
using namespace cv;
using namespace std;
Mat wiener2(Mat I, Mat image_spectrum, int noise_stddev);
Mat padd_image(Mat I);
Mat get_spectrum(Mat I);
Mat get_dft(Mat I);
Mat with_noise(Mat image, int stddev);
Mat rand_noise(Mat I, int stddev);
Mat createavg(Size imsize) ;
void shift(Mat magI);
int main(int argc, char *argv[]) {
int noise_stddev=20;
string input_filename="blur.png", output_filename="write.png"; // Have a blurred image here
cout << "noise standard deviation: " << noise_stddev << "\n";
cout << "input file: " << input_filename << "\n";
Mat I = imread(input_filename, CV_LOAD_IMAGE_GRAYSCALE);
if(I.data==NULL){
cout << "Can't open file: " << input_filename << "\n";
return 2;
}
Mat raw_sample = imread("blur.png", CV_LOAD_IMAGE_GRAYSCALE);
if(raw_sample.data==NULL){
cout << "Can't open file: sample.bmp\n";
return 3;
}
Mat padded = padd_image(I);
Mat noisy;
noisy = with_noise(padded, noise_stddev);
Mat sample(padded.rows, padded.cols, CV_8U);
resize(raw_sample, sample, sample.size());
Mat spectrum = get_spectrum(sample); //to get signal spectrum of known image
Mat enhanced = wiener2(noisy, spectrum, noise_stddev);
imshow("image 1", noisy);
imshow("image 2", enhanced);
waitKey();
}
Mat createavg(Size imsize) {
Mat kernel = Mat(5,5,CV_32FC1,Scalar(0.04));
int w = imsize.width-kernel.cols;
int h = imsize.height-kernel.rows;
int r = w/2;
int l = imsize.width-kernel.cols -r;
int b = h/2;
int t = imsize.height-kernel.rows -b;
Mat ret;
copyMakeBorder(kernel,ret,t,b,l,r,BORDER_CONSTANT,Scalar::all(0));
return ret;
}
//inputs are the blurry image with noise , the original image power spectra , and standard deviation of the noise introduced
Mat wiener2(Mat final_noise, Mat image_spectrum, int noise_stddev){
Mat padded = padd_image(final_noise);
Mat noise = rand_noise(padded, noise_stddev);
Mat noise_spectrum = get_spectrum(noise);
Scalar padded_mean = mean(padded);
Mat planes[2];
Mat complexI = get_dft(padded);
split(complexI, planes); // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
Mat factor = (noise_spectrum / image_spectrum); //calculates the signal to noise ratio
//-----------------compute the frequency domain multiplier
Mat mask = createavg(padded.size()); //creating the kernel which initally prduced the blurred image
shift(mask);// shifting the filter
Mat mplane[] = {Mat_<float>(mask), Mat::zeros(mask.size(), CV_32F)};
Mat kernelcomplex;
merge(mplane, 2, kernelcomplex);
dft(kernelcomplex, kernelcomplex); // computing dft of kernel
split(kernelcomplex, mplane);// splitting the dft of kernel ...