今天要整理记录的内容是，关于散焦图像的去模糊滤波，这部分内容参考学习于OpenCV官方文档。
首先我们需要先了解下什么是散焦图像。我们观察下面这两张图像
当我们看到这两张图像的时候，给我们的第一感觉是很模糊，就好像得了近视眼但是忘了带眼镜在看东西的感觉，虽然我不近视。。。也不知道忘带眼镜是什么感觉。。。这只是一个比方。。。
这种情况，其实就是拍摄物体的时候，没有把拍摄主体对上焦，也就叫做散焦图像。正常来说，当像距等于焦距时，一个物点经过透镜在成像面所成的像也应该是一个点，但是如果成像面与透镜之间的距离发生变化了，就导致光线在尚未汇聚到一点时就已经投射在成像面上了，此时形成的像就不再是一个点，而是一块光斑。
也就是说，当物体成像时，如果不是物点-像点映射，而是物点-像斑映射的时候，就会导致出现散焦图像，图像将会变得模糊、难以观看，也就是俗称的 “对不上焦”。上面的散焦图像中，第一张是我用手机拍摄的时候手动调节焦距使之对焦失败的结果，第二张则是从网上寻找的散焦图像。
应该说，这种散焦图像的信息，是受到很大程度的损失，我们很难将它完全得清晰化，因为它本身在拍摄的时候就已经丢失了很多信息。但是我们可以通过频域滤波，来对其进行去模糊操作，获得较为清晰的恢复图像。
首先，我们假设散焦图像的频谱为S，而原始的未失真图像频谱为U，那么我们可以把散焦图像频谱表示为：
其中 H 为点扩散函数PSF，N 为加性噪声。
点扩散函数（PSF——point spread function），对光学系统来讲，输入物为一点光源时其输出像的光场分布，称为点扩散函数，也称点扩展函数。在数学上点光源可用δ函数（点脉冲）代表，输出像的光场分布叫做脉冲响应，所以点扩散函数也就是光学系统的脉冲响应函数。上面提到的物点-像点映射，其实就是一种最理想的点扩散函数，然而在散焦图像中这种理想状态是不存在的。在散焦图像中，存在的PSF一般来说是物点-像斑映射，尤其是圆形PSF点扩散函数非常适合用来描述这种散焦图像的情况。
我们在频域中对散焦图像进行去模糊滤波，得到的并非百分百清晰的原图像，而仅仅只是通过模糊图像来估计、重建的图像频谱，我们将估计得到的图像频谱用 U’ 来表示，可以用下式表示：
其中 Hw 是在频域中使用的恢复滤波器，一般是维纳滤波器。维纳滤波器的简单表示如下：
其中，H是我们之前的PSF点扩散函数，SNR则是点扩散函数的信噪比。
所以通过以上公式我们可以得知，如果我们想要得到一张散焦图像的去模糊图像，则需要先定义一个合适的PSF点扩散函数，然后通过这个PSF点扩散函数来计算出一个维纳滤波器。再使用维纳滤波器在频域中对散焦图像进行去模糊滤波，最终得到一张相比之下较为清晰的恢复图像。
所以这一系列步骤的前提是：定义合适的PSF点扩散函数。
前面说到，对于散焦图像，圆形PSF点扩散函数是非常合适的，那么我们需要的就是人为规定一个半径R和信噪比SNR来确定一个圆形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];
}
这段代码中，输出的就是点扩散函数在空间域的表示，也就是一个圆形光斑。R是圆形光斑的半径，filterSize是要进行过滤的散焦图像的尺寸。得到的圆形PSF点扩散函数如下（规定不同参数能得到不同的PSF）：
得到PSF点扩散函数之后，我们就要把PSF代入上面的公式中去，计算出维纳滤波器，演示代码如下：
void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double snr)
{
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 += 1.0 / double(snr);
divide(planes[0], denom, output_G);
}
在这段代码中，有一个fftshift(input_h_PSF, h_PSF_shifted)函数，这个函数所起的作用，是将PSF在空间域中的对角区域互调，也就是象限互换，以保证将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);
}
重新排列象限后的PSF如下，可以看到原本中心的圆形区域分成四块分别位于四个对角：
那么到这里，我们就得到了维纳滤波器，接下来就可以通过维纳滤波器来对散焦图像进行去模糊滤波。注意我们需要先把散焦图像也变换到频域中，在与维纳滤波器频谱进行逐元素相乘，最后取相乘结果的实部，就得到了滤波后的恢复图像频谱了。代码演示如下：
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];
}
那么到这里，就分别定义好了PSF点扩散函数的获取、维纳滤波器的获取、散焦图像频域滤波的主要函数，下面就可以通过这几个函数进行散焦图像的去模糊滤波了。效果如下：
上图中，左边是散焦图像，右边是滤波后的恢复图像。可以看到，这两张图像都是模糊的，但是经过滤波后的图像相比起原图像还是有一定程度的改善，至少有些字能够看清了。
但很明显的是，恢复图像由于在频域进行滤波时选用的频域滤波器具有陡峭变化的缘故，产生了振铃效应，所以我们能够看到图像中有一些圆环状的波纹，在图像边缘处的波纹更加明显。那么我们想要减轻振铃效应带来的影响，就得从振铃效应产生的根源入手，也就是说，我们需要减少选用频域滤波器的陡峭变化，让滤波器的过渡变得更加平滑。
在OpenCV官方文档中，提出了一种削弱振铃效应的方法，那就是通过使恢复滤波器的边缘逐渐变细，从而减少恢复图像中的振铃效果，代码如下:
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);
}
我们再来看看使用了削弱振铃效果的滤波结果如何：
可以看到现在的恢复图像的清晰度有了些许提升，已经能够勉强的阅读上面的内容了。
完整的散焦图像去模糊滤波代码如下：
#include<opencv.hpp>
#include<stdlib.h>
#include<string>
#include<vector>
using namespace std;
using namespace cv;
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);
void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma = 5.0, double beta = 5.2);
int main()
{
Mat imgIn = imread("散焦5.jpg",0);
resize(imgIn, imgIn, Size(600, 800));
Mat imgOut;
//截取处理区域，对宽高取偶数
Rect roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
Mat Hw, h;
int R = 5; //圆形点扩散函数PSF的半径
int SNR =35; //圆形点扩散函数PSF的信噪比
//创建圆形点扩散函数PSF
calcPSF(h, roi.size(), R);
//创建维纳滤波器
calcWnrFilter(h, Hw, SNR);
//削弱振铃效应
edgetaper(Hw, Hw);
//进行滤波
filter2DFreq(imgIn(roi), imgOut, Hw);
//转换类型并归一化
imgOut.convertTo(imgOut, CV_8UC1);
normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
//namedWindow("input", WINDOW_FREERATIO);
imshow("input", imgIn);
//namedWindow("output", WINDOW_FREERATIO);
imshow("output", imgOut);
waitKey(0);
return 0;
}
//对于散焦图像去模糊时使用点扩散函数，如果针对运动图像去模糊可以使用带有角度和长度的点扩散函数，角度即是物体运动方向、长度取决于运动速度
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 calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double snr)
{
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 += 1.0 / double(snr);
divide(planes[0], denom, output_G);
}
//重新排列傅立叶图像的象限，使原点位于图像中心
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 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);
}
总结：
（1）对于散焦图像而言，因为当它拍摄的时候就已经失去了过多的信息，导致这种图像所包含的信息是非常不完整的，也就导致如果我们想要对散焦图像进行清晰化有着非常大的困难。
（2）对于这样的模糊图像，我们可以尝试通过频域滤波在一定程度上恢复图像信息，但远远无法与未失真图像相媲美。
（3）这种通过频域滤波的方法，需要我们寻找到一个最合适的PSF点扩散函数对应的半径R和信噪比SNR，其中半径R的影响权重更大，可以优先调节R的值。这个调整参数的过程相当考验我们的耐心。
（4）对于频域滤波器的选择和处理也很重要，要尽量避免振铃效应对图像结果的影响，尽量让滤波器的过渡更加平滑。
好的，那本次笔记到此结束，谢谢阅读~
PS：本人的注释比较杂，既有自己的心得体会也有网上查阅资料时摘抄下的知识内容，所以如有雷同，纯属我向前辈学习的致敬，如果有前辈觉得我的笔记内容侵犯了您的知识产权，请和我联系，我会将涉及到的博文内容删除，谢谢！