function newFacetNormalMatrix = computeAdditionalFacetMatrix(dec, zerotol)
% search for additional facets in pairs of SEs
%
% THIS IS NO USER FUNCTION

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2014 Alexander Reinhold
%
% This program is free software: you can redistribute it and/or modify it 
%   under the terms of the GNU General Public License as published by the 
%   Free Software Foundation, either version 3 of the License, or (at your 
%   option) any later version.
%
% This program is distributed in the hope that it will be useful, but 
%   WITHOUT ANY WARRANTY; without even the implied warranty of 
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU 
%   General Public License for more details.
%
% You should have received a copy of the GNU General Public License along 
%   with this program.  If not, see <http://www.gnu.org/licenses/>.

% the following loop goes through all pairwise sums of structuring elements
%   that are not adjacent (are not added to represent crystals / do not
%   have a unified partition in common).
% These pairs are expected to deliver a new set of facets. If no extra
%   facets are created, an error is thrown. The facets are listed in
%   newFacetNormalMatrix and a logical matrix tracks in which summations
%   which facets were created. These sets are expected to be distinct,
%   which is also explicitely checked (only pairwise).
% From this facet creation matrix the new facet-groups are identified.
[rowIndexVector, colIndexVector] = find(...
    (dec.seIncidenceMatrix + tril(ones(dec.nS) - eye(dec.nS))) == 0);
newFacetNormalMatrix = [];
% facetCreationMatrix = [];

for iPair = 1:length(rowIndexVector)
    % The two summands
    hrepSummand1 = dec.se(rowIndexVector(iPair), :) * dec.mappingReducedToFull';
    hrepSummand1 = struct('A', dec.A, 'h', hrepSummand1');
    hrepSummand2 = dec.se(colIndexVector(iPair), :) * dec.mappingReducedToFull';
    hrepSummand2 = struct('A', dec.A, 'h', hrepSummand2');
    
    % Minkowski addition of pair
    hrepSum = addMinkowskiHrep(hrepSummand1, hrepSummand2);
    
    % it might happen that the resulting shape is a flat object. In this
    %   case we have (exactly) two opposing facet normals with distance 0
    %   and the other facet normals have absolutely no meaning. The facet
    %   distance might also be set off by translation.
    if computeVolumeHrep(hrepSum, zerotol) < zerotol
        vrepSum = convertHrepToVrep(hrepSum);
        % create a matrix with each row corresponding to a vertex and each
        %   column to a facet.
        facetFilter = (hrepSum.A * vrepSum.V')';
        % The (two) coumns where all values are equal correspond to the two
        %   facet normals we search for. The first is sufficient.
        for iFacet = 1:size(facetFilter, 2)
            if all(abs(facetFilter(:,iFacet) - facetFilter(1,iFacet)) < zerotol)
                break;
            end
        end
        if iFacet == size(facetFilter, 2)
            error('elk:decomposition:numericalError', ...
                  'Could not find the facets of a flat object. This should not happen');
        end
        % obtain facet normals
        thisFacetNormalMatrix = [hrepSum.A(iFacet, :); ...
                              -1*hrepSum.A(iFacet, :)];
        
        % Filter original facets
        thisFacetNormalMatrix = eliminateRows(thisFacetNormalMatrix, ...
                                              dec.A, zerotol);
    else
        % Filter original facets
        thisFacetNormalMatrix = eliminateRows(hrepSum.A, dec.A, zerotol);
    end
 
    % throw error, if there are no additional facets
    if isempty(thisFacetNormalMatrix)
        error('elk:decomposition:numericalError', ...
            ['Structuring elements ' num2str(rowIndexVector(iPair)) ' and ' ...
             num2str(colIndexVector(iPair)) ' do not share a unified partition ' ...
             ' but also do not create additional facets by Minkowski addition.']);
    end
    
    % Filter already known facets
    [thisFacetNormalMatrix, thisCreatedFacetIndexVector] = eliminateRows(...
        thisFacetNormalMatrix, newFacetNormalMatrix, zerotol);
    
    % Append new facets
    thisCreatedFacetIndexVector = [thisCreatedFacetIndexVector, ...
          (1:size(thisFacetNormalMatrix, 1)) + size(newFacetNormalMatrix, 1)];
    newFacetNormalMatrix = [newFacetNormalMatrix; thisFacetNormalMatrix];
    
    % Enter information in facet creation matrix
%     facetCreationMatrix(thisCreatedFacetIndexVector, iPair) = 1;
end
