% ucitavanje i priprema slike
ims1 = dyaddown(im2double(rgb2gray(imread('C:\MATLAB7\toolbox\images\imdemos\football.jpg'))),'m');

[r1,c1] = size(ims1); 

% slika suma
ims2 = randn(size(ims1)) / 7;

% priprema jednodimezionalnih vektora
datalen = r1*c1;
s1 = reshape(ims1, 1, datalen); s2 = reshape(ims2, 1, datalen); S = [s1; s2];

% prikaz ulaznih slika
figure(1), imshow(ims1, [min(s1), max(s1)]), figure(2), imshow(ims2, [min(s2), max(s2)]);

% priprema jednodimezionalnih vektora
datalen = r1*c1;
s1 = reshape(ims1, 1, datalen); s2 = reshape(ims2, 1, datalen); S = [s1; s2];

% mjereni signal s aditivnim sumom
imx1 = ims1 + ims2;
x1 = s1 + s2; 

% umjetni signal s wiener2
imx2 = wiener2(imx1);
x2 = reshape(imx2, 1, datalen);

X = [x1; x2];

% Prikaz izmijesanih slika
figure(3), imshow(imx1, [min(x1), max(x1)]), figure(4), imshow(imx2, [min(x2), max(x2)]);


% PCA
Xcov = cov(X');             % kovarijancijska matrica 2 x 2
[Xv, Xd] = eig(Xcov);       % eigen dekompozicija
Vx = -diag(diag(Xd).^(-1/2)) * Xv';    % dekorelacija i skaliranje signala (var = 1)

Z = Vx * X;
z1 = Z(1,:); z2 = Z(2,:);

figure(5), scatter(Z(1,:), Z(2,:))

% formiranje i prikaz razdvojenih slika PCA metodom
imz1 = reshape(z1, r1, c1); imz2 = reshape(z2, r1, c1);
figure(6), imshow(imz1, [min(z1), max(z1)]), figure(7), imshow(imz2, [min(z2), max(z2)]);

% ICA, rotacijska metoda by Debeli

% Ovisnost kurtosis o theta 
z1kurt = zeros(360,1); z2kurt = zeros(360,1); 
theta = [-pi: 2*pi/360: pi-2*pi/360];
for i = [1:360]
    R = [cos(theta(i)) -sin(theta(i)); sin(theta(i)) cos(theta(i))];
    Zr = R * Z;
    z1kurt(i) = kurtosis(Zr(1,:))-3;
    z2kurt(i) = kurtosis(Zr(2,:))-3;
end

sumAbsKurt = abs(z1kurt)+abs(z2kurt);
figure(8), plot(theta, [z1kurt z2kurt sumAbsKurt])

% nalazimo theta za maksimalnu sumu apsolutnih iznosa kurtosisa
[maksimum, ind] = max(sumAbsKurt);
theta = theta(ind(1));

% to je dobra matrica rotacije!
R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
Zr = R * Z;
figure(9), scatter(Zr(1,:), Zr(2,:))
zr1 = -Zr(1,:); zr2 = -Zr(2,:);
imzr1 = reshape(zr1, r1, c1); imzr2 = reshape(zr2, r1, c1);
figure(10), imshow(imzr1, [min(zr1), max(zr1)]), figure(11), imshow(imzr2, [min(zr2), max(zr2)]);
figure(12), imshow(-imzr1, [min(-zr1), max(-zr1)]), figure(13), imshow(-imzr2, [min(-zr2), max(-zr2)]);

m22 = sum(ims1.*ims1.*(-imzr1).*(-imzr1));
m2i = sum(ims1.*ims1);
m2j = sum((-imzr1).*(-imzr1));
m11 = sum(ims1.*(-imzr1));
c22 = (m22-m2i*m2j-2*m11*m11)/(m2i*m2j);