Goal
在本教程中，您将学习：
什么是退化图像（degradation image model）模型
失焦图像的 PSF （PSF of an out-of-focus image）是多少
如何恢复模糊的图像（restore a blurred image）
什么是维纳滤波器（Wiener filter）
Theory
笔记
解释基于书籍[95]和[303]。 另外，您可以参考 Matlab 的教程 Image Deblurring in Matlab Image Deblurring - MATLAB & Simulink (mathworks.com) 和文章 SmartDeblurRestoration of defocused and blurred images. Yuzhikov.com 。
此页面上的失焦图像是真实世界的图像。 失焦是通过相机光学系统手动实现的。
What is a degradation image model?
什么是退化图像模型？
这是频域表示(frequency domain representation)中图像退化的数学模型：
S=H?U+N
其中 S 是模糊（退化）图像的频谱，U 是原始真实（未退化）图像的频谱，H 是点扩散函数(point spread function) (PSF) 的频率响应，N 是额外噪声的频谱。
圆形 PSF 是离焦失真的良好近似。 这样的 PSF 仅由一个参数指定 - 半径 R。在这项工作中使用圆形 PSF。
Circular point spread function
?
How to restore a blurred image?
如何恢复模糊的图像？
恢复（去模糊）的目的是获得对原始图像的估计。 频域恢复公式为：
其中U'是原始图像U的估计谱，Hw是恢复滤波器，例如Wiener filter.维纳滤波器。
Opencv 实现图像的离散傅里叶变换（DFT）、卷积运算（相关滤波） - gavanwanggw - 博客园 (cnblogs.com)
?
What is the Wiener filter?
维纳滤波器是一种恢复模糊图像的方法。 假设PSF是a real and symmetric signal一个真实的对称信号，原始真实图像的功率谱和噪声是未知的，那么简化的维纳公式为：
?
其中 SNR 是信噪比。
因此，为了通过维纳滤波器恢复失焦图像，需要知道圆形 PSF 的 SNR 和 R。
其中SNR表示信噪比，因此可以基于维纳滤波恢复离焦图像，实现图像反模糊。这个过程最终重要的两个参数，分别是半径R与信噪比SNR，在反模糊图像时候，要先尝试调整R，然后再尝试调整SNR。
Source code
您可以在 OpenCV 源代码库的 samples/cpp/tutorial_code/ImgProc/out_of_focus_deblur_filter/out_of_focus_deblur_filter.cpp 中找到源代码。
/**
* @brief You will learn how to recover an out-of-focus image by Wiener filter
* 您将学习如何通过维纳滤波器恢复失焦图像
* @author Karpushin Vladislav, karpushin@ngs.ru, https://github.com/VladKarpushin
*/
#include <iostream>
#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 |53 | radius }"
"{SNR |5200 | signal to noise ratio}"
;
int main(int argc, char *argv[])
{
help();
CommandLineParser parser(argc, argv, keys);
if (parser.has("help"))
{
parser.printMessage();
return 0;
}
//圆形 PSF 是离焦失真的良好近似
int R = parser.get<int>("R");//计算圆形点扩散函数point spread function 的半径 R
int snr = parser.get<int>("SNR");//SNR 是信噪比
string strInFileName = parser.get<String>("image");//输入图像路径
if (!parser.check())
{
parser.printErrors();
return 0;
}
Mat imgIn;
imgIn = imread(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;
//! [main]
// 它只需要处理偶数行列图像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);//圆形的PSF因为只有一个半径参数R，是一个非常好的失焦畸变近似，所以算法采用圆形的PSF。 生成PSF
calcWnrFilter(h, Hw, 1.0 / double(snr));//计算Hw 维纳滤波器 nsr=1/snr
//Hw calculation (stop)
// filtering (start) // 去模糊
filter2DFreq(imgIn(roi), imgOut, Hw);
// filtering (stop)
//! [main]
// 归一化显示
imgOut.convertTo(imgOut, CV_8U);
normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
imwrite("result.jpg", imgOut);
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;
}
//! [calcPSF] 点扩散函数图
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];
}
//! [calcPSF]
//! [fftshift]快速傅里叶变换 频率对齐函数（fftshift）
//fftshift的作用正是让正半轴部分和负半轴部分的图像分别关于各自的中心对称，
//因为直接用fft得出的数据与频率不是对应的，fftshift可以纠正过来。
//具体实现，每一行后半部分移到前半部分，每一列后半部分移到前半部分。
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);
}
//! [fftshift]
//! [filter2DFreq] 频域去模糊滤波，滤波后输出 时域图
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);//计算频域图complexI
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);//对于两张频谱图中每个元素的乘法 U′=Hw?S U'是原始图像U的估计谱，Hw是恢复滤波器
idft(complexIH, complexIH);//逆向傅里叶变换
split(complexIH, planes);//
outputImg = planes[0];//去模糊滤波后时域图
}
//! [filter2DFreq]
//! [calcWnrFilter] 计算维纳滤波器
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;//|H|^2+1/nsr
divide(planes[0], denom, output_G);//计算维纳滤波器Hw -> output_G H/(|H|^2+1/nsr)
}
//! [calcWnrFilter]
Explanation
失焦图像恢复算法包括 PSF 生成、Wiener 滤波器生成和在频域中过滤模糊图像：
??? //它只需要处理偶数图像 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)
函数 calcPSF() 根据输入参数半径 R 形成一个圆形 PSF：
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];
}
函数 calcWnrFilter() 根据上述公式合成简化的维纳滤波器 Hw：
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);
}
函数 fftshift() 重新排列 PSF。 此代码刚刚从教程离散傅里叶变换中复制：
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);
}
函数 filter2DFreq() 过滤频域中的模糊图像：
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
下面你可以看到真正的失焦图像：
?并且使用 R = 53 和 SNR = 5200 参数计算了以下结果：
?
使用了 Wiener 滤波器，并且手动选择了 R 和 SNR 的值以提供最佳的视觉效果。 我们可以看到结果并不完美，但它为我们提供了图像内容的提示。 有一些困难，文本是可读的。
笔记
参数 R 是最重要的。 所以你应该先调整R，然后再调整SNR。
有时您可以在恢复的图像中观察到振铃效应。 可以通过几种方法来减少这种影响。 例如，您可以使输入图像边缘变细。
References
- Image Deblurring in Matlab?- Image Deblurring in Matlab
- SmartDeblur?- SmartDeblur site
OpenCV 频域-时域的傅里叶变换及逆变换 C++_Dwyane05的博客-CSDN博客_c++ opencv 傅里叶变换
opencv的频域滤波 - 凤凰_1 - 博客园 (cnblogs.com)
快速傅里叶变换C++完整实现（包括FFT、IFFT & FFTShift） - SegmentFault 思否
您还可以在 YouTube 上找到一个快速演示视频。
https://youtu.be/0bEcE4B0XP4