%% shamir46share.m - Shamir (4,6)-threshold sharing
%
%     The Elapsed time ~ 2 seconds for sharing and recovering
% (1) 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;  
R=I(:,:,1); G=I(:,:,2); B=I(:,:,3); 
R0=zeros(m,n); G0=zeros(m,n); B0=zeros(m,n);  % for reconstruction
%
% (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) Image Sharing - generate 6 Shadows for (4,6)-threshold mod p=251
%
Shadow=imread('color256.tif');     % format (m2*n2)*3 size for shadows 1~6
R1=zeros(m2,n2); R2=R1; R3=R1; R4=R1; R5=R1; R6=zeros(m2,n2);
G1=zeros(m2,n2); G2=G1; G3=G1; G4=G1; G5=G1; G6=zeros(m2,n2);
B1=zeros(m2,n2); B2=B1; B3=B1; B4=B1; B5=B1; B6=zeros(m2,n2);
% 
% (3a) Prepare a perfect recovering
%
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
%
% (4) Compute the shadow images
%
v=round(250*rand(1,6));      % parameters could be randomly generated
v=sort(v);
v=[37,60,101,149,203,246];   % specified {m_i} as an example only
v2=mod(v.^2,p); v3=mod(v.^3,p);
AX=[1,1,1,1,1,1; v; v2; v3];
r=zeros(1,4); g=zeros(1,4); b=zeros(1,4);
for i=0:2:m-2
    ix=1+i/2;
    for j=0:2:n-2
        jy=1+j/2;
        r(1)=R(i+1,j+1); r(2)=R(i+1,j+2); r(3)=R(i+2,j+1); r(4)=R(i+2,j+2);
        yr=mod(r*AX,p);
        R1(ix,jy)=yr(1); R2(ix,jy)=yr(2); R3(ix,jy)=yr(3);
        R4(ix,jy)=yr(4); R5(ix,jy)=yr(5); R6(ix,jy)=yr(6);
        g(1)=G(i+1,j+1); g(2)=G(i+1,j+2); g(3)=G(i+2,j+1); g(4)=G(i+2,j+2);
        yg=mod(g*AX,p);
        G1(ix,jy)=yg(1); G2(ix,jy)=yg(2); G3(ix,jy)=yg(3);
        G4(ix,jy)=yg(4); G5(ix,jy)=yg(5); G6(ix,jy)=yg(6);
        b(1)=B(i+1,j+1); b(2)=B(i+1,j+2); b(3)=B(i+2,j+1); b(4)=B(i+2,j+2);
        yb=mod(b*AX,p);
        B1(ix,jy)=yb(1); B2(ix,jy)=yb(2); B3(ix,jy)=yb(3);
        B4(ix,jy)=yb(4); B5(ix,jy)=yb(5); B6(ix,jy)=yb(6);
    end
end
%
% (5) Print the color input image and its Shadow Color Images
%
S1=Shadow; S2=Shadow; S3=Shadow; S4=Shadow; S5=Shadow; S6=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
S6(:,:,1)=R6;  S6(:,:,2)=G6;  S6(:,:,3)=B6;  % S5: Color Shadow 6
subplot(3,3,1)
imshow(I)
xlabel('Secret')
subplot(3,3,2)
imshow(S1)
xlabel('x_1=37')
subplot(3,3,3)
imshow(S2)
xlabel('x_2=60')
subplot(3,3,4)
imshow(S3)
xlabel('x_3=101')
subplot(3,3,5)
imshow(S4)
xlabel('x_4=149')
subplot(3,3,6)
imshow(S5)
xlabel('x_5=203')
subplot(3,3,7)
imshow(S6)
xlabel('x_6=246')
imwrite(S1,'shadow1.tif'); imwrite(S2,'shadow2.tif');
imwrite(S3,'shadow3.tif'); imwrite(S4,'shadow4.tif');
imwrite(S5,'shadow5.tif'); imwrite(S6,'shadow6.tif');
%
% (6) Histogram of Shadow Images 5 and  6
%
hr=zeros(250); hg=zeros(250); hb=zeros(250);
for i=1:m2
    for j=1:n2
        k=R5(i,j)+1; hr(k)=hr(k)+1;
        k=G5(i,j)+1; hg(k)=hg(k)+1;
        k=B5(i,j)+1; hb(k)=hb(k)+1;
    end
end
subplot(3,3,8)
L=0:249;
plot(L,hr,'r-',L,hg,'g-',L,hb,'b-')
xlabel('Histogram of Shadow 5')
ylabel('No. of Pixels')
hr=zeros(250); hg=zeros(250); hb=zeros(250);
for i=1:m2
    for j=1:n2
        k=R6(i,j)+1; hr(k)=hr(k)+1;
        k=G6(i,j)+1; hg(k)=hg(k)+1;
        k=B6(i,j)+1; hb(k)=hb(k)+1;
    end
end
subplot(3,3,9)
plot(L,hr,'r-',L,hg,'g-',L,hb,'b-')
xlabel('Histogram of Shadow 6')
ylabel('No. of Pixels')