OpenCV  4.9.0-dev Open Source Computer Vision
Searching...
No Matches
Out-of-focus Deblur Filter

Prev Tutorial: Image Segmentation with Distance Transform and Watershed Algorithm
Next Tutorial: Motion Deblur Filter

Compatibility OpenCV >= 3.0

# Goal

In this tutorial you will learn:

• what a degradation image model is
• what the PSF of an out-of-focus image is
• how to restore a blurred image
• what is a Wiener filter

# Theory

Note
The explanation is based on the books [108] and [325]. Also, you can refer to Matlab's tutorial Image Deblurring in Matlab and the article SmartDeblur.
The out-of-focus image on this page is a real world image. The out-of-focus was achieved manually by camera optics.

## What is a degradation image model?

Here is a mathematical model of the image degradation in frequency domain representation:

$S = H\cdot U + N$

where $$S$$ is a spectrum of blurred (degraded) image, $$U$$ is a spectrum of original true (undegraded) image, $$H$$ is a frequency response of point spread function (PSF), $$N$$ is a spectrum of additive noise.

The circular PSF is a good approximation of out-of-focus distortion. Such a PSF is specified by only one parameter - radius $$R$$. Circular PSF is used in this work.

## How to restore a blurred image?

The objective of restoration (deblurring) is to obtain an estimate of the original image. The restoration formula in frequency domain is:

$U' = H_w\cdot S$

where $$U'$$ is the spectrum of estimation of original image $$U$$, and $$H_w$$ is the restoration filter, for example, the Wiener filter.

## What is the Wiener filter?

The Wiener filter is a way to restore a blurred image. Let's suppose that the PSF is a real and symmetric signal, a power spectrum of the original true image and noise are not known, then a simplified Wiener formula is:

$H_w = \frac{H}{|H|^2+\frac{1}{SNR}}$

where $$SNR$$ is signal-to-noise ratio.

So, in order to recover an out-of-focus image by Wiener filter, it needs to know the $$SNR$$ and $$R$$ of the circular PSF.

# Source code

You can find source code in the samples/cpp/tutorial_code/ImgProc/out_of_focus_deblur_filter/out_of_focus_deblur_filter.cpp of the OpenCV source code library.

