lab5.m

采用总变分法（TV）进行图像修复的matlab源代码。包括热扩散修复和多尺度方法。

clear all;
close all;


%% CODE FOR PROBLEM 1

%% FOR IMAGE A

A = double(imread('A_Lab5.jpg'));
figure;
imshow(uint8(A));

thresh = 0.3;
cannyA = edge(A,'canny',thresh);
prewittA = edge(A,'prewitt');
sobelA = edge(A,'sobel');
laplacianA = edge(A,'log');
zeroA = edge(A,'zerocross');

figure;
subplot(2,2,1)
imshow(uint8(cannyA));
subplot(2,2,2)
imshow(uint8(prewittA));
subplot(2,2,3)
imshow(uint8(sobelA));
subplot(2,2,4)
imshow(uint8(laplacianA));

%% FOR IMAGE B

B = double(imread('B_Lab5.jpg'));
figure;
imshow(uint8(B));

thresh = 0.3;
cannyB = edge(B,'canny',thresh);
prewittB = edge(B,'prewitt');
sobelB = edge(B,'sobel');
laplacianB = edge(B,'log');
zeroB = edge(B,'zerocross');

figure;
subplot(2,2,1)
imshow(uint8(cannyB));
subplot(2,2,2)
imshow(uint8(prewittB));
subplot(2,2,3)
imshow(uint8(sobelB));
subplot(2,2,4)
imshow(uint8(laplacianB));



%% CODE FOR PROBLEM 2
width=10;
sigma = 1/2;
[X,Y] = meshgrid(-(width)/2:(width)/2,-(width)/2:(width)/2);
K=exp(X.^2/-2/sigma).*exp(Y.^2/-2/sigma)/(pi*2*sigma);
smoothA=conv2(A,K);


figure;
imshow(uint8(smoothA));

[Ax,Ay] = gradient(smoothA);
MgradA = sqrt(Ax.^2 + Ay.^2);
figure;
imshow(MgradA);
edgesA=zeros(size(smoothA));

M1=0
M2=0
for i=1:size(smoothA,1)-1
    for j=1:size(smoothA,2)-1
        if  M1 < Ax(i+1,j+1)^2 + Ay(i+1,j+1)^2;
            M2=M1;
            M1= Ax(i+1,j+1)^2 + Ay(i+1,j+1)^2;
        else
            M1=M1;
        end
    end
end
maxM=(M1+M2)/2;

aveM = sqrt(dot(dot(Ax,Ax),dot(Ax,Ax)) + dot(dot(Ay,Ay),dot(Ay,Ay)))/(size(smoothA,1)*size(smoothA,2));

threshA = maxM/aveM;


for i=1:size(smoothA,1)
    for j=1:size(smoothA,2)
        if (Ax(i,j)^2 + Ay(i,j)^2) < threshA
            edgesA(i,j) =255;
        end
    end
end

figure;
imshow(uint8(edgesA));


%% FOR IMAGE B

width=10;
sigma = 1/2;
[X,Y] = meshgrid(-(width)/2:(width)/2,-(width)/2:(width)/2);
K=exp(X.^2/-2/sigma).*exp(Y.^2/-2/sigma)/(pi*2*sigma);
smoothB=conv2(B,K);


figure;
imshow(uint8(smoothB));

[Bx,By] = gradient(smoothB);
MgradB = sqrt(Bx .^2 + By .^2);
figure;
imshow(MgradB);
edgesB=zeros(size(smoothB));

M1=0;
M2=0;
for i=1:size(smoothB,1)-1
    for j=1:size(smoothB,2)-1
        if  M1 < Bx(i+1,j+1)^2 + By(i+1,j+1)^2;
            M2=M1;
            M1= Bx(i+1,j+1)^2 + By(i+1,j+1)^2;
        else
            M1=M1;
        end
    end
end
maxM=1/2*(M1+M2);

aveM = sqrt(dot(dot(Bx,Bx),dot(Bx,Bx)) + dot(dot(By,By),dot(By,By)))/(size(smoothB,1)*size(smoothB,2));

threshB = maxM/aveM;

edgesB=zeros(size(smoothB));
for i=1:size(smoothB,1)
    for j=1:size(smoothB,2)
        if (Bx(i,j)^2 + By(i,j)^2) < threshB
            edgesB(i,j) =255;
        end
    end
