function fpip3D(D, addpoints)
% plots the outer polygons and the inner polygons for a data matrix D with
% respect to s = 4 components
% 
% input:
% D ....... (k x n)-matrix, with
%      k is the number of spectra, 
%      n is the number of channels
% 
% optional input:
% addlines ...... use "addpoints = 1" to plot also the limiting lines,
%                 use "addpoints = 0" to skip them
% 
% examples:
% fpip3D(D)
% fpip3D(D,1)

col = [0     0.4470     0.7410;
    0.8500    0.3250    0.0980];

COL = [0.85     0.9     1;
    0.995    0.8250    0.75];

s = 4;
[U,S,V] = svds(D,s);
for i=1:s
    if min(V(:,i)) < -max(V(:,i))
        V(:,i) = -V(:,i);
        U(:,i) = -U(:,i);
    end
end
a = diag(1./U(:,1)/S(1,1))*U(:,2:4)*S(2:4,2:4);
b = diag(1./V(:,1))*V(:,2:4);

% computation of the indices of a resp. b that define the vertices of the inner polygons
kips = convhull(a);
kipc = convhull(b);

% computation of outer polygon (FIRPOL) for S using duality and IPC
FPS = [];
m = size(kipc,1);
for i=1:m
    X = b(kipc(i,:),:);
    FPS(end+1,:) = X\[-1;-1;-1];
end
kfps = convhull(FPS);

% computation of outer polygon (FIRPOL) for C using duality and IPS
FPC = [];
m = size(kips,1);
for i=1:m
    X = a(kips(i,:),:);
    FPC(end+1,:) = X\[-1;-1;-1];
end
kfpc = convhull(FPC);

% plot outer polygon and inner polygon for S
figure
xls = [min(FPS(:,1)) max(FPS(:,1))];
yls = [min(FPS(:,2)) max(FPS(:,2))];
zls = [min(FPS(:,3)) max(FPS(:,3))];

trisurf(kfps, FPS(:,1), FPS(:,2), FPS(:,3), 'facecolor', col(1,:), 'facealpha',0.2);
hold on;
trisurf(kips, a(:,1), a(:,2), a(:,3), 'facecolor', col(2,:), 'facealpha',0.7);
if exist('addpoints') == 1 && addpoints == 1
    line(a(:,1), a(:,2), a(:,3), 'color',col(2,:), 'linewidth', 1, 'Marker','+', 'Markersize',10, 'linestyle', 'none')
end
hold off;
title('outer and inner polygon spectral factor','FontWeight','Normal')
set(gca,'xlim',1.05*xls,'ylim',1.05*yls,'zlim',1.05*zls)

% plot outer polygon and inner polygon for C
figure
xlc = [min(FPC(:,1)) max(FPC(:,1))];
ylc = [min(FPC(:,2)) max(FPC(:,2))];
zlc = [min(FPC(:,3)) max(FPC(:,3))];
trisurf(kfpc, FPC(:,1), FPC(:,2), FPC(:,3), 'facecolor', col(1,:), 'facealpha',0.2);
hold on;
trisurf(kipc, b(:,1), b(:,2), b(:,3), 'facecolor', col(2,:), 'facealpha',0.7);
if exist('addpoints') == 1 && addpoints == 1
    line(b(:,1), b(:,2), b(:,3), 'color',col(2,:), 'linewidth', 1, 'Marker','+', 'Markersize',10, 'linestyle', 'none')
end
hold off;
title('outer and inner polygon concentrational factor','FontWeight','Normal')
set(gca,'xlim',1.05*xlc,'ylim',1.05*ylc,'zlim',1.05*zlc)
