2006-TOG - Removing camera shake from a single photograph

项目地址：https://cs.nyu.edu/~fergus/research/deblur.html

• 多尺度
• 变分贝叶斯
• 高斯混合模型
• 用户选取区域及模糊核方向
• 泊松重构中间图像
• Richardson-Lucy (RL)算法恢复最终图像（优势：速度快）

假设为一致性模糊（卷积操作$\otimes$的边界条件为非周期 non-periodic）：
$$\mathbf{B}=\mathbf{K}\otimes \mathbf{L}+\mathbf{N}$$

In order to estimate the latent image from such limited measurements, it is essential to have some notion of which images are apriori more likely.

先验的作用：从有限的测量值B估计L（K未知），有很多解，有必要知道那种解更有可能。

研究表明，真实场景图像虽然像素分布变化很大，但是梯度分布服从重尾分布。

the distribution of gradients has most of its mass on small values but gives signicantly more probability to large values than a aussian distribution.

重尾分布相比高斯分布在梯度值大的区域具有更高的概率，如果使用高斯分布则恢复出的图像过于平滑。

The estimation process is performed in a coarse-to-fine fashion in order to avoid local minima.

估计模糊核时通过有粗到细的方式，可以避免局部最小值，特别是针对大模糊核。

算法需要用户提供（1）模糊的图像$\mathbf{P}$，（2）模糊图像中的矩形块（局部区域$\mathbf{P}$），（3）模糊核大小的上限(以像素为单位)，以及（4）模糊核的初始方向(水平或垂直)。

由于先验是关于图像梯度的，因此算法优化需要从像素域转到梯度域。
因为卷积操作是线性操作，所以$\nabla \mathbf{P}=\nabla {\mathbf{L}}_{\mathbf{p}}\otimes \mathbf{K}$

论文使用C个分量（component）的高斯混合模型表示图像梯度先验，D个分量的混合指数模型表示模糊核先验。

论文首先通过最大后验概率MAP实现去卷积（deconvolution）, 等价于解一个正则化的最小二乘问题（a regularized-least squares problem）传送门：最小二乘法的本质是什么
但结果并未成功。一种解释是，MAP的目标函数试图最小化所有的梯度(即使是较大的梯度)，而我们期望自然图像具有一些较大的梯度。此外，作者发现MAP目标函数很容易受到局部极小值的影响。

论文采用变分贝叶斯（variational Bayesian, VB）估计一个分布$q(\mathbf{K},\nabla {\mathbf{L}}_{\mathbf{p}})$近似$p(\mathbf{K},\nabla {\mathbf{L}}_{\mathbf{p}}|\nabla \mathbf{P})$，代价函数为：
$$KL(q(\mathbf{K},\nabla {\mathbf{L}}_{\mathbf{p}},{\sigma }^{-2})||p(\mathbf{K},\nabla {\mathbf{L}}_{\mathbf{p}}|\nabla \mathbf{P}))$$
具体优化过程没看懂=_=
传送门：

1. 变分贝叶斯相关PPT
2. 变分贝叶斯初探

附:

绘制图像梯度分布代码：

f=imread('onion.png');
subplot(1,2,1);
if size(f,3) == 3
    f_gray = rgb2gray(f);
else
    f_gray = f;
end
f_gray=double(f_gray);

kx = [-1 1; 0 0];
ky = [-1 0; 1 0];

f_grad_x=conv2(f_gray,kx,'valid');%水平方向梯度
f_grad_y=conv2(f_gray,ky,'valid');%垂直方向梯度
% g=sqrt(f_grad_x.^2+f_grad_y.^2);

Ix=f_grad_x;
Value=unique(Ix(:));
Count=[hist(Ix(:),Value)]';
s=sum(Count(:));
count1=log2(Count/s);
figure,plot(Value,count1,'Color','black','LineWidth',2);
axis([-100 100 -18 0]);
title('Heavy-tailed distribution on image gradients');
xlabel('Gradient');
ylabel('Log2 probability density');
grid on;