end

figure;
imshow(uint8(edgesB));
  

%% Thinning Edges


%% Thinning Image A Edges

              
thetaA = mod(round(mod((atan2(Ay,Ax))/pi*8,8)),8);

ka=zeros(size(thetaA));
la=zeros(size(thetaA));
ka(thetaA==0)=1;
ka(thetaA==2)=-1;
ka(thetaA==3)=-1;
ka(thetaA==4)=-1;
ka(thetaA==6)=1;
ka(thetaA==7)=1;


la(thetaA==0)=-1;
la(thetaA==1)=-1;
la(thetaA==2)=-1;
la(thetaA==4)=1;
la(thetaA==5)=1;
la(thetaA==6)=1;

edgesA2=edgesA;
for i=2:size(smoothA,1)-1
    for j=2:size(smoothA,2)-1
        if edgesA(i,j) == 0;
            if (MgradA(i,j) < MgradA(i-la(i,j),j-ka(i,j))) | (MgradA(i,j) < MgradA(i+la(i,j),j+ka(i,j)));

edgesA2(i,j)=255;
        end
    end
end
end
figure;
imshow(uint8(edgesA2));


%% Thinning Image B Edges

             
thetaB = mod(round(mod((atan2(By,Bx))/pi*8,8)),8);

kb=zeros(size(thetaB));
lb=zeros(size(thetaB));
kb(thetaB==0)=1;
kb(thetaB==2)=-1;
kb(thetaB==3)=-1;
kb(thetaB==4)=-1;
kb(thetaB==6)=1;
kb(thetaB==7)=1;


lb(thetaB==0)=-1;
lb(thetaB==1)=-1;
lb(thetaB==2)=-1;
lb(thetaB==4)=1;
lb(thetaB==5)=1;
lb(thetaB==6)=1;


edgesB2=edgesB;
for i=2:size(smoothB,1)-1
    for j=2:size(smoothB,2)-1
        if edgesB(i,j) == 0;
            if (MgradB(i,j) < MgradB(i-lb(i,j),j-kb(i,j))) | (MgradB(i,j) < MgradB(i+lb(i,j),j+kb(i,j)))
                edgesB2(i,j)=255;
        end
    end
end
end
figure;
imshow(uint8(edgesB2));



%% CODE FOR PROBLEM 5

load k_coef
k_scale=[175 420 420 70 100 420 60 140 140 60];
width=10;
sigma=3;
PSF=fspecial('gaussian',width,sigma);
smoothA=imfilter(A,PSF,'conv');


figure;
for i=1:10;
    ki=eval(['k', int2str(i)]);
    k(:,:,i)=filter2(ki/k_scale(i),smoothA);
    subplot(3,4,i);
    imagesc(k(:,:,i));
end

sin_theta=k(:,:,2)./(sqrt(k(:,:,2).^2+k(:,:,3).^2));
cos_theta=k(:,:,3)./(sqrt(k(:,:,2).^2+k(:,:,3).^2));
figure
imagesc(sin_theta)
figure
imagesc(cos_theta)

a=2;
h=2*a-1;
for i=1:size(smoothA,1);
    for j=1:size(smoothA,2);
        for x=1:h;
            for y=1:h;
                fprime(i,j,x,y) = (3*k(i,j,2)*k(i,j,3)+2*k(i,j,3)*k(i,j,4)*(x-h/2)/h +k(i,j,3)*k(i,j,5)*(y-h/2)/h +3*k(i,j,3)*k(i,j,7)*(x-h/2)^2/h^2+2*k(i,j,3)*k(i,j,8)*(x-h/2)*(y-h/2)/h^2+k(i,j,3)*k(i,j,9)*(y-h/2)^2/h^2+k(i,j,2)*k(i,j,5)*(x-h/2)/h+2*k(i,j,2)*k(i,j,6)*(y-h/2)/h+k(i,j,2)*k(i,j,8)*(x-h/2)/h^2+2*k(i,j,2)*k(i,j,9)*(x-h/2)*(y-h/2)/h^2+3*k(i,j,2)*k(i,j,10)*(y-h/2)^2/h^2)/(sqrt(k(i,j,2)^2+k(i,j,3)^2));
            end
        end
    end
end


