前言：图像修复是一项非常有意义的研究工作，比如我们生活中的照片被污染，再比如名贵字画、国家文物壁画等珍贵物品被破坏，这些都需要图像修复工作来完成。

简介：整体变分(Total Variation)的方法最早是用来对受到噪声污染的图像进行降噪的，在这方面的应用最早是由L.Rudin和S.Osher等人在1992年提出的，2002年Chan等人把TV模型推广到图像修补中，并提出了基于TV模型的图像修补方法，同时说明了TV修补模型的缺点，进一步提出了CDD修补模型(curvature driven diffusions)，此修补模型改正了TV修补模型的缺陷，对图像的修补具有很好的效果。

一、TV模型介绍

如图所示：D区域是被污染区(待修复)，E是D的邻域

下面直接给出TV模型的数学公式：

                      ①

其中：u是图像中的像素点，λ为设定的参数

在该模型基础上，考虑到噪声的影响，边界E区域产生的噪声不能超过一定的范围；根据最佳猜测和贝叶斯理论，要求图像u在满足约束条件下使它的能量泛函最小，因此约束条件记做：公式②。根据拉格朗日乘数法，将①②方程转化成为一个求极值的方程，对其求导数并令其等于0，可得到如下方程：

其中：div代表散度(关于图像中的散度解释，可见此处：在图像处理中，散度 div 具体的作用是什么？)

由于图像是离散的数值，故可看做如下构成。其中：O为污染点，邻域B=(N,S,W,E)，半像素邻域B‘ =(n,s,w,e)。

因此，离散化后可得到表达式：

化简得到最终的表达式：

其中：λe(O)为中心O处的λ参数，与上λe一致；uo为O点修复后的像素，另一个为O点修复前的原始像素。将上式迭代，知道达到较好的修复效果。

到此，TV模型的理论推导已完成。接下来就是要编程实现其功能。

matlab源码实现：

1 img=double(imread(‘lena.jpg‘));2 mask=imread(‘mask.jpg‘);3 a1=find(mask>127);4 b1=find(mask<=127);5 mask(a1)=0;6 mask(b1)=255;7 [m n]=size(img);8 for i=1:m9 for j=1:n10 if mask(i,j)==0

11 img(i,j)=0;12 end13 end14 end15 imshow(img,[]); %合成的需要修复的图像16

17 lambda=0.2;18 a=0.5;%避免分母为019 imgn=img;20 for l=1:1500 %迭代次数21 for i=2:m-1

22 for j=2:n-1

23 if mask(i,j)==0 %如果当前像素是被污染的像素，则进行处理24 Un=sqrt((img(i,j)-img(i-1,j))^2+((img(i-1,j-1)-img(i-1,j+1))/2)^2);25 Ue=sqrt((img(i,j)-img(i,j+1))^2+((img(i-1,j+1)-img(i+1,j+1))/2)^2);26 Uw=sqrt((img(i,j)-img(i,j-1))^2+((img(i-1,j-1)-img(i+1,j-1))/2)^2);27 Us=sqrt((img(i,j)-img(i+1,j))^2+((img(i+1,j-1)-img(i+1,j+1))/2)^2);28

29 Wn=1/sqrt(Un^2+a^2);30 We=1/sqrt(Ue^2+a^2);31 Ww=1/sqrt(Uw^2+a^2);32 Ws=1/sqrt(Us^2+a^2);33

34 Hon=Wn/((Wn+We+Ww+Ws)+lambda);35 Hoe=We/((Wn+We+Ww+Ws)+lambda);36 How=Ww/((Wn+We+Ww+Ws)+lambda);37 Hos=Ws/((Wn+We+Ww+Ws)+lambda);38

39 Hoo=lambda/((Wn+We+Ww+Ws)+lambda);40 value = Hon*img(i-1,j)+Hoe*img(i,j+1)+How*img(i,j-1)+Hos*img(i+1,j)+Hoo*img(i,j);41 imgn(i,j)=value;42 end43 end44 end45 img=imgn;46 end47 figure;48 imshow(img)

