%% shamir46recover.m - Shamir (4,6)-threshold recovering
%
%     The Elapsed time ~ 2 seconds for sharing and recovering
% (0) Read an RGB-based input image in Tiff format, e.g., Lenna.tiff
% 
m=512; n=512; np=m*n; m2=m/2;  n2=n/2;
I=imread('Peppers.tiff');        IC=I;        % getting pixel values > 250  
R=I(:,:,1); G=I(:,:,2); B=I(:,:,3); 
R0=zeros(m,n); G0=zeros(m,n); B0=zeros(m,n);  % for reconstruction
%
% (1) Prestore original pixel values great than 250 for a perfect recovery
%
for i=1:m
    for j=1:n
        if (R(i,j)>250)  R0(i,j)=R(i,j)-250; R(i,j)=250;  end
        if (G(i,j)>250)  G0(i,j)=G(i,j)-250; G(i,j)=250;  end
        if (B(i,j)>250)  B0(i,j)=B(i,j)-250; B(i,j)=250;  end
    end
end
%
% (2) Establish the table of ax=1 mod p=251; x is in [1,250]
%
p=251;
sinv=zeros(1,250);
for i=1:p-1
    s=i; k=1;
    while (k<p)
        if (mod(s*k,p)==1)
            sinv(i)=k; k=p;
        else
            k=k+1;
        end
    end
end
%
% (3) Recovering based on shadow images 1~4 and corresponding u=v(1:4)
%
S1=imread('shadow1.tif'); S2=imread('shadow2.tif');
S3=imread('shadow3.tif'); S4=imread('shadow4.tif');
R1=S1(:,:,1); G1=S1(:,:,2); B1=S1(:,:,3);         % S1: color shadow 1
R2=S2(:,:,1); G2=S2(:,:,2); B2=S2(:,:,3);         % S2: color shadow 2
R3=S3(:,:,1); G3=S3(:,:,2); B3=S3(:,:,3);         % S1: color shadow 1
R4=S4(:,:,1); G4=S4(:,:,2); B4=S4(:,:,3);         % S2: color shadow 2
%
% (3a) need to call adj.m to compute the adjoint of Vandermonde matrix A
%
u=[37, 60, 101, 149];         % v(1:6)=(37, 60, 101, 149, 203, 246)
u2=mod(u.^2,p); u3=mod(u.^3,p);
A=[1,1,1,1; u; u2; u3];
ti=mod(round(det(A)),p);
if (ti<0) t=ti+p; else t=ti; end;
s=sinv(t);
C=mod(s*adj(A),p); y=zeros(1,4);
for ix=1:m2
    i=2*(ix-1);
    for jy=1:n2
        j=2*(jy-1);
        y(1)=R1(ix,jy); y(2)=R2(ix,jy); y(3)=R3(ix,jy); y(4)=R4(ix,jy);
        r=mod(y*C,p);
        R(i+1,j+1)=r(1); R(i+1,j+2)=r(2); R(i+2,j+1)=r(3); R(i+2,j+2)=r(4);
        y(1)=G1(ix,jy); y(2)=G2(ix,jy); y(3)=G3(ix,jy); y(4)=G4(ix,jy);
        g=mod(y*C,p);
        G(i+1,j+1)=g(1); G(i+1,j+2)=g(2); G(i+2,j+1)=g(3); G(i+2,j+2)=g(4);
        y(1)=B1(ix,jy); y(2)=B2(ix,jy); y(3)=B3(ix,jy); y(4)=B4(ix,jy);
        b=mod(y*C,p);
        B(i+1,j+1)=b(1); B(i+1,j+2)=b(2); B(i+2,j+1)=b(3); B(i+2,j+2)=b(4);
    end
end
%
% (3a) IC holds a perfect revealed image
%
for i=1:m
    for j=1:n
            R0(i,j)=R0(i,j)+mod(R(i,j),p);
            G0(i,j)=G0(i,j)+mod(G(i,j),p); 
            B0(i,j)=B0(i,j)+mod(B(i,j),p);
    end
end
subplot(2,4,1)
imshow(S1)
xlabel('x_1=37')
subplot(2,4,2)
imshow(S2)
xlabel('x_2=60')
subplot(2,4,3)
imshow(S3)
xlabel('x_3=101')
subplot(2,4,4)
imshow(S4)
xlabel('x_4=149')
% Reconstructed image IC
IC(:,:,1)=R0; IC(:,:,2)=G0; IC(:,:,3)=B0; imwrite(IC,'tRecov.tif');
subplot(2,4,8)
imshow(IC)
xlabel('Recovered Image')
%
% (4) Histogram of Shadow Images 1,2,3
%
hr=zeros(250); hg=zeros(250); hb=zeros(250);
for i=1:m2
    for j=1:n2
        k=R1(i,j)+1; hr(k)=hr(k)+1;
        k=G1(i,j)+1; hg(k)=hg(k)+1;
        k=B1(i,j)+1; hb(k)=hb(k)+1;
    end
end
subplot(2,4,5)
L=0:249;
plot(L,hr,'r-',L,hg,'g-',L,hb,'b-')
xlabel('Histogram of Shadow 1')
ylabel('No. of Pixels')
hr=zeros(250); hg=zeros(250); hb=zeros(250);
for i=1:m2
    for j=1:n2
        k=R2(i,j)+1; hr(k)=hr(k)+1;
        k=G2(i,j)+1; hg(k)=hg(k)+1;
        k=B2(i,j)+1; hb(k)=hb(k)+1;
    end
end
subplot(2,4,6)
plot(L,hr,'r-',L,hg,'g-',L,hb,'b-')
xlabel('Histogram of Shadow 2')
ylabel('No. of Pixels')
hr=zeros(250); hg=zeros(250); hb=zeros(250);
for i=1:m2
    for j=1:n2
        k=R3(i,j)+1; hr(k)=hr(k)+1;
        k=G3(i,j)+1; hg(k)=hg(k)+1;
        k=B3(i,j)+1; hb(k)=hb(k)+1;
    end
end
subplot(2,4,7)
plot(L,hr,'r-',L,hg,'g-',L,hb,'b-')
xlabel('Histogram of Shadow 3')
ylabel('No. of Pixels')