for i=1:size(smoothA,1);
    for j=1:size(smoothA,2);
        for x=1:h;
            for y=1:h;
                fdblprime(i,j,x,y) = 2*(3*k(i,j,7)*sin_theta(i,j)^2 + 2*k(i,j,8)*sin_theta(i,j)*cos_theta(i,j) + k(i,j,9)*cos_theta(i,j)^2)*(x-h/2)/h + 2*(k(i,j,8)*sin_theta(i,j)^2 + 2*k(i,j,9)*sin_theta(i,j)*cos_theta(i,j) + 3*k(i,j,10)*cos_theta(i,j)^2)*(y-h/2)/h + 2*(k(i,j,4)*sin_theta(i,j)^2 + k(i,j,5)*sin_theta(i,j)*cos_theta(i,j) + k(i,j,6)*cos_theta(i,j)^2);
            end
        end
    end
end

%f(i,j) = k(i,j,1)+k(i,j,2)*i + k(i,j,3)*j + k(i,j,4)*i^2 + k(i,j,5)*i*j + k(i,j,6)*j^2 + k(i,j,7)*+k(i,j,8) k(i,j,9) k(i,j,10);

cubic_edges = zeros(size(smoothA));

for i=1:size(smoothA,1)-1;
    for j=1:size(smoothA,2)-1;
        for m=-a/2:a/2;
            for n=-a/2:a/2;
                if (fdblprime(i,j,a+m,a+n)/fdblprime(i,j,a,a))<0 & abs(fprime(i,j,x,y))>2;
                    cubic_edges(i,j)=255;
                end
            end
        end
    end
end
figure;
imshow(cubic_edges);



% 
% 
% 
% [dx,dy] = gradient(A);
% f1=dx.*cos_theta+dy.*sin_theta;
% cf1=f1>0;
% cf2=zeros(size(cf1));
% Af=6*(3*k(:,:,7).*sin_theta.^2+k(:,:,8).*sin_theta.^2.*cos_theta+k(:,:,9).*sin_theta.*cos_theta.^2+k(:,:,1).*cos_theta.^3);
% Bf=2*(k(:,:,4).*sin_theta.^2+k(:,:,5).*sin_theta.*cos_theta+k(:,:,6).*cos_theta.^2);
% for r=0:.01:0.7
% f2=Af*r+Bf;
% subplot(3,1,1)
% imagesc(f2)
% drawnow
% subplot(3,1,2)
% cf2=cf2|(abs(f2)<0.5).*cf1;
% imagesc(abs(f2)<0.5)
% drawnow
% 
% subplot(3,1,3)
% imagesc(cf2)
% colormap gray
% drawnow
% 
% end
% figure
% imagesc((cf2>0).*cf1)
% colormap gray
% 
% figure
% subplot(2,1,1)
% imagesc(Af)
% 
% subplot(2,1,2)
% imagesc(Bf)
% 
% 
% f2dblprime = simplify(2*(3*k(i,j,7)*sin_theta^2 + 2*k(i,j,8)*sin_theta*cos_theta + k(i,j,9)*cos_theta^2)*x + 2*(k(i,j,8)*sin_theta^2 + 2*k(i,j,9)*sin_theta*cos_theta + 3*k(i,j,10)*cos_theta^2)*y + 2*(k(i,j,4)*sin_theta^2 + k(i,j,5)*sin_theta*cos_theta + k(i,j,6)*cos_theta^2));
% x = -(y*k(i,j,8)*k(i,j,2)^2+2*y*k(i,j,9)*k(i,j,2)*k3+3*y*k(i,j,10)*k3^2+k(i,j,4)*k(i,j,2)^2+k(i,j,5)*k(i,j,2)*k3+k(i,j,6)*k3^2)/(3*k(i,j,7)*k(i,j,2)^2+2*k(i,j,8)*k(i,j,2)*k3+k(i,j,9)*k3^2);
% fprime = (3*k(i,j,2)*k3+2*k3*k4*x+k3*k5*y+3*k3*k(i,j,7)*x^2+2*k3*k(i,j,8)*x*y+k3*k(i,j,9)*y^2+k(i,j,2)*k5*x+2*k(i,j,2)*k(i,j,6)*y+k(i,j,2)*k(i,j,8)*x^2+2*k(i,j,2)*k(i,j,9)*x*y+3*k(i,j,2)*k(i,j,10)*y^2)/(sqrt(k(i,j,2)^2+k3^2));
% 
%