% Nearest Feature Plane Classifier-NFP
function [NFPCrate]=NFPclassifier(features,test_features,trnum,tenum,classnum)

% features       the matrix that training samples projected on feature subspace(特征維數*訓練樣本數，列矢量)
% test_features  the matrix that test samples projected on feature subspace(特征維數*測試樣本數，列矢量) 
% trnum          the number of training samples of each class
% tenum          the number of test samples of each class
% classnum       the number of classes

% NFPCrate       the output correct classification rate of each class and the total
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% decompose the features and test_features according to the class flag
for i=1:classnum
    subd(:,:,i) =features(:,(i-1)*trnum+1:i*trnum);      % 訓練樣本集
    subdd(:,:,i)=test_features(:,(i-1)*tenum+1:i*tenum); % 測試樣本集
end

% number of NFP of each class
total=trnum*(trnum-1)*(trnum-2)/6;
% compute the parameter U (斜率) and the position P (垂點)
for m=1:classnum                     % 類別數
    for t=1:tenum                    % 測試樣本數/類
        k=1;        
        for c=1:classnum             % 類別數
            for i=1:trnum            % 訓練樣本數/類
                for j=i+1:trnum      % 訓練樣本數/類
                    for h=j+1:trnum  % 訓練樣本數/類
                        Wijh=[subd(:,i,c),subd(:,j,c),subd(:,h,c)];
                        u=size(Wijh,2);
                        P=Wijh*inv(Wijh'*Wijh+(1e-3)*eye(u))*Wijh'*subdd(:,t,m);  % 垂足                        
                        dis(t,k,m)=norm(subdd(:,t,m)-P);            % 計算測試點到垂足的距離
                        k=k+1;       % 計算k條特征線(NFP)
                    end
                end
            end
        end       
    end    
end

% compute the correct recognition rate and the according mean value
CCrate=zeros(1,classnum);
for c=1:classnum
    dist=dis(:,:,c)';
    [Y,I]=min(dist);
    k=0;
    for i=1:tenum
        if ceil(I(i)/total)==c
            k=k+1;
        end
    end
    CCrate(c)=k/tenum;            
end
rtmean=sum(CCrate)/classnum;

% output
NFPCrate=zeros(1,classnum+1);
NFPCrate=[CCrate,rtmean];