function esDec = computeExtendedEmbedding(esDec, zerotol)
% compute extended embedding from H- and original representation
%
% * H-representation given by esDec.A
% * original representation is given by esDec.dataOriginal (as
%     decomposition)
%
% THIS IS NO USER FUNCTION

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

%% Lift structuring elements to the extended H-representation
% therefore, we cycle through all structuring elements and build a matrix
%   that gives the SEs in the extended H-representation.
nExtraFacet = size(esDec.originalData.newFacetNormalMatrix, 1);
seFacetDistanceMatrix = [...
    esDec.originalData.mappingReducedToFull * esDec.originalData.se'; ...
    zeros(nExtraFacet, esDec.originalData.nS)];
for iSe = 1:esDec.originalData.nS
    % generate se in hrep
    thisHrep = struct('A', esDec.originalData.A, 'h', ...
                       esDec.originalData.mappingReducedToFull * ...
                       esDec.originalData.se(iSe, :)');
    % compute support values for additional faces
    seFacetDistanceMatrix((esDec.originalData.nH+1):end, iSe) = ...
        computeSupportValueHrep(thisHrep, ...
                                esDec.originalData.newFacetNormalMatrix);
end

%% Find new constrained H-representation (extended embedding)
% Essentially, we search a linear subspace in the new h-space. With SVD
%   this is easily calculated giving any orthonormal basis. Since it is not
%   guaranteed to obtain a proper representation, this mapping is
%   sufficient.
[u s v] = svd(seFacetDistanceMatrix);
esDec.nC = sum(diag(s) > zerotol);
esDec.mappingReducedToFull = u(:, 1:esDec.nC);
esDec.mappingFullToReduced = u(:, 1:esDec.nC)';


%% Adding some more (already known information)

% copy unchanged information
esDec.nDim = esDec.originalData.nDim;
esDec.nS = esDec.originalData.nS;

% merge facet information
esDec.constraintMatrix = eye(esDec.nH) - ...
    esDec.mappingReducedToFull * esDec.mappingFullToReduced;

% adopt (real-valued) structuring elements for additional facets. note,
%   that the real-scaling (mean width) is still valid because we've only
%   added facets that do not exist on those shapes.
esDec.se = (esDec.mappingFullToReduced * seFacetDistanceMatrix)';

% scale to close information
esDec.centerRay = mean(esDec.se, 1)';
esDec.seScaledToClose = scaleToClose(esDec.se, esDec);

% validity cone - note that the valdity cone is not the true cone(!) it is
%   chosen so that the resulting representation survives validateDec
coneVrep.V = closeCone(esDec.seScaledToClose);
coneHrep = convertVrepToHrep(coneVrep);
esDec.validityCone = openCone(coneHrep);
esDec.validityCone.seIndexVector = 1:esDec.nS;
esDec.validityCone.volume = computeVolumeVrep(coneVrep);
esDec.confinementCone = esDec.validityCone;

% partition information - see validity cone: volumeScaled is wrong and the
%   partitioning is NOT complete
esDec = appendUnifiedPartition(esDec);
% this partition should be valid
esDec = appendUnifiedPartition(esDec, 1:esDec.nS, zerotol, 1, 1);
esDec.nPU = 1;
esDec.nPS = 1;

% the output representation is not necessarily proper; in fact, we dont
%   even know the validity and confinement cones; the same applies for
%   singular/ordinary
esDec.isProper = 0;
esDec.isOrdinary = 0;
esDec.isComplete = 0;

% cleanup (so far)
esDec = cleanupDecStructure(esDec);