opencv源码实现：

1 #include

2 #include

3 #include

4 #include

5 #include

6 #include

7 #include

8

9 using namespacecv;10

11 int main(void)12 {13 //读取原始图像及掩模图像

14 IplImage *src_uint8 = cvLoadImage("src.jpg", CV_LOAD_IMAGE_GRAYSCALE);15 IplImage *mask = cvLoadImage("mask.jpg", CV_LOAD_IMAGE_GRAYSCALE);16 //合成需要修复的图像

17 int M = mask->height;18 int N = mask->width;19 inti, j;20 CvMat *src = cvCreateMat(M, N, CV_32FC1);//存放浮点图像

21 cvConvert(src_uint8, src);22 for (i = 0; i < M; i++)23 {24 for (j = 0; j < N; j++)25 {26 if ((mask->imageData + i * mask->widthStep)[j] < 0)//理解此处判别条件，根据情况自行更改

27 {28 ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] = 0.0;29 }30 if (((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] < 0)31 {32 ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] += 256.0;33 }34 }35 }36 cvConvert(src, src_uint8);37 cvShowImage("需要修复的图像", src_uint8);38 cvWaitKey(0);39

40 double t = getTickCount();//当前滴答数

41 float lambda = 0.2;42 float delta = 0.5;43 floatUO, UN, UW, US, UE, UNE, UNW, USW, USE;44 floatUn, Ue, Uw, Us;45 floatWn, We, Ww, Ws;46 floatHon, Hoe, How, Hos;47 floatHoo;48 int iteration = 500;49 while(iteration)50 {51 for (i = 1; i < M - 1; i++)52 {53 for (j = 1; j < N - 1; j++)54 {55 if (((char *)(mask->imageData + i * mask->widthStep))[j] < 0)//坏损区

56 {57 UO = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j];58 UN = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i-1)))[j];59 US = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i+1)))[j];60 UE = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j+1];61 UW = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j-1];62

63 UNE = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i-1)))[j+1];64 UNW = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i-1)))[j-1];65 USE = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i+1)))[j+1];66 USW = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i+1)))[j-1];67

68 Un = sqrt((UO - UN) * (UO - UN) + ((UNW - UNE) / 2.0) * ((UNW - UNE) / 2.0));69 Ue = sqrt((UO - UE) * (UO - UE) + ((UNE - USE) / 2.0) * ((UNE - USE) / 2.0));70 Uw = sqrt((UO - UW) * (UO - UW) + ((UNW - USW) / 2.0) * ((UNW - USW) / 2.0));71 Us = sqrt((UO - US) * (UO - US) + ((USW - USE) / 2.0) * ((USW - USE) / 2.0));72

73 Wn = 1.0/sqrt(Un * Un + delta *delta);74 We = 1.0/sqrt(Ue * Ue + delta *delta);75 Ww = 1.0/sqrt(Uw * Uw + delta *delta);76 Ws = 1.0/sqrt(Us * Us + delta *delta);77

78 Hon = Wn/(Wn+We+Ww+Ws+lambda);79 Hoe = We/(Wn+We+Ww+Ws+lambda);80 How = Ww/(Wn+We+Ww+Ws+lambda);81 Hos = Ws/(Wn+We+Ww+Ws+lambda);82

83 Hoo = lambda/(Wn+We+Ww+Ws+lambda);84 ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j]=(Hon*UN+Hoe*UE+How*UW+Hos*US+Hoo*UO);85 }86 }87 }88 iteration--;89 }90 cvConvert(src, src_uint8);91 t = ((double)getTickCount() - t)/getTickFrequency();92 printf("算法用时：%f秒\n", t);93 cvShowImage("修复结果", src_uint8);94 cvWaitKey(0);95 }

由于迭代次数和浮点数的运算，使得算法时间较长，效果如下，仔细观察可以看出仍有细节处修复效果不是很理想。在TV模型之后，又出现了许多改进的TV模型，在速度和效果上都比理想，此处不深入探讨。

原文：http://www.cnblogs.com/hxjbc/p/6675901.html