#include <iostream>
using namespace cv;
using namespace std;
void help();
void calcPSF(Mat& outputImg, Size filterSize, int R);
void fftshift(const Mat& inputImg, Mat& outputImg);
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H);
void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr);
const String keys =
"{help h usage ? | | print this message }"
"{image |original.jpg | input image name }"
"{SNR |100 | signal to noise ratio}"
;
int main(int argc, char *argv[])
{
help();
CommandLineParser parser(argc, argv, keys);
if (parser.has("help"))
{
parser.printMessage();
return 0;
}
int R = parser.get<int>("R");
int snr = parser.get<int>("SNR");
string strInFileName = parser.get<String>("image");
if (!parser.check())
{
parser.printErrors();
return 0;
}
Mat imgIn;
if (imgIn.empty()) //check whether the image is loaded or not
{
cout << "ERROR : Image cannot be loaded..!!" << endl;
return -1;
}
Mat imgOut;
// it needs to process even image only
Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
//Hw calculation (start)
Mat Hw, h;
calcPSF(h, roi.size(), R);
calcWnrFilter(h, Hw, 1.0 / double(snr));
//Hw calculation (stop)
// filtering (start)
filter2DFreq(imgIn(roi), imgOut, Hw);
// filtering (stop)
imgOut.convertTo(imgOut, CV_8U);
normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
imshow("Original", imgIn);
imshow("Debluring", imgOut);
imwrite("result.jpg", imgOut);
waitKey(0);
return 0;
}
void help()
{
cout << "2018-07-12" << endl;
cout << "DeBlur_v8" << endl;
cout << "You will learn how to recover an out-of-focus image by Wiener filter" << endl;
}
void calcPSF(Mat& outputImg, Size filterSize, int R)
{
Mat h(filterSize, CV_32F, Scalar(0));
Point point(filterSize.width / 2, filterSize.height / 2);
circle(h, point, R, 255, -1, 8);
Scalar summa = sum(h);
outputImg = h / summa[0];
}
void fftshift(const Mat& inputImg, Mat& outputImg)
{
outputImg = inputImg.clone();
int cx = outputImg.cols / 2;
int cy = outputImg.rows / 2;
Mat q0(outputImg, Rect(0, 0, cx, cy));
Mat q1(outputImg, Rect(cx, 0, cx, cy));
Mat q2(outputImg, Rect(0, cy, cx, cy));
Mat q3(outputImg, Rect(cx, cy, cx, cy));
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}
void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
{
Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI, DFT_SCALE);
Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
Mat complexH;
merge(planesH, 2, complexH);
Mat complexIH;
mulSpectrums(complexI, complexH, complexIH, 0);
idft(complexIH, complexIH);
split(complexIH, planes);
outputImg = planes[0];
}
void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr)
{
Mat h_PSF_shifted;
fftshift(input_h_PSF, h_PSF_shifted);
Mat planes[2] = { Mat_<float>(h_PSF_shifted.clone()), Mat::zeros(h_PSF_shifted.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI);
split(complexI, planes);
Mat denom;
pow(abs(planes[0]), 2, denom);
denom += nsr;
divide(planes[0], denom, output_G);
}
Designed for command line parsing.
Definition utility.hpp:820
Template matrix class derived from Mat.
Definition mat.hpp:2230
n-dimensional dense array class
Definition mat.hpp:812
Creates a full copy of the array and the underlying data.
MatSize size
Definition mat.hpp:2160
void copyTo(OutputArray m) const
Copies the matrix to another one.
int cols
Definition mat.hpp:2138
bool empty() const
Returns true if the array has no elements.
int rows
the number of rows and columns or (-1, -1) when the matrix has more than 2 dimensions
Definition mat.hpp:2138
void convertTo(OutputArray m, int rtype, double alpha=1, double beta=0) const
Converts an array to another data type with optional scaling.
Template class for 2D rectangles.
Definition types.hpp:444
Size_< _Tp > size() const
size (width, height) of the rectangle
Template class for specifying the size of an image or rectangle.
Definition types.hpp:335
_Tp height
the height
Definition types.hpp:363
_Tp width
the width
Definition types.hpp:362
void split(const Mat &src, Mat *mvbegin)
Divides a multi-channel array into several single-channel arrays.
void mulSpectrums(InputArray a, InputArray b, OutputArray c, int flags, bool conjB=false)
Performs the per-element multiplication of two Fourier spectrums.
void divide(InputArray src1, InputArray src2, OutputArray dst, double scale=1, int dtype=-1)
Performs per-element division of two arrays or a scalar by an array.
Scalar sum(InputArray src)
Calculates the sum of array elements.
void merge(const Mat *mv, size_t count, OutputArray dst)
Creates one multi-channel array out of several single-channel ones.
void idft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Calculates the inverse Discrete Fourier Transform of a 1D or 2D array.
void dft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Performs a forward or inverse Discrete Fourier transform of a 1D or 2D floating-point array.
void pow(InputArray src, double power, OutputArray dst)
Raises every array element to a power.
std::string String
Definition cvstd.hpp:151
#define CV_8U
Definition interface.h:73
#define CV_32F
Definition interface.h:78
softfloat abs(softfloat a)
Absolute value.
Definition softfloat.hpp:444
@ circle
Definition gr_skig.hpp:62
int main(int argc, char *argv[])
Definition highgui_qt.cpp:3
"black box" representation of the file storage associated with a file on disk.
Definition core.hpp:102
STL namespace.

# Explanation

An out-of-focus image recovering algorithm consists of PSF generation, Wiener filter generation and filtering a blurred image in frequency domain:

// it needs to process even image only
Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
//Hw calculation (start)
Mat Hw, h;
calcPSF(h, roi.size(), R);
calcWnrFilter(h, Hw, 1.0 / double(snr));
//Hw calculation (stop)
// filtering (start)
filter2DFreq(imgIn(roi), imgOut, Hw);
// filtering (stop)

A function calcPSF() forms a circular PSF according to input parameter radius $$R$$:

void calcPSF(Mat& outputImg, Size filterSize, int R)
{
Mat h(filterSize, CV_32F, Scalar(0));
Point point(filterSize.width / 2, filterSize.height / 2);
circle(h, point, R, 255, -1, 8);
Scalar summa = sum(h);
outputImg = h / summa[0];
}

A function calcWnrFilter() synthesizes the simplified Wiener filter $$H_w$$ according to the formula described above:

void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr)
{
Mat h_PSF_shifted;
fftshift(input_h_PSF, h_PSF_shifted);
Mat planes[2] = { Mat_<float>(h_PSF_shifted.clone()), Mat::zeros(h_PSF_shifted.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI);
split(complexI, planes);
Mat denom;
pow(abs(planes[0]), 2, denom);
denom += nsr;
divide(planes[0], denom, output_G);
}

A function fftshift() rearranges the PSF. This code was just copied from the tutorial Discrete Fourier Transform:

void fftshift(const Mat& inputImg, Mat& outputImg)
{
outputImg = inputImg.clone();
int cx = outputImg.cols / 2;
int cy = outputImg.rows / 2;
Mat q0(outputImg, Rect(0, 0, cx, cy));
Mat q1(outputImg, Rect(cx, 0, cx, cy));
Mat q2(outputImg, Rect(0, cy, cx, cy));
Mat q3(outputImg, Rect(cx, cy, cx, cy));
Mat tmp;
q0.copyTo(tmp);
q3.copyTo(q0);
tmp.copyTo(q3);
q1.copyTo(tmp);
q2.copyTo(q1);
tmp.copyTo(q2);
}

A function filter2DFreq() filters the blurred image in the frequency domain:

void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
{
Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
Mat complexI;
merge(planes, 2, complexI);
dft(complexI, complexI, DFT_SCALE);
Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
Mat complexH;
merge(planesH, 2, complexH);
Mat complexIH;
mulSpectrums(complexI, complexH, complexIH, 0);
idft(complexIH, complexIH);
split(complexIH, planes);
outputImg = planes[0];
}

# Result

Below you can see the real out-of-focus image:

And the following result has been computed with $$R$$ = 53 and $$SNR$$ = 5200 parameters:

The Wiener filter was used, and values of $$R$$ and $$SNR$$ were selected manually to give the best possible visual result. We can see that the result is not perfect, but it gives us a hint to the image's content. With some difficulty, the text is readable.

Note
The parameter $$R$$ is the most important. So you should adjust $$R$$ first, then $$SNR$$.
Sometimes you can observe the ringing effect in a restored image. This effect can be reduced with several methods. For example, you can taper input image edges.

You can also find a quick video demonstration of this on YouTube.