function shape = obtainShapeFromFunctional(collectiveStruct, bnd, lambda, ...
    vertexIndex)
% extract shape information from collective structure
%
% This does also include roughness information.
%
% 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/>.

%% identify precalculated inner base shape
lambdaVector = [collectiveStruct(:).lambda];
baseIndex = find(lambdaVector <= lambda, 1, 'last');

%% shape information
baseShape = collectiveStruct(baseIndex);
shape.hrep = struct('A', baseShape.hrep.A, 'h', baseShape.hrep.h - ...
    (lambda - baseShape.lambda));
shape.lambda = lambda;
shape.edgeLengthVector = baseShape.edgeLengthFunctional(shape.hrep.h);
shape.vrep = baseShape.vrepFunctional(shape.hrep.h);


if ~exist('vertexIndex', 'var') || isempty(vertexIndex) || vertexIndex == 0
    %% overall distance
    shape.distanceVector = computeDistanceFromHrep2D(shape.hrep, shape.vrep, ...
        [bnd.cartX; bnd.cartY]);
else
    %% evaluating a specific edge
    % In this case, we only consider points that belong to a specific
    %   vertex of the initial shape. In doing so, we are able to fit this
    %   particular vertex.
    % The core is the proper selection of points for roughness calculation.
    %   Using too few point, especially when the amount increases later on,
    %   is a bad choice because roughness will be initially extremely
    %   small. A good choice is to identify the sample point right above
    %   the vertex, meaning: close to the angular bisecting line of the
    %   original vertex. Because the initial vertex vanishes with
    %   increasing lambda, this angular bisector is always shited to the
    %   currently supporting vertex. We then consider simply the 10% sample
    %   points around.
    % Each vertex itself features only a limited portion of the whole
    %   shape, given the opening angle of the adjacent edges. This opening
    %   angle is tracked, also for later evaluation (weighting of the 
    %   obtained lambda values).
    
    % The original edges creating the vertex are, with index and normal 
    %   vectors:
    originalEdgeFilter = collectiveStruct(1).incidenceMatrix(:, vertexIndex);
    originalEdgeNormalMatrix = collectiveStruct(1).hrep.A(originalEdgeFilter, :)';
    
    % we generate some angle information: angle of normals, opening angle
    %   and angle of the angle bisector (meanAngle)
    normalAngleVector = convertCartToPolar(originalEdgeNormalMatrix(1,:), ...
        originalEdgeNormalMatrix(2, :));
    if diff(normalAngleVector) < pi
        meanAngle = mean(normalAngleVector);
        shape.angleRatio = (max(normalAngleVector) ...
                          - min(normalAngleVector)) / 2 / pi;
    else
        meanAngle = mod(0.5*max(normalAngleVector) + ...
                        0.5*(min(normalAngleVector)+2*pi), 2*pi);
        shape.angleRatio = (min(normalAngleVector) ...
                          - max(normalAngleVector) + 2*pi) / 2 / pi;
    end
    
    % The supporting vertex for these edges is computed, by:
    scalarProductMatrix = shape.vrep.V * originalEdgeNormalMatrix;
    vertexIndex = find(...
        (abs(scalarProductMatrix(:, 1) - max(scalarProductMatrix(:, 1))) < 1e-10) & ...
        (abs(scalarProductMatrix(:, 2) - max(scalarProductMatrix(:, 2))) < 1e-10) );
    if numel(vertexIndex) ~= 1
        error('elk:boundary:numericalError', ['Supporting vertex for ' ...
            'initially adjacent edges is not found. This should not happen.']);
    end
    supportingVertex = shape.vrep.V(vertexIndex, :)';
    
    % the polar angle of the bisecting point on the shape is:
    bisectPoint = supportingVertex + lambda*computeNormalFromAngle(meanAngle);
    bisectAngle = convertCartToPolar(bisectPoint(1), bisectPoint(2));
    
    % the shape point sitting on the angular bisecting line is:
    closenessVector = abs(bnd.polarAngle - bisectAngle);
    bisectIndex = find(closenessVector == min(closenessVector), 1);
    
    % these are the indices of the neighbouring points:
    halfRange = round(length(bnd.polarAngle)/10/2);
    sampleIndexVector = mod(...
        ((bisectIndex-halfRange):(bisectIndex+halfRange))-1, ...
        length(bnd.polarAngle)) + 1;
    
    % get distance from base hrep
    shape.distanceVector = computeDistanceFromHrep2D(shape.hrep, shape.vrep, ...
        [bnd.cartX(sampleIndexVector); bnd.cartY(sampleIndexVector)]);
end

%% roughness information from distances
shape.roughness = max([0 (shape.distanceVector - lambda)]) ...
    + max([0 (lambda - shape.distanceVector)]);
shape.squareRoughness = sqrt(max([0 (shape.distanceVector - lambda)])^2 ...
    + max([0 (lambda - shape.distanceVector)])^2);
shape.meanSquare = mean((shape.distanceVector - lambda).^2);