function result = fitShapeBnd(bnd, varargin)
% fit a rounded shape to boundary data
% 
% Syntax: result = fitShapeBnd(bnd)
%
% This function takes a boundary curve (generated by obtainBnd) and fits a
%   rounded shape. The resulting structure contains the final result but
%   also the original boundary, characteristic values and intermediate
%   steps of the fitting process. 
%
% The returned structure contains (despite others) the fields:
%     lambda - the estimated additive sphericity in units of the input data
%     circumfence - circumfence in units of the input data
%     meanWidth - mean width (=mean Feret diameter) in units of the input 
%       data
%     roughness - the maximal deviation of input points from the estimated
%       shape
%     shapeStruct - the estimated shape as a structure array for each edge
%       or vertex, respectively. It has the following fields:
%     >.edgeAngle - the angle of each edge in radian
%     >.edgeDistance - the distance of the edge from the origin
%     >.edgeCart - a matrix containing the start point (first column) and
%       end point (second column) of the edge
%     >.vertexOrigin - the center of the vertex circle
%     >.vertexRadius - the radius of the vertex radius
%     >.vertexStartAngle - start angle for vertex circle
%     >.vertexEndAngle - end angle for vertex circle 
%
% Optional parameters in the form fitShapeBnd(bnd, 'option', value) are the
%   following:
%     'tolAngle' - accuracy for edge angle in radian. Default: 2/180*pi
%     'minDiffAngle' - minimum difference for the angle of two edges. 
%       Default:, 15/180*pi
%     'tolDistance' - accuracy for edge distance in pixel. Default: 1
%     'tolLambda' - accuracy for the fitted radius. Default: 1
%     'lambdaLocal' - switch for local (1) or global (0) approaches. 
%       Default: 0
%     'lambdaError' - selected objective method with 'meanSq' for mean 
%       squares and 'sqRough' for  // meanSq
%   The following activate additional filters based on values of the 
%   boundary data (see obtianBnd). Particles are not processed if the 
%   filter criteria is violated.
%     'minConvexity' - minimum convexity ratio - Default: 0.9
%     'minDiameter' - minimum diameter - Default: 100
%     'maxSphericity' - maximum sphericity - Default: 0 (filter off)
%     'maxSphericity2' - maximum sphericity2 - Default: 0 (filter off)
%     'minBrightFraction' - minimum bright fraction (bubble detection).
%       Default: 0 (filter off)
%     'maxAngleGap' - maximum gap in polar angles of input data
%       (convexity). Default: 0 (filter off)
%
% See also: obtainBnd, viewBoundary, viewBoundaryShape

% 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/>.

%! DEV: The tolerances were tweaked and denoted by new values:
%        new tolerances - angle: 2/180*pi, tolDistance 1, tolLambda 1
%   but the server images went missing so that it could not be published
%   with these values. The old varlues are:
%        old tolerances - angle 0.05; distance 2, tolLambda 2
%/
optionStruct = mapOptionStruct(varargin, ...
    'tolAngle', 2/180*pi, ...
    'minDiffAngle', 15/180*pi, ...
    'tolDistance', 1, ...
    'tolLambda', 1, ...
    'lambdaLocal', 0, ...
    'lambdaError', 'meanSq', ... % sqHeight // meanSq
    'minConvexity', 0.9, ...
    'minDiameter', 100, ...
    'maxSphericity', 0, ...
    'maxSphericity2', 0, ...
    'minBrightFraction', 0, ...
    'maxAngleGap', 0, ...
    'debugHough', 0, ...
    'debugEdgeResult', 0, ...
    'debugLambda', 0);

%% ignore particles
result = handleIgnoreFilter(bnd, optionStruct);
if ~isempty(result)
    return
end

%% identify edges
if optionStruct.debugHough
    [hrep debug.edge.hough] = identifyEdges(bnd, optionStruct);
else
    hrep = identifyEdges(bnd, optionStruct);
end
debug.edge.hrep = hrep;

%% debug code
if optionStruct.debugEdgeResult
    figure,hold on;
    viewBoundary(bnd);
    viewPolytope(hrep, 'shade', 0, 'edgeColor', [1 0 0]);
    title('complete set of edges, lambda = 0')
end

%% Precompute edge reduction
% the following struct contains all information of subsequently reduced
%   edges for P(h-lambda).
collectiveStruct = computeCollectiveShapeStruct(hrep, ...
    1e-5*optionStruct.tolDistance);
% the next variable denotes where no inner polytope exists anymore:
seperatingLambda = collectiveStruct(end).lambda;
debug.lambda.collectiveStruct = collectiveStruct;

