| OpenCV 5.0.0alpha Open Source Computer Vision

Loading...
Searching...
No Matches
Out-of-focus Deblur Filter

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

Original author | Karpushin Vladislav
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 [109] and [328]. 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.

Circular point spread function

## 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>
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/imgcodecs.hpp"
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 }"
"{R |5 | radius }"
"{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");
samples::addSamplesDataSearchSubDirectory("doc/tutorials/imgproc/out_of_focus_deblur_filter/images");
if (!parser.check())
{
parser.printErrors();
return 0;
}
Mat imgIn;
imgIn = imread(samples::findFile( strInFileName ), IMREAD_GRAYSCALE);
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);
}
cv::CommandLineParser
Designed for command line parsing.
Definition utility.hpp:890
cv::Mat_
Template matrix class derived from Mat.
Definition mat.hpp:2518
cv::Mat
n-dimensional dense array class
Definition mat.hpp:951
cv::Mat::clone
CV_NODISCARD_STD Mat clone() const
Creates a full copy of the array and the underlying data.
cv::Mat::size
MatSize size
Definition mat.hpp:2448
cv::Mat::copyTo
void copyTo(OutputArray m) const
Copies the matrix to another one.
cv::Mat::cols
int cols
Definition mat.hpp:2425
cv::Mat::empty
bool empty() const
Returns true if the array has no elements.
cv::Mat::rows
int rows
the number of rows and columns or (-1, -1) when the matrix has more than 2 dimensions
Definition mat.hpp:2425
cv::Mat::convertTo
void convertTo(OutputArray m, int rtype, double alpha=1, double beta=0) const
Converts an array to another data type with optional scaling.
cv::Point_< int >
cv::Rect_
Template class for 2D rectangles.
Definition types.hpp:447
cv::Rect_::size
Size_< _Tp > size() const
size (width, height) of the rectangle
cv::Scalar_< double >
cv::Size_
Template class for specifying the size of an image or rectangle.
Definition types.hpp:338
cv::Size_::height
_Tp height
the height
Definition types.hpp:366
cv::Size_::width
_Tp width
the width
Definition types.hpp:365
cv::split
void split(const Mat &src, Mat *mvbegin)
Divides a multi-channel array into several single-channel arrays.
cv::mulSpectrums
void mulSpectrums(InputArray a, InputArray b, OutputArray c, int flags, bool conjB=false)
Performs the per-element multiplication of two Fourier spectrums.
cv::divide
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.
cv::sum
Scalar sum(InputArray src)
Calculates the sum of array elements.
cv::merge
void merge(const Mat *mv, size_t count, OutputArray dst)
Creates one multi-channel array out of several single-channel ones.
cv::idft
void idft(InputArray src, OutputArray dst, int flags=0, int nonzeroRows=0)
Calculates the inverse Discrete Fourier Transform of a 1D or 2D array.
cv::dft
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.
cv::String
std::string String
Definition cvstd.hpp:151
CV_8U
#define CV_8U
Definition interface.h:76
CV_32F
#define CV_32F
Definition interface.h:81
cv::datasets::circle
@ circle
Definition gr_skig.hpp:62
highgui.hpp
main
int main(int argc, char *argv[])
Definition highgui_qt.cpp:3
imgcodecs.hpp
imgproc.hpp
cv
Definition core.hpp:107
std
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.

# References

- Image Deblurring in Matlab - Image Deblurring in Matlab
- SmartDeblur - SmartDeblur site