function simcaplot(result); 

%SIMCAPLOT plots a scatter plot with the boundaries defined by the SIMCA method. 
% It is based on the results from a simca analysis (see rsimca.m or csimca.m).
% 
% For technical reasons, 6 different groups can be plotted (with different symbols).
% In case there are more groups, please adapt lines 15--17.
%
% I/O: simcaplot(result)
%
% Written by K. Vanden Branden on 01/08/04
% Last update: 21/09/2004

p = size(result.x,2);
markings={'bx';'ro';'kdiamond';'y+';'g.';'m*'};
linecolor={'b-','r-','k-','y-','g-','m-'};
colorellipse={'b','r','k','y','g','m'};
if p>3
    error('The dimension of the dataset is larger than 3.')
elseif length(result.avemisclasprob)>6
    error(['Only 6 groups can be drawn with different symbols. Please adapt the code',... 
            ' of simcaplot.m if you want to draw more groups.'])
else
    legstr=[];
    nClass = length(result.groupmisclasprob);
    for iClass = 1:nClass
        out = result.pca{iClass};
        groupi = result.x(find(result.group == iClass),:);
        %legstr=[legstr; sprintf(['Group',num2str(iClass)])]; 
        if out.k == 1
            [out.lbord,out.rbord]=confreg(out);
            bord=[out.lbord,out.rbord];
        elseif out.k == 2
            bord=ellipse([0 0],diag(out.L(1:2)));
        elseif out.k == 3
            bord=ellips3D(out.M, out.P*diag(out.L)*out.P');
        end
        if p == 2
            xlabel('x');ylabel('y');
            if size(bord) == [1,4]
                plot(groupi(:,1),groupi(:,2),markings{iClass});hold on;
                plot([bord(1,1),bord(1,3)],[bord(1,2),bord(1,4)],linecolor{iClass});
            elseif size(bord)==[802,2]
                plot(groupi(:,1),groupi(:,2),markings{iClass});hold on;
                newbord=bord*out.P'+repmat(out.M,802,1);
                line(newbord(:,1),newbord(:,2));
                
            end
        elseif p == 3
            xlabel('x');ylabel('y');zlabel('z')
            if size(bord) == [1,6]
                plot3(groupi(:,1),groupi(:,2),groupi(:,3),markings{iClass});hold on;grid on;
                plot3([bord(1,1),bord(1,4)],[bord(1,2),bord(1,5)],[bord(1,3),bord(1,6)],linecolor{iClass})
            elseif size(bord)==[802,2]
                plot3(groupi(:,1),groupi(:,2),groupi(:,3),markings{iClass});hold on;grid on;
                newbord=bord*out.P'+repmat(out.M,802,1);
                plot3(newbord(:,1),newbord(:,2),newbord(:,3),linecolor{iClass});
            else 
                plot3(groupi(:,1),groupi(:,2),groupi(:,3),markings{iClass});hold on;grid on;
                mesh(bord.x,bord.y1,bord.z1,'Edgecolor',colorellipse{iClass});alpha(0.2);
            end
        end
    end
    title(result.class)
    %legend(legstr,0)
    hold off
end


%------------------
function [lbord,rbord,outl]=confreg(result);

n1=size(result.T,1);p=size(result.P,1);
rbord=result.M+sqrt(result.L)*result.cutoff.sd*result.P';
lbord=result.M-sqrt(result.L)*result.cutoff.sd*result.P';

%----------------------
function coord=ellipse(mean,covar)

% Determines the coordinates of some points that lie on the 97.5 % tolerance ellipse.

deter=covar(1,1)*covar(2,2)-covar(1,2)^2;
ylimit=sqrt(7.37776*covar(2,2));
y=-ylimit:0.005*ylimit:ylimit;
sqtdi=sqrt(deter*(ylimit^2-y.^2))/covar(2,2);
sqtdi([1,end])=0;
b=mean(1)+covar(1,2)/covar(2,2)*y;
x1=b-sqtdi;
x2=b+sqtdi;
y=mean(2)+y;
coord=[x1,x2([end:-1:1]);y,y([end:-1:1])]';

%----------------
function out = ellips3D(mu, sigma)

%This function determines the 3D ellipsoid.

dist = sqrt(chi2inv(0.975,3));

A = inv(sigma);
[ev,ew] = eig(A);
A = ev'*A*ev;

c1 = 0.1;
zlimit = dist/sqrt(A(3,3));
z = -zlimit:(c1*zlimit):zlimit;
c2 = (1/c1*2 + 1)*2; % is length(z)*2
c3 = c2/2; % is  length(z)
z1 = ones(1,c2)*z(1);
x1 = zeros(c3,c3);
x2 = zeros(c3,c3);
y = zeros(c3,c3);
y1 = zeros(c3,c2);

for i = 2:(c3-1)
    z1 = [z1 ones(1,c2)*z(i)];
    D = dist^2 - A(3,3)*z(i)^2;
    ylimit = sqrt(D)/sqrt(A(2,2));
    y(i,:) = -ylimit:(c1*ylimit):ylimit;
    y1(i,:) = [y(i,:) ylimit:(-c1*ylimit):-ylimit];
    
    for j = 2:(c3-1)
        x1(i,j) = sqrt((D - A(2,2)*y(i,j)^2)/A(1,1));
        x2(i,j) = -sqrt((D - A(2,2)*y(i,j)^2)/A(1,1));
    end
end

x = [x1,x2];
z1 = [z1 ones(1,c2)*z(c3)];
z1 = vec2mat(z1,c2);
grotex = vec2mat([mat2vec(x) mat2vec(y1) mat2vec(z1)],c2*c3);
groottex = ev*grotex;

groottex(1,:) = groottex(1,:) + mu(1);
groottex(2,:) = groottex(2,:) + mu(2);
groottex(3,:) = groottex(3,:) + mu(3);

x = vec2mat(groottex(1,:),c2);
y1 = vec2mat(groottex(2,:),c2);
z1 = vec2mat(groottex(3,:),c2);
out.x=x;
out.y1=y1;
out.z1=z1;


%-----------------
function vec = mat2vec(mat)

nkolom = size(mat,2);
nrij = size(mat,1);

vec = [];

for i=1:nkolom
   hulpvec = wkeep(mat,[nrij,1],[1,1]);
   vec = [vec hulpvec'];
   mat = mat(:,2:(nkolom - (i-1)));
end

%----------------------

function mat = vec2mat(vec,ncol)

nrow = length(vec)/ncol;

for i = 1:nrow
   mat(i,:)  = wkeep(vec,ncol,'l');
   vec = vec((ncol + 1):length(vec));
end