clear
clc
close all;
%m1,m2均為10個樣本
% m1=[-0.4,0.58,0.089;-0.31,0.27,-0.04;0.38,0.055,-0.035;-0.15,0.53,0.011;-0.35,0.47,0.034;0.17,0.69,0.1;-0.011,0.55,-0.18];
% m2=[0.83,1.6,-0.014;1.1,1.6,0.48;-0.44,-0.41,0.32;0.047,-0.45,1.4;0.28,0.35,3.1;-0.39,-0.48,0.11;0.34,-0.079,0.14];
% a) 對類別m1,m2,利用fisher原則計算投影方向w
m1=[0 0 0; 1 0 0; 1 0 2; 1 1 0; 0.5 0 0.5;1 0.5 0;1 0 0.5];
m2=[0 0 1; 0 1 0; 0 1 1; 1 2 1;0 0.5 1;0 1 0.5; 0.5 1 2];
u1=mean(m1);%求均值
u2=mean(m2);
%計算類內散度Si和總類內散度Sw
S1=(m1(1,:)-u1)'*(m1(1,:)-u1);
S2=(m2(1,:)-u2)'*(m2(1,:)-u2);
for i=2:7
    S1=S1+(m1(i,:)-u1)'*(m1(i,:)-u1);
    S2=S2+(m2(i,:)-u2)'*(m2(i,:)-u2);
end
Sw=S1+S2;
%變換向量
w=inv(Sw)*(u1-u2)';
% b) 畫出投影前的兩類樣本點和投影線.(即特征向量：各類方差)
x1=m1(1:7,1);y1=m1(1:7,2);z1=m1(1:7,3);
x2=m2(1:7,1);y2=m2(1:7,2);z2=m2(1:7,3);
subplot(1,2,1);plot3(x1,y1,z1,'r*',x2,y2,z2,'b*'); %plot3為三維線圖
title('原樣本分布圖');%畫出原兩類樣本點
hold on
x=-2:0.1:2;
y=w(2,1)/w(1,1)*x;
z=w(3,1)/w(1,1)*x;
plot3(x,y,z);%畫出投影線
legend('第一類樣本','第二類樣本','投影線',2);
%rotate3D on       %鼠標移動可以旋轉三維坐標

% c) 畫出投影后的兩類樣本點和投影線
subplot(1,2,2);
x=-2:0.1:2;
y=w(2,1)/w(1,1)*x;
z=w(3,1)/w(1,1)*x;
plot3(x,y,z);%畫出投影線
hold on
% 對w1中的點投影
for i=1:7
    x1=(m1(i,1)+w(2,1)/w(1,1)*m1(i,2)+w(3,1)/w(1,1)*m1(i,3))/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
    y1=w(2,1)/w(1,1)*x1;
    z1=w(3,1)/w(1,1)*x1;
    plot3(x1,y1,z1,'r*');  
    X1(i,1)=x1;
    X1(i,2)=y1;
    X1(i,3)=z1;
end

% 對w2中的點投影
for i=1:7
    x2=(m2(i,1)+w(2,1)/w(1,1)*m2(i,2)+w(3,1)/w(1,1)*m2(i,3))/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
    y2=w(2,1)/w(1,1)*x2;
    z2=w(3,1)/w(1,1)*x2;
    plot3(x2,y2,z2,'b*');  title('投影后的樣本分布圖');
    X2(i,1)=x2;
    X2(i,2)=y2;
    X2(i,3)=z2;
end

for i=1:7
    X11(i)=w'*X1(i,:)';
    X21(i)=w'*X2(i,:)';
end
M1=mean(X11);
M2=mean(X21);
F=(M1+M2)/2;%分界閾值點
a=input('請輸入測試樣本點坐標  x=');
b=input('請輸入測試樣本點坐標  y=');
c=input('請輸入測試樣本點坐標  z=');
plot3(a,b,c,'go');
a1=(a+w(2,1)/w(1,1)*b+w(3,1)/w(1,1)*c)/(1+(w(2,1)/w(1,1))^2+(w(3,1)/w(1,1))^2);
b1=w(2,1)/w(1,1)*a1;
c1=w(3,1)/w(1,1)*a1;
L=[a1 b1 c1];
exam=w'*L';
if exam>F
   plot3(a1,b1,c1,'go'); 
   fprintf('該點屬于第一類\n');
else
   plot3(a1,b1,c1,'go');
   fprintf('該點屬于第二類\n');
end
rotate3D on    %鼠標移動可以旋轉三維坐標