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

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

目标

在本教程中，你将学习

什么是动态模糊图像的 PSF
如何还原动态模糊图像

理论

有关退化图像模型理论和维纳滤波理论，请参阅教程 “去失焦模糊滤波器”。本页仅考虑线性运动模糊失真。本页中的运动模糊图像是真实世界中的图像。模糊是由移动的主体造成的。

什么是运动模糊图像的 PSF？

线性运动模糊失真的点扩散函数（PSF）是一条线段。这种 PSF 由两个参数指定： LEN 是模糊的长度，THETA 是运动的角度。

在这里插入图片描述

线性运动模糊变形的点扩散函数 \

如何还原模糊图像？

本页使用维纳滤波器作为复原滤波器，详情请参考教程 “去失焦模糊滤波器”。为了合成运动模糊情况下的维纳滤波器，需要指定 PSF 的信噪比（SNR）、LEN 和 THETA。

源代码

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

#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 len, double theta);
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);
void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma = 5.0, double beta = 0.2);
const String keys =
"{help h usage ? | | print this message }"
"{image |input.png | input image name }"
"{LEN |125 | length of a motion }"
"{THETA |0 | angle of a motion in degrees }"
"{SNR |700 | 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 LEN = parser.get<int>("LEN");
 double THETA = parser.get<double>("THETA");
 int snr = parser.get<int>("SNR");
 string strInFileName = parser.get<String>("image");
 if (!parser.check())
 {
 parser.printErrors();
 return 0;
 }
 Mat imgIn;
 imgIn = imread(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(), LEN, THETA)；
 calcWnrFilter(h, Hw, 1.0 / double(snr))；
 //计算 Hw（停止）
 imgIn.convertTo(imgIn, CV_32F)；
 edgetaper(imgIn, imgIn)；
 // 滤波（开始）
 filter2DFreq(imgIn(roi), imgOut, Hw)；
 // 滤波（停止）
 imgOut.convertTo(imgOut, CV_8U);
 normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
 imwrite("result.jpg", imgOut);
 return 0;
}
void help()
{
 cout << "2018-08-14" << endl;
 cout << "Motion_deblur_v2" << endl;
 cout << "You will learn how to recover an image with motion blur distortion using a Wiener filter" << endl;
}
void calcPSF(Mat& outputImg, Size filterSize, int len, double theta)
{
 Mat h(filterSize, CV_32F, Scalar(0));
 Point point(filterSize.width / 2, filterSize.height / 2);
 ellipse(h, point, Size(0, cvRound(float(len) / 2.0)), 90.0 - theta, 0, 360, Scalar(255), FILLED);
 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);
}
void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma, double beta)
{
 int Nx = inputImg.cols;
 int Ny = inputImg.rows;
 Mat w1(1, Nx, CV_32F, Scalar(0));
 Mat w2(Ny, 1, CV_32F, Scalar(0));
 float* p1 = w1.ptr<float>(0);
 float* p2 = w2.ptr<float>(0);
 float dx = float(2.0 * CV_PI / Nx);
 float x = float(-CV_PI);
 for (int i = 0; i < Nx; i++)
 {
 p1[i] = float(0.5 * (tanh((x + gamma / 2) / beta) - tanh((x - gamma / 2) / beta)));
 x += dx;
 }
 float dy = float(2.0 * CV_PI / Ny);
 float y = float(-CV_PI);
 for (int i = 0; i < Ny; i++)
 {
 p2[i] = float(0.5 * (tanh((y + gamma / 2) / beta) - tanh((y - gamma / 2) / beta)));
 y += dy;
 }
 Mat w = w2 * w1;
 multiply(inputImg, w, outputImg);
}

说明

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

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

函数 calcPSF() 根据输入参数 LEN 和 THETA（单位：度）形成 PSF：

void calcPSF(Mat& outputImg, Size filterSize, int len, double theta)
{
 Mat h(filterSize, CV_32F, Scalar(0))；
 Point point（filterSize.宽/2，filterSize.高/2）；
 ellipse(h, point, Size(0, cvRound(float(len) / 2.0)), 90.0 - theta, 0, 360, Scalar(255), FILLED)；
 Scalar summa = sum(h)；
 outputImg = h / summa[0]；
}

函数 edgetaper() 将输入图像的边缘变细，以减少还原图像中的振铃效应：

void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma, double beta)
{
 int Nx = inputImg.cols；
 int Ny = inputImg.rows；
 Mat w1(1, Nx, CV_32F, Scalar(0))；
 Mat w2(Ny, 1, CV_32F, Scalar(0))；
 float* p1 = w1.ptr<float>(0)；
 float* p2 = w2.ptr<float>(0)；
 float dx = float(2.0 * CV_PI / Nx)；
 float x = float(-CV_PI)；
 for (int i = 0; i < Nx; i++)
 {
 p1[i] = float(0.5 * (tanh((x + gamma / 2) / beta) - tanh((x - gamma / 2) / beta)))；
 x += dx；
 }
 float dy = float(2.0 * CV_PI / Ny)；
 float y = float(-CV_PI)；
 for (int i = 0; i < Ny; i++)
 {
 p2[i] = float(0.5 * (tanh((y + gamma / 2) / beta) - tanh((y - gamma / 2) / beta)))；
 y += dy；
 }
 Mat w = w2 * w1；
 multiply(inputImg, w, outputImg)；
}