function tensor = computeMixedVolumeTensor(hMatrix, handleH, dim, ...
    incidenceMatrix, verbose)
% calculate the coefficients for measure calculation in S-representation
%
% IMPORTANT: this function uses the sum of hC-vectors instead of the
%   Minkoski sum. This is valid as long as the tensor of mixed volume is
%   applied to unified partitions. It CANNOT be applied to arbitrary
%   combinations of structuring elements.
%
% THIS FUNCTION IS NO USER FUNCTION.

% 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/>.
if ~exist('verbose', 'var') || isempty(verbose)
    verbose = 0;
end

if verbose
    disp('Computing mixed volumes..');
    fprintf('......');
end
switch dim
    case 1
        tensor = zeros([1 size(hMatrix, 1)]);
        for iSe = 1:size(hMatrix, 1)
            tensor(iSe) = handleH(hMatrix(iSe, :)');
        end
        
    case 2
        tensor = nan([size(hMatrix, 1) size(hMatrix, 1)]);
        for iSe = 1:size(hMatrix, 1)
            tensor(iSe,iSe) = handleH(hMatrix(iSe, :)');
        end
        % Measure for all combinations of two polytopes
        % V(lA*A + lB*B) = lA^2*V(A) + lB^2*V(B) + 2*lA*lB V(A, B)
        %
        % V(A,B) = (V_AB - V(A) - V(B)) / 2
        % but as V(A,B) is applied it would be multiplied by 2, hence, we
        % omit this.
        for iSe = 1:size(hMatrix, 1)
            for jSe = (iSe+1):size(hMatrix, 1)
                if incidenceMatrix(iSe, jSe) == 0
                    continue
                end
                tensor(iSe,jSe) = (handleH(...
                                     hMatrix(iSe, :)' + hMatrix(jSe, :)')...
                                   - tensor(iSe,iSe) ...
                                   - tensor(jSe,jSe));
            end
        end
    case 3
        % initialise 3D matrix
        tensor = nan([size(hMatrix, 1), ...
                        size(hMatrix, 1), ...
                        size(hMatrix, 1)]);
        
        % Measure for all structuring elements, which can be calculated
        % directly.
        for iSe = 1:size(hMatrix, 1)
            if verbose
                fprintf('\b\b\b\b\b\b%5.1f%%', ...
                    (iSe/size(hMatrix, 1))*100);
            end
            tensor(iSe, iSe, iSe) = handleH(hMatrix(iSe, :)');
        end
        
        % Measure for all combinations of two polytopes
        % V(lA A + lB B) = lA^3 * V(A) + 3*lA^2*lB V(A, A, B) + ...
        %                  lB^3 * V(B) + 3*lA*lB^2 V(A, B, B)
        %
        % V(A + B)  = V(A) +   V(B) + 3*V(A, A, B) +  3*V(A, B, B)
        % V(A + 2B) = V(A) + 8*V(B) + 6*V(A, A, B) + 12*V(A, B, B)
        %
        % V(A + 2B) - 2*V(A + B) = -V(A) + 6*V(B) + 6*V(A, B, B)
        % V(A, B, B) = (V(A + 2B) - 2*V(A + B) + V(A) - 6 V(B)) / 6
        %
        % V(A, A, B) = (V(A + B) - V(A) - V(B) - 3*V(A, B, B)) / 3
        for iSe = 1:size(hMatrix, 1)
            if verbose
                fprintf('\b\b\b\b\b\b%5.1f%%', ...
                       (iSe/size(hMatrix, 1))*100);
            end
            for jSe = (iSe+1):size(hMatrix, 1)
                if incidenceMatrix(iSe, jSe) == 0
                    continue
                end
                V_AB = handleH(hMatrix(iSe,:)' + hMatrix(jSe,:)');
                V_A2B = handleH(hMatrix(iSe,:)' + 2*hMatrix(jSe,:)');
                tensor(iSe,jSe,jSe) = (V_A2B - 2*V_AB ...
                                 + tensor(iSe,iSe,iSe) ...
                                 -6*tensor(jSe,jSe,jSe) ...
                                )/2;
                tensor(iSe,iSe,jSe) = (V_AB ...
                                 - tensor(iSe,iSe,iSe) ...
                                 - tensor(jSe,jSe,jSe) ...
                                 - tensor(iSe,jSe,jSe) ...
                                );
            end
        end
        
        % Measure for all combinations of three polytopes
        % V(A+B+C) = V(A) + V(B) + V(C) + ...
        %            3*V(A,B,B) + 3*V(A,C,C) + 3*V(B,C,C) +
        %            6*V(A,B,C)
        % V(A,B,C) = (V(A+B+C) - V(A) - V(B) - V(C) ...
        %             -3*V(A,B,B) -3*V(A,C,C) -3*V(B,C,C)) / 6
        for iSe = 1:size(hMatrix, 1)
            if verbose
                fprintf('\b\b\b\b\b\b%5.1f%%', ...
                    (iSe/size(hMatrix, 1))*100);
            end
            for jSe = (iSe+1):size(hMatrix, 1)
                if incidenceMatrix(iSe, jSe) == 0
                    continue
                end
                for c = (jSe+1):size(hMatrix, 1)
                    if incidenceMatrix(iSe, c) == 0 || ...
                       incidenceMatrix(jSe, c) == 0
                        continue
                    end
                    V_ABC = handleH(hMatrix(iSe,:)' + hMatrix(jSe,:)' + hMatrix(c,:)');
                    tensor(iSe,jSe,c) = (V_ABC ...
                                     -tensor(iSe,iSe,iSe) ...
                                     -tensor(jSe,jSe,jSe) ...
                                     -tensor(c,c,c) ...
                                     -tensor(iSe,jSe,jSe) ...
                                     -tensor(iSe,c,c) ...
                                     -tensor(jSe,c,c) ...
                                     -tensor(iSe,iSe,jSe) ...
                                     -tensor(iSe,iSe,c) ...
                                     -tensor(jSe,jSe,c) ...
                                    );
                end
            end
        end
end
if verbose, fprintf('\b\b\b\b\b\b');end
end