%DE-BLURING color

%Este progrma realiza el debluring de una imagen a color de 256x256,
%obteniendo el parametro de regularizacion automaticamente mediante GCV

% Nota: para que este programa funcione es necesario tener las funciones:
% blur
% GCVahi

%Andr閟 Hoyos Idrobo     marzo de 2009

clear all;
close all;
clc;

%Carga de imagen 

x=double(imread('House1.png'));
x=x(1:256,1:256,:);
[nn,mm,kk]=size(x);

%R es el porcentaje de ruido
R=0.0001;
%Vector alfa
alfa=0.000001:0.001:0.03;

t=20;
ss=2.5;

tipo='exp';
[F]=blur(nn,mm,kk,t,ss,tipo);


G=zeros(nn,mm,kk);

for iu=1:3
G(:,:,iu)=F*x(:,:,iu);
end

G=cat(3,G(:,:,1),G(:,:,2),G(:,:,3));

ruido=R*255*randn(size(G));
bruido1=(G+ruido);

b1=bruido1(:,:,1);
b2=bruido1(:,:,2);
b3=bruido1(:,:,3);

%Se llama la funcion GCVahi
%Para R
[V_alfa1,U1,s1,V1,iT,E11]=GCVahi(F,b1,alfa);
%Para G
[V_alfa2,U2,s2,V2,iT,E12]=GCVahi(F,b2,alfa);
%Para B
[V_alfa3,U3,s3,V3,iT,E13]=GCVahi(F,b3,alfa);
 

KK1=V_alfa1(1,2:iT);
dmin1=(min(KK1));
df1=find(V_alfa1==dmin1);
alfa11=abs(alfa(df1));

KK2=V_alfa2(1,1:iT);
dmin2=(min(KK2));
df2=find(V_alfa2==dmin2);
alfa12=abs(alfa(df2));

KK3=V_alfa3(1,2:iT);
dmin3=(min(KK3));
df3=find(V_alfa3==dmin3);
alfa13=abs(alfa(df3));

H1=s1'*s1;
pr1=(alfa11)*E11; 
kj1=s1*inv(H1+pr1);
xref1=uint8((V1*kj1*U1')*b1);

H2=s2'*s2;
pr2=(alfa12)*E12; 
kj2=s2*inv(H2+pr2);
xref2=uint8((V2*kj2*U2')*b2);

H3=s3'*s3;
pr3=(alfa13)*E13; 
kj3=s3*inv(H3+pr3);
xref3=uint8((V3*kj3*U3')*b3);

Xref=cat(3,xref1,xref2,xref3);
d1=num2str(alfa11,'%10.3e\n');
d2=num2str(alfa12,'%10.3e\n');
d3=num2str(alfa13,'%10.3e\n');

L=num2str(t,'%10.2f');
rho=num2str(ss,'%10.2f');
rui=num2str(R*100);

subplot(1,3,1)
imshow(uint8(x)),title('Imagen Original');
subplot(1,3,2)
imshow(uint8(bruido1)),title('Imagen adquirida');
xlabel(['L=',L,' Rho=',rho,' Ruido=',rui,'%']);
subplot(1,3,3)
imshow(uint8(Xref)),title('Regularizacion optima por Tikhonov');
xlabel(['alfa R=',d1 ;'alfa G=',d2;'alfa B=',d3]);

figure
plot(alfa,V_alfa1,'ob',alfa,V_alfa2,'*r',alfa,V_alfa3,'ks'),grid;
h=legend('Valfa R','Valfa G','Valfa B');

title('GCV');
xlabel('alfa');
ylabel('GCV(alfa)');