图像修复项目《问题一》

首先，我们项目研究的时关于图像修复的快速算法，当然是基于前人基础的提出改进的算法。现在研究的一篇论文是Mingqiang Zhu的一篇论文《An Efficient Primal-Dual Hybrid Gradient Algorithm For Total Variation Image Restoration 》，因此本系列的项目都是围绕这篇论文展开的。

这篇论文主要研究的是原始对偶混合梯度算法的全变差图像复原问题。

我们知道全变分模型已经在国内有了一个好的图像修复效果，最基础的问题就是去噪。

1、背景

直接看以下最小化问题：

 

 

 

 

 

 

 

 

2、全变分离散化

那么，假设图像域为方形的，现在在这里面定义一个n*n的像素格，每个元素的索引使用（i, j），其中i=1,2,...,n，j=1,2,...,n。以这个n*n的矩阵来代表一个图像。u(i,j)表示函数u在像素（i,j）上的值。全变分定义为图像梯度幅值之和：

其中:

且符号||.||是一个欧式范式，在二维图像中指代的就是二范式。二范式的定义可以自行百度，从空间中两个点之间的距离来理解是最简单的二范式。

为了在矩阵代数中描述问题我们需要对图像矩阵u（对应加噪图像f）拉直为y向量（对应z向量），具体的方法是将二维中的矩阵元素（i,j）转换到向量中的第（j-1）*n+i个元素。也就是说第i列，第j行中的元素转换为向量中第l=（j-1）*n+i个元素。

举个例子比较好理解：

这里面的mod指的是整除的意思。

同样上面的例子，在这里应用就是考虑对应以下四种情况（注意：实际情况是u是未知图像，因此y也是未知，因此我们只能令A(l)为特定的矩阵才能满足以上条件）：

（1）l=1时，1%3!=0且1〈=9-3，此时A(1).T*y=(y(2)-y(1), y(4)-y(1)).T：

（2）l=3时，3%3==0且3〈=9-3，此时A(3).T*y=(0, y(6)-y(3)).T：

（3）l=7时，7%3！=0且7〉9-3，此时A(7).T*y=((y(8)-y(7)), 0).T：

（4）l=9时，9%3==0且9〉9-3，此时A(9).T*y=(0, 0).T：

那么，很明显，A(9)是一个零矩阵，所以如果A默认是稀疏矩阵的话就不需要考虑第四种情况啦。

现在问题来了，如果一张图很小，就只有上面例子中的3行3列的矩阵像素组成，那么就有9个A(l)，每一个A(l)都是2行9列大小的矩阵，也就是说我们现在要做的就是将这9个矩阵存储在一个大矩阵A中，A就是2*9行9列的矩阵。根据以上算法轻松制订matlab代码：

function A=DisGra(n)  
N=n^2;
B=sparse(2*N,N); % 初始化A的转置B为零矩阵
for l=1:N-n  % 前两个情况
    i=2*l-1;  % i指代A中A(l)所在的行
    if (rem(l,n)~=0)  % 如果l不被n整除
      B(i,l)=-1;
      B(i,l+1)=1;
      B(i+1,l)=-1;
      B(i+1,l+n)=1;
    else
       B(i+1,l)=-1;
       B(i+1,l+n)=1;
    end
end
for l=N-n+1:N-1  % 第三个情况
    i=2*l-1;
    B(i,l)=-1;
    B(i,l+n)=1;
end
A=B';

运行脚本如下：

tic
A = DisGra(256);
toc

得到的结果是：

Elapsed time is 25.362021 seconds when n=256

哈哈，不用我翻译吧，想想，一般的普通经典的图片如我的头像lena.jpg都是512行512列的矩阵，也就是说一般的图片n=512花费的时间太长咱们不考虑，为了方便我们取小图片n=256花费的时间都是25秒，想想半分钟计算一个矩阵也是蛮长时间的，那么以上算法有没有改进的空间呢？

 

公式（1）就可以写成以下形式：

 

3、对偶问题

对偶算法在TV图像修复中扮演着很重要的角色，尤其是在当前快速发展的数值算法中往往利用对偶算法能计算出最优解。本节中，我们要解决离散原问题ROF模型（对应公式（5）的）的对偶问题（以下问题直接取上面的例子解说，也就是n=3，N=9）。

由柯西不等式原则可知存在向量b满足：

使用以上公式（6）同理与可将y的离散TV模型写成：

其中，由于A(l)是9行两列的矩阵，那么A=[A(1),A(2),...,A(9)]是一个9行18列的矩阵；且：