function paintFacet(patchHandleVector, facetNormalMatrix, colorMatrix, varargin)
% paints facets of a polytope
%
% Syntax: paintFacet(patchHandleVector, facetNormalMatrix, colorMatrix,)
%
% A typical call could be:
%     paintFacet(get(gca, 'Children), cdef.A, [1 0 0])
%   to paint all facets of the crystal definition in red for a polytope
%   that is just viewed in a new figure.
% The rows of colorMatrix define the colors of facets in the order
% specified in facetNormalMatrix
%
% See also: viewPolytope

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2012 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/>.


%% make color matrix if only a single color is specified
if numel(colorMatrix) == 1
    colorMatrix(size(facetNormalMatrix, 1), :) = colorMatrix;
elseif (size(colorMatrix, 1) == 3 && size(colorMatrix, 2) == 1) || ...
       (size(colorMatrix, 1) == 1 && size(colorMatrix, 2) == 3)
    tmp = colorMatrix;
    colorMatrix(1:size(facetNormalMatrix, 1), 1) = tmp(1);
    colorMatrix(1:size(facetNormalMatrix, 1), 2) = tmp(2);
    colorMatrix(1:size(facetNormalMatrix, 1), 3) = tmp(3);
elseif size(colorMatrix, 1) == size(facetNormalMatrix, 1) && ...
       size(colorMatrix, 2) == 3
    %everything OK
else
    error('elk:polytope:wrongInput', ...
        'the color matrix is wrong, try "help paintFacet"');
end

% ensure only patches are in the handles
patchHandleVector(~strcmp(get(patchHandleVector, 'type'), 'patch')) = [];

patchX = [];
patchY = [];
patchZ = [];
for iPatch = 1:length(patchHandleVector)
    patchX = [patchX; get(patchHandleVector(iPatch), 'XData')];
    patchY = [patchY; get(patchHandleVector(iPatch), 'YData')];
    patchZ = [patchZ; get(patchHandleVector(iPatch), 'ZData')];
end

origin = [(max(patchX) + min(patchX))/2; ...
          (max(patchY) + min(patchY))/2; ...
          (max(patchZ) + min(patchZ))/2];
      
%% cycle through facets and move the patches coordinates so that 
% origin is in the middle of the polytope
for iPatch = 1:length(patchHandleVector)     
    patchX = get(patchHandleVector(iPatch), 'XData'); 
    patchX = patchX - origin(1) * ones(size(patchX));
    patchY = get(patchHandleVector(iPatch), 'YData');
    patchY = patchY - origin(2) * ones(size(patchY));
    patchZ = get(patchHandleVector(iPatch), 'ZData');
    patchZ = patchZ - origin(3) * ones(size(patchZ));
    
    % distances of each point from each facet
    idx = facetNormalMatrix*[patchX';patchY';patchZ'];
    
    % This checks if all distances for one facet have the same positive
    % value, meaning they are on that facet.
    % The value should be positive as the origin is in the middle of
    % the polytope. Therefore, if the values are the same and positive, 
    % idx will be zero everywhere for that facet, if all values are same
    % and negative, they will have twice the original value
    idx = idx - abs(idx(:,1)) * ones([1, size(idx,2)]);
    
    % find the rows with all values equal to zero
    idx = abs(idx) < 1e-6;
    idx = prod(double(idx), 2);
    
    % take first valid facet
    idx = find(idx, 1, 'first');
    
    % color patch
    if isempty(idx),continue;end
    set(patchHandleVector(iPatch), 'FaceColor', colorMatrix(idx, :), varargin{:});
end