
function [A,b]=determine_facets_planes(vertices,facets)

%--------------------------------------------------------------------------
% Object:
% Determining convex polyhedron equations A*u <= b.
% Each facet suitably determines one of these equations.
%
% A convex polyhedron has points "u=[x; y; z]" that satisfy the equation
% "A*u <= b" for a matrix "A" and a vector "b". We determine such "A" and
% "b", knowing that the k-th facet lies on a plane with equation
% "A(k,:)*u=b(k)".
%
% As procedure we
% 1. compute the equation of each facet;
% 2. determine an other "test" facet;
% 3. see if any the vertices of the test facet satisfy "A(k,:)*u < b(k)"
%    and if not modify "A(k,:)" into "-A(k,:)", "b(k)" into "b(k)".
%--------------------------------------------------------------------------
% Input:
% 
% vertices: The vertices are stored in "vertices" and is an "m x 3" matrix.
%  The k-th vertex has coordinates "vertices(k,:)".
%
% facets: The variable "facets" can be:
%
% 1. A "cell array":
%    In this case the i-th cell is a vector of positive integers,
%    representing the indices of vertices of the i-th facet, taking into
%    account the variable "vertices". The orientation of the facet is
%    counterclockwise (w.r.t to outward normal).
%
% 2. A "matrix" with dimension "m x 3":
%    In this case the i-th row is a vector of positive integers,
%    representing the indices of vertices of the triangular facet, taking
%    into account the variable "vertices". The orientation of the
%    triangular facet is counterclockwise (w.r.t to outward normal).
%--------------------------------------------------------------------------
% Output:
% 
% A,b: matrix and vector so that the k-th facet is contained in the plane
%     whose points X=(x(1); x(2); x(3)) satisfy the equation 
%
%                              A(k,:)*X <= b(k)
%--------------------------------------------------------------------------

tf = iscell(facets);
number_of_facets=length(facets);
A=[]; b=[];

for k=1:number_of_facets
    
    % Row vector of indices of vertices describing the k-th facet.
    kt=setdiff(1:number_of_facets,k);
    
    if tf == 1 % facets are stored in cells
        facetL=facets{k}; facet_test_L=facets{kt};
    else % facets are stored in m x 3 matrix (triangular facets)
        facetL=facets(k,:); 
        facet_test_L=facets(kt,:);
    end
    
    % Vertices of the facet ordered counterclockwise w.r.t. the outward
    % normal.
    XV=vertices(facetL',1); YV=vertices(facetL',2); ZV=vertices(facetL',3);
    
    % ..... Computing plane equation ..... 
    
    P1=[XV(1) YV(1) ZV(1)]; P1=P1';
    P2=[XV(2) YV(2) ZV(2)]; P2=P2';
    P3=[XV(3) YV(3) ZV(3)]; P3=P3';
    
    % Equation of the plane containing the facets.
    % See e.g.:
    % https://math.stackexchange.com/questions/2686606/...
    % equation-of-a-plane-passing-through-3-points
    %
    % Here "a" is a column vector and "b" is a scalar.
    
    aL=cross(P2-P1,P3-P1); bL=P1'*aL;
    
    XV_test=vertices(facet_test_L',1);
    YV_test=vertices(facet_test_L',2);
    ZV_test=vertices(facet_test_L',3);
    
    % Warning:
    % This check is a little dangerous due to possible machine errors.
    % Must be improved in next releases.
    residual_test_vertices0=bL-aL'*[XV'; YV'; ZV'];
    
    residual_test_vertices=bL-aL'*[XV_test'; YV_test'; ZV_test'];
    [max_distance_test_vertex,index]=max(abs(residual_test_vertices));
    
    convex_check=residual_test_vertices(index);
    
    
    if convex_check < 0
        % warning('Modifying plane equation.');
        aL=-aL; bL=-bL;
    else
        if convex_check == 0
            warning('Troubles in determining polyhedron equation Ax<=b.');
        end
    end
        
    A=[A; aL']; b=[b; bL];
    
end