散焦模糊去除滤波

目标

在本教程中，你将学到：

• 什么是图像退化模型（degradation image model）；
• 什么是散焦图像的点扩散函数（PSF，point spread function）；
• 如何恢复一幅模糊图像；
• 什么是维纳滤波器（Wiener filter）。

原理

Note

本说明基于书籍 Gonzalez1987 与 gruzman。你也可以参阅 Matlab 的教程 [Image Deblurring in Matlab] 以及文章 [SmartDeblur]。

Note

本页中的散焦图像是一幅真实世界图像。散焦效果是通过相机光学系统手动实现的。

什么是图像退化模型？

下面是频域表示下的图像退化数学模型：

$S = H\cdot U + N$

其中 $S$ 是模糊（退化）图像的频谱，$U$ 是原始真实（未退化）图像的频谱，$H$ 是点扩散函数（PSF）的频率响应，$N$ 是加性噪声的频谱。

圆形 PSF 是对散焦失真的一种良好近似。这样的 PSF 只由一个参数——半径 $R$ 指定。本工作使用的就是圆形 PSF。

如何恢复一幅模糊图像？

恢复（去模糊）的目标是获得原始图像的一个估计。频域中的恢复公式为：

$U' = H_w\cdot S$

其中 $U'$ 是原始图像 $U$ 的估计的频谱，$H_w$ 是恢复滤波器，例如维纳滤波器。

什么是维纳滤波器？

维纳滤波器是恢复模糊图像的一种方法。假设 PSF 是实对称信号，且原始真实图像与噪声的功率谱均未知，那么简化后的维纳公式为：

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

其中 $SNR$ 是信噪比（signal-to-noise ratio）。

因此，为了用维纳滤波器恢复一幅散焦图像，需要知道圆形 PSF 的 $SNR$ 和 $R$。

源代码

你可以在 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/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);
}

说明

一个散焦图像恢复算法由 PSF 生成、维纳滤波器生成以及在频域中对模糊图像进行滤波三部分组成：

// 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() 按照上面描述的公式合成简化的维纳滤波器 $H_w$：

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 进行重新排列。这段代码是从离散傅里叶变换（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);
}

函数 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$ 的值是手动选择的，以获得尽可能好的视觉效果。可以看到，结果并不完美，但它让我们大致看出了图像的内容。尽管有些吃力，文字已经可以辨认了。

Note

参数 $R$ 是最重要的。因此你应该先调整 $R$，然后再调整 $SNR$。

Note

有时你会在恢复后的图像中观察到振铃效应（ringing effect）。这种效应可以用几种方法来减轻，例如对输入图像的边缘做渐晕（taper）处理。

你还可以在 YouTube 上找到一个快速视频演示。

参考

• Image Deblurring in Matlab — Matlab 中的图像去模糊
• SmartDeblur — SmartDeblur 网站