【圖像修複】淺析：圖像修複中的TV模型

前言：圖像修複是一項非常有意義的研究工作，比如我們生活中的照片被污染，再比如名貴字畫、國家文物壁畫等珍貴物品被破壞，這些都需要圖像修複工作來完成。

簡介：整體變分（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源碼實作：

img=double(imread('lena.jpg'));
mask=imread('mask.jpg');
a1=find(mask>127);
b1=find(mask<=127);
mask(a1)=0;
mask(b1)=255;
[m n]=size(img);
for i=1:m
    for j=1:n
        if mask(i,j)==0
           img(i,j)=0; 
        end
    end
end
imshow(img,[]);     %合成的需要修複的圖像

lambda=0.2;
a=0.5;%避免分母為0
imgn=img;
for l=1:1500         %疊代次數
    for i=2:m-1
        for j=2:n-1
            if mask(i,j)==0     %如果目前像素是被污染的像素，則進行處理
                Un=sqrt((img(i,j)-img(i-1,j))^2+((img(i-1,j-1)-img(i-1,j+1))/2)^2);
                Ue=sqrt((img(i,j)-img(i,j+1))^2+((img(i-1,j+1)-img(i+1,j+1))/2)^2);
                Uw=sqrt((img(i,j)-img(i,j-1))^2+((img(i-1,j-1)-img(i+1,j-1))/2)^2);
                Us=sqrt((img(i,j)-img(i+1,j))^2+((img(i+1,j-1)-img(i+1,j+1))/2)^2);

                Wn=1/sqrt(Un^2+a^2);
                We=1/sqrt(Ue^2+a^2);
                Ww=1/sqrt(Uw^2+a^2);
                Ws=1/sqrt(Us^2+a^2);

                Hon=Wn/((Wn+We+Ww+Ws)+lambda);
                Hoe=We/((Wn+We+Ww+Ws)+lambda);
                How=Ww/((Wn+We+Ww+Ws)+lambda);
                Hos=Ws/((Wn+We+Ww+Ws)+lambda);

                Hoo=lambda/((Wn+We+Ww+Ws)+lambda);
                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);
                imgn(i,j)= value;
            end
        end
    end
    img=imgn; 
end
figure;
imshow(img)
           

opencv源碼實作：

#include <iostream>
#include <stdlib.h>
#include <cv.h>
#include <math.h>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>

using namespace cv;

int main(void)
{
	//讀取原始圖像及掩模圖像
	IplImage *src_uint8 = cvLoadImage("src.jpg", CV_LOAD_IMAGE_GRAYSCALE);
	IplImage *mask = cvLoadImage("mask.jpg", CV_LOAD_IMAGE_GRAYSCALE);
	//合成需要修複的圖像
	int M = mask->height;
	int N = mask->width;
	int i, j;
	CvMat *src = cvCreateMat(M, N, CV_32FC1);//存放浮點圖像
	cvConvert(src_uint8, src);
	for (i = 0; i < M; i++)
	{
		for (j = 0; j < N; j++)
		{
			if ((mask->imageData + i * mask->widthStep)[j] < 0)//了解此處判别條件，根據情況自行更改
			{
				((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] = 0.0;
			}
			if (((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] < 0)
			{
				((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j] += 256.0;
			}
		}
	}
	cvConvert(src, src_uint8);
	cvShowImage("需要修複的圖像", src_uint8);
	cvWaitKey(0);

	double t = getTickCount();//目前滴答數
	float lambda = 0.2;
	float delta = 0.5;
	float UO, UN, UW, US, UE, UNE, UNW, USW, USE;
	float Un, Ue, Uw, Us;
	float Wn, We, Ww, Ws;
	float Hon, Hoe, How, Hos;
	float Hoo;
	int iteration = 500;
	while(iteration)
	{
		for (i = 1; i < M - 1; i++)
		{
			for (j = 1; j < N - 1; j++)
			{
				if (((char *)(mask->imageData + i * mask->widthStep))[j] < 0)//壞損區
				{
					UO = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j];
					UN = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i-1)))[j];
					US = ((float*)(void*)(src->data.ptr + (size_t)src->step*(i+1)))[j];
					UE = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j+1];
					UW = ((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j-1];

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

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

					Wn = 1.0/sqrt(Un * Un + delta * delta);
					We = 1.0/sqrt(Ue * Ue + delta * delta);
					Ww = 1.0/sqrt(Uw * Uw + delta * delta);
					Ws = 1.0/sqrt(Us * Us + delta * delta);

					Hon = Wn/(Wn+We+Ww+Ws+lambda);
					Hoe = We/(Wn+We+Ww+Ws+lambda);
					How = Ww/(Wn+We+Ww+Ws+lambda);
					Hos = Ws/(Wn+We+Ww+Ws+lambda);

					Hoo = lambda/(Wn+We+Ww+Ws+lambda);
					((float*)(void*)(src->data.ptr + (size_t)src->step*i))[j]=(Hon*UN+Hoe*UE+How*UW+Hos*US+Hoo*UO);
				}
			}
		}
		iteration--;
	}
	cvConvert(src, src_uint8);
	t = ((double)getTickCount() - t)/getTickFrequency();
	printf("算法用時：%f秒\n", t);
	cvShowImage("修複結果", src_uint8);
	cvWaitKey(0);
}
           

由于疊代次數和浮點數的運算，使得算法時間較長，效果如下，仔細觀察可以看出仍有細節處修複效果不是很理想。在TV模型之後，又出現了許多改進的TV模型，在速度和效果上都比理想，此處不深入探讨。

本文所需的實驗材料都已打包上傳，點選此處下載下傳：http://download.csdn.net/detail/hujingshuang/8500553

如有疑問請加群：數字圖像處理群：389011389

特别注意：

1、本部落格例程僅做學習交流用，切勿用于商業用途。

2、歡迎交流，轉載請注明出處：http://blog.csdn.net/hujingshuang/article/details/44257179