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

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

s = 3;
[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:3)*S(2:3,2:3);
b = diag(1./V(:,1))*V(:,2:3);

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

% computation of the inner polygons
IPS = a(cips,:);
IPC = b(cipc,:);

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

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

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

if exist('addlines') == 1 && addlines == 1
    ylin = 5*yls;
    for i=1:1:size(D,2)
        xlin(i,:) = -V(i,1)/V(i,2)-ylin*V(i,3)/V(i,2);
        if i==1
            plot(xlin(i,:), ylin,'color',[.7 .7 .7]);
        else
            line(xlin(i,:), ylin,'color',[.7 .7 .7]);
        end
    end
    line(FPS(:,1), FPS(:,2),'color',col(1,:), 'linewidth', 1)
    line(IPS(:,1), IPS(:,2),'color',col(2,:), 'linewidth', 1)
    line(a(:,1), a(:,2),'color',col(2,:), 'linewidth', 0.5, 'Marker','+', 'Markersize',5, 'linestyle', 'none')
else
    plot(FPS(:,1), FPS(:,2),'color',col(1,:), 'linewidth', 1)
    line(IPS(:,1), IPS(:,2),'color',col(2,:), 'linewidth', 1)
end
title('outer and inner polygon spectral factor','FontWeight','Normal')
set(gca,'xlim',1.05*xls,'ylim',1.05*yls)

% plot outer polygon and inner polygon for C
figure
xlc = [min(FPC(:,1)) max(FPC(:,1))];
ylc = [min(FPC(:,2)) max(FPC(:,2))];
if exist('addlines') == 1 && addlines == 1
    ylin = 5*ylc;
    for i=1:1:size(D,1)
        xlin(i,:) = -U(i,1)*S(1,1)/U(i,2)/S(2,2)-ylin*U(i,3)*S(3,3)/U(i,2)/S(2,2);
        if i==1
            plot(xlin(i,:), ylin,'color',[.7 .7 .7]);
        else
            line(xlin(i,:), ylin,'color',[.7 .7 .7]);
        end
    end
    line(FPC(:,1), FPC(:,2),'color',col(1,:), 'linewidth', 1)
    line(IPC(:,1), IPC(:,2),'color',col(2,:), 'linewidth', 1)
    line(b(:,1), b(:,2),'color',col(2,:), 'linewidth', 0.5, 'Marker','+', 'Markersize',5, 'linestyle', 'none')
else
    plot(FPC(:,1), FPC(:,2),'color',col(1,:), 'linewidth', 1)
    line(IPC(:,1), IPC(:,2),'color',col(2,:), 'linewidth', 1)
end
title('outer and inner polygon concentrational factor','FontWeight','Normal')
set(gca,'xlim',1.05*xlc,'ylim',1.05*ylc)