%% Identify lambda
% The shape model is P(h) + lam*B and for P(h) the edge angles are already 
%   determined by the previous step. The previous step also identified edge
%   distances so that we might rewrite the model as (P(hI-lam) + lam*B).
% The roughness is defined by two limiting shapes. A lower bound shape:
%     (P(hI-lam) + lam*B) - gL*B
%   uses the maximum gL so that no sample point is contained in that shape.
%   And an upper bound shape:
%     (P(hI-lam) + lam*B) + gU*B
%   uses a minimum gU so that all points are contained in that shape. Note
%   that these decisions can be made based on the (signed) distance of
%   all points from the polytope P(hI-lam). Signed means using negative
%   values for inner points and positive values for outer points.
% Rouhness is, however not an ideal objective function. Usually, the
%   roughness does not change if lam is chosen in the interval between 
%   (lam-gL) and (lam+gU) so that both extreme cases are possible. To avoid
%   this, we use a squarred roughness:
%     squareRoughness = sqrt(gU^2 + gL^2)
%   so that lambda should lie in the middle of (lam-gL) and (lam+gU).

%% all vertices at once
% This old code fragment fits the whole shape as one. It was found that
%   this method often fits to the smallest lambda that is featured by an
%   vertex. As soon as lambda exceeds this value, gU rises. Even if a lot
%   more vertices feature a larger lambda value, this will not be taken
%   into account.

% getObjectFunction = @(lam)(obtainShapeFromFunctional(...
%     collectiveStruct, bnd, lam));
% objectiveFunction = @(shape)(shape.squareRoughness);
% 
% shape = lineSearchGolden(getObjectFunction, objectiveFunction, ...
%         0, seperatingLambda, optionStruct.tolLambda);
    
if optionStruct.lambdaLocal
    %% seperate vertices
    % in this approach, we consider each vertex of the initial
    %   H-representation.
    nVertex = size(collectiveStruct(1).vrep.V, 1);
    lambdaVector = nan(1, nVertex);
    angleRatioVector = nan(1,nVertex);
    for iVertex = 1:nVertex
        getObjectFunction = @(lam)(obtainShapeFromFunctional(...
            collectiveStruct, bnd, lam, iVertex));
        if strcmpi(optionStruct.lambdaError, 'sqHeight') || ...
           strcmpi(optionStruct.lambdaError, 'sqRough')
            % Well, square heigt was named square roughness in the early
            %   stages.
            objectiveFunction = @(shape)(shape.squareRoughness);
        elseif strcmpi(optionStruct.lambdaError, 'meanSq')
            objectiveFunction = @(shape)(shape.meanSquare);
        end
        
        shape = lineSearchGolden(getObjectFunction, objectiveFunction, ...
            0, seperatingLambda, optionStruct.tolLambda);
        
        lambdaVector(iVertex) = shape.lambda;
        %     angleRatioVector(iVertex) = shape.angleRatio;
        angleRatioVector(iVertex) = 1/nVertex;
        debug.lambda.shape(iVertex) = shape;
        
        %! Debug:
        %     figure,hold on
        %     h = viewBoundaryShape(collectiveStruct(1));set(h, 'color', [0 1 0]);
        %     viewBoundary(bnd),viewBoundaryShape(shape)
        %     plot(collectiveStruct(1).vrep.V(iVertex, 1), ...
        %          collectiveStruct(1).vrep.V(iVertex, 2), 'm.', 'markerSize', 20);
        %/
    end
    % we use a mean from the lambdas, weighted by the corresponding opening
    %   angle:
    thisLambda = angleRatioVector * lambdaVector';
    shape = obtainShapeFromFunctional(collectiveStruct, bnd, thisLambda);
    % safe local fitting information
    debug.lambdaVector = lambdaVector;
    debug.angleRatioVector = angleRatioVector;
else
    %% global lambda approach
    getObjectFunction = @(lam)(obtainShapeFromFunctional(...
        collectiveStruct, bnd, lam, 0));
    if strcmpi(optionStruct.lambdaError, 'sqHeight') || ...
       strcmpi(optionStruct.lambdaError, 'sqRough')
        objectiveFunction = @(shape)(shape.squareRoughness);
    elseif strcmpi(optionStruct.lambdaError, 'meanSq')
        objectiveFunction = @(shape)(shape.meanSquare);
    end
    
    shape = lineSearchGolden(getObjectFunction, objectiveFunction, ...
        0, seperatingLambda, optionStruct.tolLambda);
end

result = mergeStruct(shape, bnd);
result.debug = debug;
result.flagIgnore = 0;
result.circumfence = sum(result.edgeLengthVector) + 2*pi*result.lambda;
result.meanWidth = result.circumfence / pi;
result.shapeArea = computeVolumeVrep(result.vrep) + ...
    result.lambda*result.meanWidth + pi*result.lambda^2;
