%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PCA, DCT and LDA % Z. Li, 2009.02.08 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mean, variance, covariance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n=1000; x0=randn(n, 2) + repmat([3 5],n, 1); A0 = [2.0 1.0; 1.0 0.8]; x = x0*A0; % plot plot(x(:,1), x(:, 2), '.'); hold on; xlabel('x1'); ylabel('x2'); grid on; % compute mean: mx = mean(x); plot(mx(1), mx(2), '+r'); % compute variance: v1 = var(x(:,1)); v2=var(x(:,2)); str=sprintf('\n variances = [%1.2f %1.2f]', v1, v2); title(str); % covariance of x: n x 2: % cv=cov(x) : 2x2 matrix: % cv(1,1) = var(x1); cv(2,2)=var(x2); % cv(1,2) = cv(2,1) cov(x) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % pca of x %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % remove mean: x = x - repmat(mx, n, 1); [A1, L1]=eig(cov(x)); a1 = A1(1,:); a2 = A1(2,:); fprintf('\n eigen values:'); fprintf('%1.2f ', diag(L1)); px = [0 a1(1)]; py=[0, a1(2)]; plot(px, py, '-m'); hold on; px = [0 a2(1)]; py=[0, a2(2)]; plot(px, py, '-r'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Eigenface %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% facepath='D:/zli/teaching/comp435/matlab/data/att_faces/'; nSubj=40; nPic=10; % icon size h=24; w=20; %count t=1; for j=1:nSubj for k=1:nPic fim_name=sprintf('%ss%d/%d.pgm', facepath, j, k); im = imread(fim_name); icn = imresize(im, [h, w], 'bilinear'); faces(t,1:w*h) = (icn(:))'; t=t+1; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % pca of faces %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% faces = double(faces); mfaces = mean(faces); mfaces2d = reshape(mfaces, h, w); S = cov(faces); [EigFaces, L]=eig(S); [EigVal, indx] = sort(diag(L), 'descend'); figure; stem(EigVal(1:200), '.'); figure; m = 8; for k=1:m subplot(2, 4,k); eigf = EigFaces(:, indx(k)); eigf2 = reshape(eigf, h, w); if (0) colormap('gray'); imagesc(eigf2); axis off; else stem(eigf, '.'); end str=sprintf('f_%d', k'); title(str); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % projection examples %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % projection x = faces*EigFaces(:, indx([1:3])); subj = [5 12 20 32]; pics =[1: 5]; styls{1}='+r'; styls{2}='*r'; styls{3}='om'; styls{4}='xm'; figure; plot(x(:,1), x(:,2), '.'); hold on; for k=1:length(subj) offs = 10*(subj(k)-1); plot(x(offs+1:offs+5, 1), x(offs+1:offs+5, 2), styls{k}); end figure; for j=1:length(subj) for k=1:length(pics) fim_name=sprintf('%ss%d/%d.pgm', facepath, subj(j), 1); f = imread(fim_name); subplot(2,2,j); imshow(f); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DCT projection examples %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure; for k=1:4 subplot(4,1,k); stem(faces(k, 1:16), '.'); end figure; for k=1:4 subplot(4,1,k); stem(dct(faces(k, 1:16), 16), '.'); end figure; im1 = imread('lena.bmp'); im1 = imresize(im1, [24, 24], 'bilinear'); im1 = double(im1); im1dct = dct2(im1, 24, 24); subplot(1,2,1); colormap('gray'); imagesc(im1); subplot(1,2,2); colormap('gray'); imagesc(im1dct); return; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % function myLDA() % z.li, 10-20-2006 % Fisher Linear Discriminant Analysis % function dependency: % - n/a % input: % X - n x d data % Lbl - n x 1, lables % A - subspace % Y - projection of Y % minCnt - use only labels with minCnt members % output: % A - local subspace % nClass - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %function [A, nClass]=myLDA(X, Lbl) %function [A, nClass]=myLDA(X, Lbl) dbg='n'; if dbg=='y' X = [ 60 61 65 92 95 86 82 70 79 216 220 222 62 55 54 71 54 54 179 207 191 211 213 211 87 82 85 71 62 77 80 164 208 218 219 219 50 49 51 49 48 60 78 175 172 189 194 195 105 110 79 84 75 83 96 106 172 217 218 218 60 74 77 124 143 205 208 208 207 208 209 211 103 126 108 96 93 93 105 129 209 227 226 226 104 134 159 198 194 155 216 217 218 219 218 219 66 72 96 124 185 184 207 207 207 207 208 209 78 58 60 51 184 144 75 172 202 206 205 204 67 67 56 70 81 55 71 86 60 87 34 200 76 88 162 180 194 195 198 197 197 198 199 200 72 49 49 48 47 47 48 47 59 62 59 171 95 120 112 147 116 112 133 147 156 220 227 227 90 80 99 189 206 210 210 211 213 213 213 213 75 72 69 57 46 45 47 54 46 49 50 112 85 69 72 67 57 56 52 54 46 49 50 58 54 48 50 173 183 175 176 176 177 177 177 177 64 69 54 51 49 100 176 179 178 179 179 179 74 69 82 90 101 149 128 205 220 228 226 218 116 180 189 197 186 194 198 197 197 198 198 199 184 188 188 189 191 192 192 191 194 194 194 195 71 65 58 71 60 122 178 179 180 177 181 182 197 199 199 200 199 200 201 201 201 200 200 200 71 81 142 176 183 184 185 187 188 189 190 191 33 45 35 34 34 49 86 83 83 81 81 80 84 66 66 68 63 96 97 72 171 226 227 228 53 61 68 65 77 100 135 194 205 205 205 207 99 126 160 122 162 184 195 197 197 198 198 199 93 84 78 80 86 139 117 158 226 220 222 221 ]; Lbl = [ 1 1 13 8 14 14 14 14 9 1 9 9 9 13 13 8 1 9 13 14 1 8 14 8 8 14 13 8 14 13 ]; end %function [A, nClass]=myLDA(X, Lbl) % var [n kDim]=size(X); uLabels = sort(unique(Lbl)); nClass = length(uLabels); % LDA mX = mean(X); mm = zeros(kDim, kDim); % compute between class variance Sb and within class variance Sw for k=1:nClass indx = find(Lbl==uLabels(k)); mXk = mean(X(indx,:)); % m x kDim mm = mm + length(indx)*mXk'*mXk; end Sw = X'*X - mm; Sb = mm - n*mX'*mX; Sw = (Sw+Sw')/2; Sb = (Sb+Sb')/2; option = struct('disp',0); [A, egv]=eigs(Sb, Sw, nClass-1, 'la', option); for k=1:nClass-1 A(:,k) = A(:,k)./norm(A(:,k)); end % proj return;