clear?all
clc
close?all

um=100;
ustd=80;

train_sample=[];%所有訓練圖像
ClassNum=40;?%?40?persons
tol_num=10;
train_samplesize=5;?%?Num?of?pict?for?training
train=1:train_samplesize;
test=train_samplesize+1:tol_num;
M=200;
p=40;?%%PCA降維維數
%每類訓練樣本、測試樣本個數
train_num=length（train）;
test_num=length（test）;
rows=112;
cols=92;
%讀入訓練樣本，測試樣本
[train_sampletrain_label]=readsample（ClassNumtrainrowscols）;
[test_sampletest_label]=readsample（ClassNumtestrowscols）;
icol=cols;
irow=rows;
%?read?and?show?image
S=[];????%?img?matrix
figure（1）;
for?i=1:40
????for?j=1:5
str=strcat（‘[ORL]/s‘num2str（i）‘/‘num2str（j）‘.bmp‘）;????%?concatenates?two?strings?that?form?the?name?of?the?image
eval（‘img=imread（str）;‘）;
%?subplot（ceil（sqrt（M））ceil（sqrt（M））i）
%?imshow（img）
%?if?i==3
%?title（‘Training?set‘‘fontsize‘18）
%?end
%?drawnow;
[irow?icol]=size（img）;????%?get?the?number?of?rows?（N1）?and?columns?（N2）
temp=reshape（img‘irow*icol1）;????%?creates?a?（N1*N2）x1?vector
S=[S?temp];?%?S?is?a?N1*N2xM?matrix?after?finishing?the?sequence
????end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%?Here?we?change?the?mean?and?std?of?all?images.?We?normalize?all?images.
%?This?is?done?to?reduce?the?error?due?to?lighting?conditions?and?background.
for?i=1:size（S2）
temp=double（S（:i））;
m=mean（temp）;
st=std（temp）;
S（:i）=（temp-m）*ustd/st+um;%
end

%?show?normalized?images
figure（2）;
for?i=1:M
str=strcat（int2str（i）‘.jpg‘）;
img=reshape（S（:i）icolirow）;
img=img‘;
eval（‘imwrite（imgstr）‘）;?
%?subplot（ceil（sqrt（M））ceil（sqrt（M））i）
%?imshow（img）
%?drawnow;
if?i==3
title（‘Normalized?Training?Set‘‘fontsize‘18）
end
end


%?mean?image
m=mean（S2）;??%?obtains?the?mean?of?each?row?instead?of?each?column
tmimg=uint8（m）;?%?converts?to?unsigned?8-bit?integer.?Values?range?from?0?to?255
img=reshape（tmimgicolirow）;?%?takes?the?N1*N2x1?vector?and?creates?a?N1xN2?matrix
img=img‘;?
figure（3）;
imshow（img）;
title（‘Mean?Image‘‘fontsize‘18）

%?Change?image?for?manipulation
dbx=[];????%?A?matrix
for?i=1:M
temp=double（S（:i））;
dbx=[dbx?temp];
end

%Covariance?matrix?C=A‘A?L=AA‘
A=dbx‘;
L=A*A‘;
%?vv?are?the?eigenvector?for?L
%?dd?are?the?eigenvalue?for?both?L=dbx‘*dbx?and?C=dbx*dbx‘;
[vv?dd]=eig（L）;
%?Sort?and?eliminate?those?whose?eigenvalue?is?zero
v=[];
d=[];
for?i=1:size（vv2）
if（dd（ii）>1e-4）
v=[v?vv（:i）];
d=[d?dd（ii）];
end
end

%sort?will?return?an?ascending?sequence
[B?index]=sort（d）;
ind=zeros（size（index））;
dtemp=zeros（size（index））;
vtemp=zeros（size（v））;
len=length（index）;
for?i=1:len
dtemp（i）=B（len+1-i）;
ind（i）=len+1-index（i）;
vtemp（:ind（i））=v（:i）;
end
d=dtemp;
v=vtemp;


%Normalization?of?eigenvectors
for?i=1:size（v2）?%access?each?column
kk=v（:i）;
temp=sqrt（sum（kk.^2））;
v（:i）=v（:i）./temp;
end

%Eigenvectors?of?C?matrix
u=[];
for?i=1:p?%size（v2）
temp=sqrt（d（i））;
u=[u?（dbx*v（:i））./temp];
end

%Normalization?of?eigenvectors
for?i=1:p%size（u2）
kk=u（:i）;
temp=sqrt（sum（kk.^2））;
u（:i）=u（:i）./temp;
end


%?show?eigenfaces
figure（4）;
for?i=1:p%size（u2）
img=??reshape（u（:i）icolirow）;
img