%% crtshare.m - CRT-Based (k,n)=(3,5)-threshold sharing
%
%  Elapsed time is 11.591415 sectonds
%  in my PC running XP (based on tic; crtshare; 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 perfect reconstruction
% (2) Image Sharing - generate 5 Shadows for (3,5)-thresholding
%
Shadow=imread('Peppers.tiff');        % (m2*n2)*3 size for shadows 1~5
R1=zeros(m,n,'uint32'); R2=R1; R3=R1; R4=R1; R5=R1;
G1=zeros(m,n,'uint32'); G2=G1; G3=G1; G4=G1; G5=G1;
B1=zeros(m,n,'uint32'); B2=B1; B3=B1; B4=B1; B5=B1;
%
%  (2a) Shift 1 bit to the right for (R,G,B) and 
%       preserve the least significant bits
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) Compute the shadow images
%
m0=128; S=[247,251, 253, 254,255]; %(m0,m1,m2,m3,m4,m5)in IS3C/2016 paper
m0=uint32(128); S=uint32(S);
Low=m0*S(4)*S(5); Up=S(1)*S(2)*S(3);
Lb=round(Low/m0+2);  Ub=round(Up/m0-3);
Len=Ub-Lb;  Rt=randi([Lb, Ub],m,n);
Rt=uint32(Rt); y=uint32(0);
for i=1:m
    for j=1:n
        r=R(i,j);  g=G(i,j);  b=B(i,j); t=Rt(i,j);
        y=t*m0+uint32(r);
        R1(i,j)=mod(y,S(1)); R2(i,j)=mod(y,S(2)); R3(i,j)=mod(y,S(3));
        R4(i,j)=mod(y,S(4)); R5(i,j)=mod(y,S(5));
        y=t*m0+uint32(g);
        G1(i,j)=mod(y,S(1)); G2(i,j)=mod(y,S(2)); G3(i,j)=mod(y,S(3));
        G4(i,j)=mod(y,S(4)); G5(i,j)=mod(y,S(5));
        y=t*m0+uint32(b);
        B1(i,j)=mod(y,S(1)); B2(i,j)=mod(y,S(2)); B3(i,j)=mod(y,S(3));
        B4(i,j)=mod(y,S(4)); B5(i,j)=mod(y,S(5));
    end
end
% R1=mod(T,S(1));
% R2=mod(T,S(2));
% R3=mod(T,S(3));
% R4=mod(T,S(4));
% R5=mod(T,S(5));
%
% (4) Print the color input image and its Shadow Color Images
%
S1=Shadow; S2=Shadow; S3=Shadow; S4=Shadow; S5=Shadow;
S1(:,:,1)=R1;  S1(:,:,2)=G1;  S1(:,:,3)=B1;  % S1: Color Shadow 1
S2(:,:,1)=R2;  S2(:,:,2)=G2;  S2(:,:,3)=B2;  % S2: Color Sahdow 2
S3(:,:,1)=R3;  S3(:,:,2)=G3;  S3(:,:,3)=B3;  % S3: Color Shadow 3
S4(:,:,1)=R4;  S4(:,:,2)=G4;  S4(:,:,3)=B4;  % S4: Color Shadow 4
S5(:,:,1)=R5;  S5(:,:,2)=G5;  S5(:,:,3)=B5;  % S5: Color Shadow 5
subplot(2,4,1)
imshow(I)
title('Secret Image')
subplot(2,4,2)
imshow(S1)
title('Shadow 1')
subplot(2,4,3)
imshow(S2)
title('Shadow 2')
subplot(2,4,4)
imshow(S3)
title('Shadow 3')
subplot(2,4,5)
imshow(S4)
title('Shadow 4')
subplot(2,4,6)
imshow(S5)
title('Shadow 5')
imwrite(S1,'shadow1.tif'); imwrite(S2,'shadow2.tif');
imwrite(S3,'shadow3.tif'); imwrite(S4,'shadow4.tif');
imwrite(S5,'shadow5.tif');
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,7)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Shadow Image 2')
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,8)
X=0:255;
plot(X,hr,'r-',X,hg,'g-',X,hb,'b-')
title('For 512 \times 512 Shadow Image 5')
%
% (5) Recovering - from S3~S6 Shadow images to solve the CRT problem
%
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=mod(s+u,M);
            %while(s>M) s=s-M; end
        % end
        s=sum(A.*(Y.*Z));
        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=mod(s+u,M);
            %while(s>M) s=s-M; end
        % end
        s=sum(A.*(Y.*Z));
        s=mod(s,M);
        B(i,j)=mod(s,m0);
    end
end