%% crtrecover.m - Recover CRT-Based (k,n)=(3,5)-threshold sharing
%
%  Elapsed time 16.291997 Seconds 
%  in my PC running XP            (tic; crtrecover; toc)
% (1) Read an RGB-based input image in Tiff format, e.g., Lenna.tiff
% 
m=512; n=512; np=m*n;
I=imread('Lenna.tiff');       IC=I;          % (m * n) * 3 = (R, G, B)*3  
R=I(:,:,1); G=I(:,:,2); B=I(:,:,3); 
R0=zeros(m,n,'uint8'); G0=zeros(m,n,'uint8'); B0=zeros(m,n,'uint8'); 
% 
%  for a perfect reconstruction 
% (2) Read in 3 Shadow Images
%
S3=imread('Shadow3.tif');        % (m2*n2)*3 size for shadows 3,4,5
S4=imread('Shadow4.tif');
S5=imread('Shadow5.tif');
R3=S3(:,:,1);  G3=S3(:,:,2);  B3=S3(:,:,3);  % S3: Color Shadow 3
R4=S4(:,:,1);  G4=S4(:,:,2);  B4=S4(:,:,3);  % S1: Color Shadow 4
R5=S5(:,:,1);  G5=S5(:,:,2);  B5=S5(:,:,3);  % S2: Color Sahdow 5
%
R0=randi([0,1],m,n); G0=randi([0,1],m,n); B0=randi([0,1],m,n);
%
%  (2a) Shift 1 bit to the right for (R,G,B) and 
%       preserve the least significant bits
%  Guarantee a perfect reconstruction, otherwise very minor errors
for i=1:m
    for j=1:n
        r0=bitget(R(i,j),1);
        g0=bitget(G(i,j),1);
        b0=bitget(B(i,j),1);
        R(i,j)=round((R(i,j)-r0)/2);  R0(i,j)=r0;
        G(i,j)=round((G(i,j)-g0)/2);  G0(i,j)=g0;
        B(i,j)=round((B(i,j)-b0)/2);  B0(i,j)=b0;
    end
end
%
% (3) Print 3 color shadow images
%
subplot(2,4,1)
imshow(S3)
title('Shadow 3')
subplot(2,4,2)
imshow(S4)
title('Shadow 4')
subplot(2,4,3)
imshow(S5)
title('Shadow 5')
hr=zeros(256); hg=zeros(256); hb=zeros(256);
for i=1:m
    for j=1:n
        r=1+R3(i,j);
        hr(r)=hr(r)+1;
        g=1+G3(i,j);
        hg(g)=hg(g)+1;
        b=1+B3(i,j);
        hb(b)=hb(b)+1;
    end
end 
subplot(2,4,5)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Shadow Image 3')
hr=zeros(256); hg=zeros(256); hb=zeros(256);
for i=1:m
    for j=1:n
        r=1+R4(i,j);
        hr(r)=hr(r)+1;
        g=1+G4(i,j);
        hg(g)=hg(g)+1;
        b=1+B4(i,j);
        hb(b)=hb(b)+1;
    end
end 
subplot(2,4,6)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Shadow Image 4')
hr=zeros(256); hg=zeros(256); hb=zeros(256);
for i=1:m
    for j=1:n
        r=1+R5(i,j);
        hr(r)=hr(r)+1;
        g=1+G5(i,j);
        hg(g)=hg(g)+1;
        b=1+B5(i,j);
        hb(b)=hb(b)+1;
    end
end 
subplot(2,4,7)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Shadow Image 5')
%
% (4) Recover secret image from Shadow Images 3,4,5 based on CRT
%
m0=128; S=[247,251, 253, 254,255]; %(m0,m1,m2,m3,m4,m5)in IS3C/2016 paper
m0=uint32(128); S=uint32(S);
M=uint32(S(3)*S(4)*S(5));
Z=uint32([S(4)*S(5), S(3)*S(5), S(3)*S(4)]); 
Y=uint32([127, 253, 128]);
A=uint32([255,255,255]);
u=uint32(0);  s=uint32(0);
for i=1:m
    for j=1:n
        A(1)=R3(i,j); A(2)=R4(i,j); A(3)=R5(i,j);
        s=0;
        % for k=1:3
            %u=mod(A(k)*mod(Y(k)*Z(k),M),M);
            %while(u>M) u=u-M; end
            %s=mod(s+u,M);
            %while(s>M) s=s-M; end
        % end
        s=sum(A.*(Y.*Z));
        s=mod(s,M);
        R(i,j)=mod(s,m0);
        
        A(1)=G3(i,j); A(2)=G4(i,j); A(3)=G5(i,j);
        s=0;
        for k=1:3
            u=mod(A(k)*mod(Y(k)*Z(k),M),M);
            while(u>M) u=u-M; end
            s=s+u;
            while(s>M) s=s-M; end
        end
        s=mod(s,M);
        G(i,j)=mod(s,m0);
        
        A(1)=B3(i,j); A(2)=B4(i,j); A(3)=B5(i,j);
        s=0;
        for k=1:3
            u=mod(A(k)*mod(Y(k)*Z(k),M),M);
            while(u>M) u=u-M; end
            s=s+u;
            while(s>M) s=s-M; end
        end
        s=mod(s,M);
        B(i,j)=mod(s,m0);
    end
end
%
% (5) IC holds a perfect revealed image (with the intensity values >250)
%
for i=1:m
    for j=1:n
            R0(i,j)=R0(i,j)+2*R(i,j);
            G0(i,j)=G0(i,j)+2*G(i,j); 
            B0(i,j)=B0(i,j)+2*B(i,j);
    end
end
IC(:,:,1)=R0; IC(:,:,2)=G0; IC(:,:,3)=B0; imwrite(IC,'Recov345.tif');
subplot(2,4,4)
imshow(IC)
title('Recovering Image by Shadows 3\sim5')
hr=zeros(256); hg=zeros(256); hb=zeros(256);
for i=1:m
    for j=1:n
        r=1+R0(i,j);
        hr(r)=hr(r)+1;
        g=1+G0(i,j);
        hg(g)=hg(g)+1;
        b=1+B0(i,j);
        hb(b)=hb(b)+1;
    end
end 
subplot(2,4,8)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Recovered Secret Image')