图像处理：其他——失焦去模糊滤镜 OpenCV v4.8.0

上一个教程 ：使用距离变换和分水岭算法进行图像分割

下一个教程 ：动态模糊滤波器

原作者	Karpushin Vladislav
兼容性	OpenCV >= 3.0

目标

在本教程中，你将学习

• 什么是图像劣化模型
• 什么是失焦图像的 PSF
• 如何还原模糊图像
• 什么是维纳滤波器(Wiener Filter)

理论知识

注释
本解释基于书籍 [98] 和 [311]。此外，您还可以参考 Matlab 的教程 Image Deblurring in Matlab 和文章 SmartDeblur。 本页上的失焦图像是一幅真实世界的图像。离焦是由相机光学器件手动实现的。

什么是图像劣化模型？

以下是以频域表示的图像劣化数学模型：
$$S=H\cdot U+N$$
其中 S 是模糊（降级）图像的频谱，U 是原始真实（未降级）图像的频谱，H 是点扩散函数（PSF）的频率响应，N 是加性噪声的频谱。

圆形 PSF 是焦外失真的良好近似值。这种 PSF 只需一个参数–半径 R。

圆点扩散函数

如何还原模糊图像？

还原（去毛刺）的目的是获得原始图像的估计值。频域复原公式为
$$U^{\prime }=H_w\cdot S$$
其中 U′ 是原始图像 U 的估计频谱，Hw 是修复滤波器，例如维纳滤波器。

什么是维纳滤波器？

维纳滤波器是一种修复模糊图像的方法。假设 PSF 是一个真实的对称信号，原始真实图像的功率谱和噪声都不已知，那么简化的维纳公式为

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

其中 SNR 为信噪比。

因此，要通过维纳滤波器恢复焦外图像，需要知道圆形 PSF 的 SNR 和 R。

源代码

您可以在 OpenCV 源代码库的 samples/cpp/tutorial_code/ImgProc/out_of_focus_deblur_filter/out_of_focus_deblur_filter.cpp 中找到源代码。

#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()) //检查图片是否已加载
 {
 cout << "ERROR : Image cannot be loaded..!!" << endl;
 return -1;
 }
 Mat imgOut;
 // 只需要处理偶数图像
 Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2)；
 //计算 Hw（开始）
 Mat Hw, h；
 calcPSF(h, roi.size(), R)；
 calcWnrFilter(h, Hw, 1.0 / double(snr))；
 //Hw 计算（停止）
 // 滤波（开始）
 filter2DFreq(imgIn(roi), imgOut, Hw)；
 // 滤波（停止）
 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);
}

说明

失焦图像恢复算法包括 PSF 生成、维纳滤波器生成和频域模糊图像滤波：

 // 只需处理偶数图像
 Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2)；
 //计	算 Hw（开始）
 Mat Hw, h；
 calcPSF(h, roi.size(), R)；
 calcWnrFilter(h, Hw, 1.0 / double(snr))；
 //Hw 计算（停止）
 // 滤波（开始）
 filter2DFreq(imgIn(roi), imgOut, Hw)；
 // 滤波（停止）

函数 calcPSF() 根据输入参数半径 R 形成圆形 PSF：

void calcPSF(Mat& outputImg, Size filterSize, int R)
{
 Mat h(filterSize, CV_32F, Scalar(0))；
 Point point（filterSize.宽/2，filterSize.高/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]；
}

结果

下面是真正的焦外图像：

失焦图像

下面是使用 R = 53 和 SNR = 5200 参数计算得出的结果：

恢复（去模糊）后的图像

我们使用了维纳滤波器，并手动选择了 R 值和 SNR 值，以获得最佳视觉效果。我们可以看到，结果并不完美，但它为我们提供了图像内容的提示。虽然有些困难，但文字还是可以阅读的。

注意事项
参数 R 是最重要的。因此应首先调整 R，然后再调整 SNR。
有时您会在修复后的图像中观察到振铃效应。有几种方法可以减少这种效果。例如，可以将输入图像的边缘变细。

您也可以在 YouTube 上找到这方面的快速视频演示。

参考资料

• Matlab 中的图像去模糊 - Matlab 中的图像去模糊
• SmartDeblur - SmartDeblur 网站