function CS6640_Week11
%

rng('default');

% A6: Shape params
paper1 = imread('paper1.jpg');
paper2 = imread('paper2.jpg');
paper1g = rgb2gray(paper1);
paper1t = paper1g>150;
paper1ed = bwmorph(bwmorph(paper1t,'erode',2),'dilate',2);
paper2g = rgb2gray(paper2);
paper2t = paper1g>130;
paper2ed = bwmorph(bwmorph(paper2t,'erode',2),'dilate',2);
stats1 = regionprops(paper1ed,'all');
stats2 = regionprops(paper2ed,'all');
[stats1.Perimeter,stats2.Perimeter,stats1.PerimeterOld,stats2.PerimeterOld]
stats1CI = regionprops(stats1.ConvexImage,'all')

% Explore dimensionality
pts1 = randn(20,2);
pts2 = randn(20,2) + [10*ones(20,1),zeros(20,1)];
figure(1)
clf
plot(pts1(:,1),pts1(:,2),'ko');
hold on
plot(-5,-10,'w.');
plot(12,10,'w.');
plot(pts2(:,1),pts2(:,2),'bo');
pts = [pts1;pts2];
[cidx,ctrs] = kmeans(pts,2);
plot(ctrs(:,1),ctrs(:,2),'r*');
pts10 = [pts,zeros(40,8)];
[cidx10,ctrs10] = kmeans(pts10,2);
ctrs10

% Use eigen stuff to reduce dimensionality
C = cov(pts10)
[V,D] = eig(C)
ptsr = (V'*pts10')';
ptsd = [ptsr(:,10),ptsr(:,9)];
plot(-ptsd(:,1),-ptsd(:,2),'r*');

% Create some points
pts = [-5,-5; -3,-3; -1,-1; 1,1; 3,3; 5,5]
num_pts = length(pts(:,1));
figure(1);
clf
plot(pts(:,1),pts(:,2),'ko');
hold on
axis equal
plot(-10,-10,'w.');
plot(10,10,'w.');
C = cov(pts)

% Add noise
ptsn = pts + 0.5*randn(6,2)
ptsnm = mean(ptsn)
pts0 = ptsn;
pts0(:,1) = pts0(:,1) - ptsnm(1);
pts0(:,2) = pts0(:,2) - ptsnm(2);
plot(pts0(:,1),pts0(:,2),'r*');
C0 = cov(pts0)

% find first principal axis
[p1,s1] = CV_total_LS(pts0(:,1),pts0(:,2))
line1 = [0,0,0; -p1(2) p1(1) 0];

% find projection of points onto line
p_pts = zeros(num_pts,2);
n_pts = zeros(num_pts,2);
dists = zeros(num_pts,1);
for k = 1:num_pts
    [dd,pp] = CS6640_cv_signed_dist_pt_line([ptsn(k,:) 0],line1);
    dists(k) = dd;
    p_pts(k,:) = pp(1:2);
    n_pts(k,:) = pts0(k,:) - p_pts(k,:);
end
plot(n_pts(:,1),n_pts(:,2),'k*');
plot(n_pts(:,1),n_pts(:,2),'go');

% find second principal axis
[p2,s2] = CV_total_LS(n_pts(:,1),n_pts(:,2))
x = [-5:0.1:5];
y1 = (-p1(3)-p1(1)*x)./p1(2);
y2 = (-p2(3)-p2(1)*x)./p2(2);
plot(x,y1,'k');
plot(x,y2,'k');

% Find principal axes as eigenectors
[V0,D0] = eigs(C0)
a3 = V0(1,1);
b3 = V0(2,1);
a4 = V0(1,2);
b4 = V0(2,2);
y3 = -a3*x/b3;
y4 = -a4*x/b4;
plot(x,y3,'r');
plot(x,y4,'r');

% Eigenvector matrix is a rotation
ptsr = (V0'*pts0')';
plot(ptsr(:,1),ptsr(:,2),'k+');

% from summary pp. 253-254

% principal axes are mutually orthogonal
dot(V0(:,1),V0(:,2))

% covariance matrix of new points uncorrelated
Cr = cov(ptsr)

% principal axes maximize variance
v = zeros(360,1);
for t = 0:259
    T = [cosd(t),sind(t); -sind(t),cosd(t)];
    ptst = (T*pts0')';
    v(t+1) = var(ptst(:,1));
end
figure(3);
clf
plot(v);
[vals,indexes] = sort(v,'descend');
[indexes(1),vals(1)]

% eigenvalues directly specify variance
D0 % note: max(v) above equals max(diag(D0))

% A6: PCA kmeans
% Dimension reduction
C = cov(pts)
[V,D] = eigs(C)
ptsr = (V'*pts')'
figure(1)
clf
plot(pts(:,1),pts(:,2),'ko');
hold on
plot(-8,-8,'w.');
plot(8,8,'w.');
plot(ptsr(:,1),ptsr(:,2),'r*');

% Reduce dimensions
pts_red = ptsr(:,1)

map1 = imread('map1.jpg');
im1 = map1(201:500,601:900,:);
im1g = rgb2gray(im1);
T = CS6640_FFT_texture(double(im1g));
size(T)

% Example: use 10 measures twice
T1 = T(:,16:25);
T2 = [T1,T1];
C = cov(T2);
[V,D] = eig(C);
indexes = find(diag(D)<0.1)

L = CS6640_Laws(double(im1g));
Lr = reshape(L,300*300,16);
size(Lr)
RGB = reshape(im1,300*300,3);
size(RGB)
TT = [T,Lr,RGB];
size(TT)
[cidx1,ctrs1] = kmeans(double(TT),6);
[cidx2,ctrs2] = CS6640_PCA_kmeans(double(TT),0.5,6);

imc1 = reshape(cidx1,300,300);
imc2 = reshape(cidx2,300,300);
figure(1)
clf
subplot(2,3,1);
combo(im1g,imc1==1);
subplot(2,3,2);
combo(im1g,imc1==2);
subplot(2,3,3);
combo(im1g,imc1==3);
subplot(2,3,4);
combo(im1g,imc1==4);
subplot(2,3,5);
combo(im1g,imc1==5);
subplot(2,3,6);
combo(im1g,imc1==6);

figure(2)
clf
subplot(2,3,1);
combo(im1g,imc2==1);
subplot(2,3,2);
combo(im1g,imc2==2);
subplot(2,3,3);
combo(im1g,imc2==3);
subplot(2,3,4);
combo(im1g,imc2==4);
subplot(2,3,5);
combo(im1g,imc2==5);
subplot(2,3,6);
combo(im1g,imc2==6);

size(ctrs1)
size(ctrs2)

tch = 0;
