function fig10_15

%  image separation: includes figure 10.16

tic

% read in original images
x1=imread('flowers.jpg');
x1 = rgb2gray(x1);
I1 = double(x1);
size1=size(I1)

x2=imread('fox2.jpg');
x2 = rgb2gray(x2);
I2B = double(x2);
I2 = I2B(1:720,:);
size2=size(I2)

figure(1)
colormap(gray)
imshow(I1,[],'Border','tight')
%imagesc(I1)
axis image
axis off

figure(2)
colormap(gray)
imshow(I2,[],'Border','tight')
%imagesc(I2)
axis image
axis off

MM=size1(1); NN=size1(2);
n=MM*NN

% mixing matrix
M11=0.75; M12=0.25;
M21=0.9; M22=0.1;
M=[[M11 M12];[M21 M22]];

%%% calculate and show the mixed images
III1=M11*I1+M12*I2;
III2=M21*I1+M22*I2;

figure(3)
colormap(gray)
imshow(III1,[],'Border','tight')
%imagesc(III1)
axis image
axis off

figure(4)
colormap(gray)
imshow(III2,[],'Border','tight')
%imagesc(III2)
axis image
axis off

%%% center image data
xM1=sum(III1,'all')/n;
xM2=sum(III2,'all')/n;
X1=III1-xM1;
X2=III2-xM2;

%%% compute product matrix
X11=sum(X1.*X1,'all')/n;
X12=sum(X1.*X2,'all')/n;
X22=sum(X2.*X2,'all')/n;

XX=[ [X11 X12];[X12 X22]]

%%% find factorization of inner product matrix
[M1,D1]=eig(XX);

if D1(1,1) >= D1(2,2)
    d1=D1(1,1); d2=D1(2,2);
    Q1=M1(:,1)/norm(M1(:,1),2);
    Q2=M1(:,2)/norm(M1(:,2),2);
else
    d2=D1(1,1); d1=D1(2,2);
    Q2=M1(:,1)/norm(M1(:,1),2);
    Q1=M1(:,2)/norm(M1(:,2),2);
end

Q=[Q1 Q2];
D5=zeros(2,2); D5(1,1)=1/sqrt(d1); D5(2,2)=1/sqrt(d2);
QD=D5*Q';

%%% find y_i's
Y1=QD(1,1)*X1+QD(1,2)*X2;
Y2=QD(2,1)*X1+QD(2,2)*X2;

%%% compute kurtosis as function of angle
ntheta=360*4
%ntheta=360*16
theta=linspace(0,2*pi,ntheta);
for ia=1:ntheta
    V=[[ cos(theta(ia)) -sin(theta(ia))];[ sin(theta(ia)) cos(theta(ia))]];
    K1=V(1,1)*Y1+V(1,2)*Y2;
    K2=V(2,1)*Y1+V(2,2)*Y2;
    KKs1=sum(K1.^4,'all')/n;
    KKs2=sum(K2.^4,'all')/n;
    K(ia)=abs(KKs1-3)+abs(KKs2-3);
    %K(ia)=(KKs1-3)^2+(KKs2-3)^2;
end

% locate max points and then pick one for separation
iz=0;
for ia=2:ntheta-1
    if K(ia-1) < K(ia) && K(ia+1) < K(ia)
        iz=iz+1;
        thetaMM(iz)=theta(ia);
    end
end
thetaMM

%theta0=thetaMM(2)
theta0=thetaMM(4)

figure(5)
clf
%get(gcf);
set(gcf,'Position', [25 1115 658 230])

co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

plot(theta,K,'LineWidth',1.4)

%axis([0 2*pi 1.5 3])
xlabel('\theta-axis')
ylabel('K -axis')
set(gca,'xtick',[0 pi/2 pi 3*pi/2 2*pi])
set(gca,'XTickLabel',{'0';'\pi/2';'\pi';'3\pi/2';'2\pi'})
grid on
box on
set(gca,'FontSize',16,'FontWeight','bold')

%%% find mixing matrix and then separate images
V=[[ cos(theta0) -sin(theta0)];[ sin(theta0) cos(theta0)]];
DD=zeros(2,2); DD(1,1)=sqrt(d1); DD(2,2)=sqrt(d2);
A=Q*DD*V'

II1=V(1,1)*Y1+V(1,2)*Y2;
II2=V(2,1)*Y1+V(2,2)*Y2;

figure(6)
clf
colormap(gray)
imshow(II1,[],'Border','tight')
%imagesc(II1)
axis image
axis off

figure(7)
clf
colormap(gray)
imshow(II2,[],'Border','tight')
%imagesc(II2)
axis image
axis off

toc





















