www.gusucode.com > 非纹理的曲率驱动扩散的图像修复matlab源码程序 > CDDinpainting/L3_29.m

    function [B] = L3_29()
R=imread('L4_6.jpg'); %读入待修复图像
R=double(R);        %把图像值uint8型转为double型
[nrow ncol]=size(R);  %读图像大小，nrow为图像高度，ncol为图像宽度
% R_mask=imread('mask.jpg');%读模版
mask=double(R>80 & R<150);  %把模版转换成0、1阵列，待修复为1，其他为0；同时把数值转为double型，使之可以参与MATLAB里的数值计算。
%系数设定
p=1;
dt=1;%时间步长
I=R;
tic   %计时开始
for n=1:2000  %迭代次数，其值可以任选
    for i=3:nrow-2
        for j=3:ncol-2  %选择开始点和终点，以消除在图像边缘，进行中心差分法时，数值溢出；因为待修复区在图像内部，去除图像边缘点，不影响结果
            if mask(i,j)==1
            %求J_e^1
            deta_I1=((I(i,j+1)-I(i,j))^2+(0.25*(I(i-1,j+1)+I(i-1,j)-I(i+1,j+1)-I(i+1,j)))^2+eps)^(1/2); 
            k1_x=(I(i,j+2)+I(i,j)-2*I(i,j+1))/((I(i,j+2)+I(i,j)-2*I(i,j+1))^2+(I(i-1,j+1)+I(i+1,j+1)-2*I(i,j+1))^2+eps)^(1/2)...
                -((I(i,j+1)+I(i,j-1)-2*I(i,j))/((I(i,j+1)+I(i,j-1)-2*I(i,j))^2+(I(i-1,j)+I(i+1,j)-2*I(i,j))^2+eps)^(1/2));
            k1_y=0.25*((I(i-2,j+1)+I(i,j+1)-2*I(i-1,j+1))/((I(i-1,j+2)+I(i-1,j)-2*I(i-1,j+1))^2+(I(i-2,j+1)+I(i,j+1)-2*I(i-1,j+1))^2+eps)^(1/2)...
                 +(I(i-2,j)+I(i,j)-2*I(i-1,j))/((I(i-1,j+1)+I(i-1,j-1)-2*I(i-1,j))^2+(I(i-2,j)+I(i,j)-2*I(i-1,j))^2+eps)^(1/2)...
                 -(I(i,j+1)+I(i+2,j+1)-2*I(i+1,j+1))/((I(i+1,j+2)+I(i+1,j)-2*I(i+1,j+1))^2+(I(i,j+1)+I(i+2,j+1)-2*I(i+1,j+1))^2+eps)^(1/2)...
                 -(I(i,j)+I(i+2,j)-2*I(i+1,j))/((I(i+1,j+1)+I(i+1,j-1)-2*I(i+1,j))^2+(I(i,j)+I(i+2,j)-2*I(i+1,j))^2+eps)^(1/2));
             k1=abs(k1_x+k1_y);
             Ve_1=(I(i,j+1)-I(i,j))*(k1^p)/deta_I1;
             %求J_w^1
             deta_I2=((I(i,j)-I(i,j-1))^2+(0.25*(I(i-1,j-1)+I(i-1,j)-I(i+1,j-1)-I(i+1,j)))^2+eps)^(1/2);
             k2_x=(I(i,j+1)+I(i,j-1)-2*I(i,j))/((I(i,j+1)+I(i,j-1)-2*I(i,j))^2+(I(i-1,j)+I(i+1,j)-2*I(i,j))^2+eps)^(1/2)...
                 -(I(i,j)+I(i,j-2)-2*I(i,j-1))/((I(i,j)+I(i,j-2)-2*I(i,j-1))^2+(I(i+1,j-1)+I(i-1,j-1)-2*I(i,j-1))^2+eps)^(1/2);
             k2_y=0.25*((I(i-2,j-1)+I(i,j-1)-2*I(i-1,j-1))/((I(i-2,j-1)+I(i,j-1)-2*I(i-1,j-1))^2+(I(i-1,j)+I(i-1,j-2)-2*I(i-1,j-1))^2+eps)^(1/2)...
                   +(I(i-2,j)+I(i,j)-2*I(i-1,j))/((I(i-2,j)+I(i,j)-2*I(i-1,j))^2+(I(i-1,j-1)+I(i-1,j+1)-2*I(i-1,j))^2+eps)^(1/2)...
                   -(I(i,j-1)+I(i+2,j-1)-2*I(i+1,j-1))/((I(i,j-1)+I(i+2,j-1)-2*I(i+1,j-1))^2+(I(i+1,j)+I(i+1,j-2)-2*I(i+1,j-1))^2+eps)^(1/2)...
                   -(I(i,j)+I(i+2,j)-2*I(i+1,j))/((I(i,j)+I(i+2,j)-2*I(i+1,j))^2+(I(i+1,j+1)+I(i+1,j-1)-2*I(i+1,j))^2+eps)^(1/2));
             k2=abs(k2_x+k2_y);
             Vw_1=(I(i,j)-I(i,j-1))*(k2^p)/deta_I2;
             %求J_n^2
             deta_I3=((0.25*(I(i-1,j+1)+I(i,j+1)-I(i-1,j-1)-I(i,j-1)))^2+(I(i,j)-I(i-1,j))^2+eps)^(1/2);
             k3_x=0.25*((I(i-1,j+2)+I(i-1,j)-2*I(i-1,j+1))/((I(i-1,j+2)+I(i-1,j)-2*I(i-1,j+1))^2+(I(i-2,j+1)+I(i,j+1)-2*I(i-1,j+1))^2+eps)^(1/2)...
                  +(I(i,j+2)+I(i,j)-2*I(i,j+1))/((I(i,j+2)+I(i,j)-2*I(i,j+1))^2+(I(i-1,j+1)+I(i+1,j+1)-2*I(i,j+1))^2+eps)^(1/2)...
                  -(I(i-1,j)+I(i-1,j-2)-2*I(i-1,j-1))/((I(i-1,j)+I(i-1,j-2)-2*I(i-1,j-1))^2+(I(i,j-1)+I(i-2,j-1)-2*I(i-1,j-1))^2+eps)^(1/2)...
                  -(I(i,j)+I(i,j-2)-2*I(i,j-1))/((I(i,j)+I(i,j-2)-2*I(i,j-1))^2+(I(i+1,j-1)+I(i-1,j-1)-2*I(i,j-1))^2+eps)^(1/2));
             k3_y= (I(i-2,j)+I(i,j)-2*I(i-1,j))/((I(i-2,j)+I(i,j)-2*I(i-1,j))^2+(I(i-1,j+1)+I(i-1,j-1)-2*I(i-1,j))^2+eps)^(1/2)...
                     -(I(i-1,j)+I(i+1,j)-2*I(i,j))/((I(i-1,j)+I(i+1,j)-2*I(i,j))^2+(I(i,j+1)+I(i,j-1)-2*I(i,j))^2+eps)^(1/2);
             k3=abs(k3_x+k3_y);
             Vn_2=(I(i-1,j)-I(i,j))*(k3^p)/deta_I3;
             %求J_s^2;
             deta_I4=(((I(i+1,j+1)+I(i,j+1)-I(i+1,j-1)-I(i,j-1))*0.25)^2+(I(i+1,j)-I(i,j))^2+eps)^(1/2);
             k4_x=0.25*((I(i+1,j+2)+I(i+1,j)-2*I(i+1,j+1))/((I(i+1,j+2)+I(i+1,j)-2*I(i+1,j+1))^2+(I(i,j+1)+I(i+2,j+1)-2*I(i+1,j+1))^2+eps)^(1/2)...
                  +(I(i,j+2)+I(i,j)-2*I(i,j+1))/((I(i,j+2)+I(i,j)-2*I(i,j+1))^2+(I(i-1,j+1)+I(i+1,j+1)-2*I(i,j+2))^2+eps)^(1/2)...
                  -(I(i+1,j)+I(i+1,j-2)-2*I(i+1,j-1))/((I(i+1,j)+I(i+1,j-2)-2*I(i+1,j-1))^2+(I(i,j-1)+I(i+2,j-1)-2*I(i+1,j-1))^2+eps)^(1/2)...
                  -(I(i,j)+I(i,j-2)-2*I(i,j-1))/((I(i,j)+I(i,j-2)-2*I(i,j-1))^2+(I(i+1,j-1)+I(i-1,j-1)-2*I(i,j-1))^2+eps)^(1/2));
             k4_y=(I(i-1,j)+I(i+1,j)-2*I(i,j))/((I(i-1,j)+I(i+1,j)-2*I(i,j))^2+(I(i,j+1)+I(i,j-1)-2*I(i,j))^2+eps)^(1/2)...
                   -(I(i,j)+I(i+2,j)-2*I(i+1,j))/((I(i,j)+I(i+2,j)-2*I(i+1,j))^2+(I(i+1,j+1)+I(i+1,j-1)-2*I(i+1,j))^2+eps)^(1/2);
              k4=abs(k4_x+k4_y);
              Vs_2=(I(i,j)-I(i+1,j))*(k4^p)/deta_I4;
              %求delta_J
              deta_J=Ve_1-Vw_1+Vn_2-Vs_2;
              %迭代公式
              I(i,j)=I(i,j)+dt*deta_J; 
        end
    end
  end 
end
toc  %计时结束，输出循环所用时间
figure,imshow(uint8(R));%输出原待修复图
figure,imshow(mask);    %输出模版
figure,imshow(uint8(I)); %输出修复后图像