您的位置: 首页  >  文章  >  【图像修复】浅析：图像修复中的TV模型

【图像修复】浅析：图像修复中的TV模型

分类: 文章 • 2023-03-30 16:48:08

转自：https://blog.****.net/hujingshuang/article/details/44257179，感谢原创作者！

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

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

一、TV模型介绍

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型                      ①

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

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

化简得到最终的表达式：

【图像修复】浅析：图像修复中的TV模型

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

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

matlab源码实现：

[plain] view plain copy
  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:m  
  9.     for j=1:n  
  10.         if mask(i,j)==0  
  11.            img(i,j)=0;   
  12.         end  
  13.     end  
  14. end  
  15. imshow(img,[]);     %合成的需要修复的图像  
  16.   
  17. lambda=0.2;  
  18. a=0.5;%避免分母为0  
  19. 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.             end  
  43.         end  
  44.     end  
  45.     img=imgn;   
  46. end  
  47. figure;  
  48. imshow(img)  
opencv源码实现：

[cpp] view plain copy
  1. #include <iostream>  
  2. #include <stdlib.h>  
  3. #include <cv.h>  
  4. #include <math.h>  
  5. #include <opencv2/highgui/highgui.hpp>  
  6. #include <opencv2/core/core.hpp>  
  7. #include <opencv2/imgproc/imgproc.hpp>  
  8.   
  9. using namespace cv;  
  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.     int i, 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.     float UO, UN, UW, US, UE, UNE, UNW, USW, USE;  
  44.     float Un, Ue, Uw, Us;  
  45.     float Wn, We, Ww, Ws;  
  46.     float Hon, Hoe, How, Hos;  
  47.     float Hoo;  
  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模型，在速度和效果上都比理想，此处不深入探讨。

【图像修复】浅析：图像修复中的TV模型【图像修复】浅析：图像修复中的TV模型

本文所需的实验材料都已打包上传，点击此处下载：http://download.****.net/detail/hujingshuang/8500553

如有疑问请加群：数字图像处理群：389011389

特别注意：

1、本博客例程仅做学习交流用，切勿用于商业用途。

2、欢迎交流，转载请注明出处：http://blog.****.net/hujingshuang/article/details/44257179


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

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

一、TV模型介绍

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型                      ①

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

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

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

【图像修复】浅析：图像修复中的TV模型

化简得到最终的表达式：

【图像修复】浅析：图像修复中的TV模型

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

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

matlab源码实现：

[plain] view plain copy
  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:m  
  9.     for j=1:n  
  10.         if mask(i,j)==0  
  11.            img(i,j)=0;   
  12.         end  
  13.     end  
  14. end  
  15. imshow(img,[]);     %合成的需要修复的图像  
  16.   
  17. lambda=0.2;  
  18. a=0.5;%避免分母为0  
  19. 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.             end  
  43.         end  
  44.     end  
  45.     img=imgn;   
  46. end  
  47. figure;  
  48. imshow(img)  
opencv源码实现：

[cpp] view plain copy
  1. #include <iostream>  
  2. #include <stdlib.h>  
  3. #include <cv.h>  
  4. #include <math.h>  
  5. #include <opencv2/highgui/highgui.hpp>  
  6. #include <opencv2/core/core.hpp>  
  7. #include <opencv2/imgproc/imgproc.hpp>  
  8.   
  9. using namespace cv;  
  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.     int i, 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.     float UO, UN, UW, US, UE, UNE, UNW, USW, USE;  
  44.     float Un, Ue, Uw, Us;  
  45.     float Wn, We, Ww, Ws;  
  46.     float Hon, Hoe, How, Hos;  
  47.     float Hoo;  
  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模型，在速度和效果上都比理想，此处不深入探讨。

【图像修复】浅析：图像修复中的TV模型【图像修复】浅析：图像修复中的TV模型

本文所需的实验材料都已打包上传，点击此处下载：http://download.****.net/detail/hujingshuang/8500553

如有疑问请加群：数字图像处理群：389011389

特别注意：

1、本博客例程仅做学习交流用，切勿用于商业用途。

2、欢迎交流，转载请注明出处：http://blog.****.net/hujingshuang/article/details/44257179

相关推荐