function varargout = identifyEdges(bnd, optionStruct)
% identify edges of a boundary
%
% 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/>.

% The main idea comes the Hough transform, an algorithm to identify lines
%   in a binary image. The idea is that every line has two parameters:
%   orientation and distance from origin. The Hough transfers an image with
%   pixely in x/y direction into an image with pixels according to
%   angles/distances of the edges. Hence, it discretizes the parameter
%   space. In principle, it takes each pixel and direction, determines the
%   required distance of the line and increases the corresponding bin of
%   the Hough-image.
% In this approach, we will do the same for every sample point. Bright
%   points (maxima) in the Hough-image will appear for a line that fits the
%   boundary curve well, because lot's of points contribute to these
%   specific parameters.
% To allow a reasonable detection of proper parameters, we also include the
%   understanding of roughnes and minimal edge lengths. If we assume a
%   certain allowed roughness, then much more points can be assigned to a
%   point in Hough-space. This increases the imporance of the maxima
%   compared to the lower values and allows to use the minimal edge length
%   in relation to the total number of pixels in the image.

%% Hough grid
% basic grid spacing
dMin = min(bnd.polarDistance)*0.95;
dMax = max(bnd.polarDistance)*1.05;
deltaAngle = optionStruct.tolAngle;
deltaDistance = optionStruct.tolDistance;
% grid points
houghAngleVector = 0:deltaAngle:(2*pi);
houghDistanceVector = dMin:deltaDistance:dMax;
% grid point boundaries for distance
distanceLowerBoundVector = houghDistanceVector - deltaDistance/2;
distanceUpperBoundVector = houghDistanceVector + deltaDistance/2;
% obtain grid (for later plotting)
[dGrid phiGrid] = meshgrid(houghDistanceVector, houghAngleVector);
% initialize hough image matrix
houghMatrix = dGrid*0;

%% evaluate Hough values
% normals according to angle grid points
normalMatrix = computeNormalFromAngle(houghAngleVector);
% determine distances according to angles
distanceMatrix = normalMatrix' * [bnd.cartX; bnd.cartY];
% fill bins by cycling through available distances
for iDist = 1:length(houghDistanceVector)
    thisPlus = sum(...
        (distanceMatrix < distanceUpperBoundVector(iDist)) & ...
        (distanceMatrix >= distanceLowerBoundVector(iDist)), 2);
    houghMatrix(:, iDist) = houghMatrix(:, iDist) + thisPlus;
end
% we standardize the values to the ratio of accumulated points by total
%   points
houghMatrix = houghMatrix/length(bnd.cartX);

if optionStruct.debugHough > 1
    % hough plot
    figure
    contourf(phiGrid/2/pi*360, dGrid, ...
        houghMatrix * 100)
    colorbar
    title('accumulated counts in percent of total pixels')
    
    % edge angle plot
    figure, hold on
    plot(houghAngleVector/2/pi*360, max(houghMatrix*100, [], 2))
    title('maximum hough value in percent of total pixels')
end

%% find maxima with angle-max from hough - sliding maxima
[angleRatingVector distanceIndexVector] = max(houghMatrix, [], 2);
edgeAngleIndexVector = identifyMaximaSlidingWindow(angleRatingVector, ...
    round(optionStruct.minDiffAngle / deltaAngle));
% remove obviously short edges
% edgeAngleIndexVector(...
%     angleRatingVector(edgeAngleIndexVector) < ...
%     0.3*max(angleRatingVector)) = [];

% assign H-representation
% houghAngleVector(edgeAngleIndexVector)/pi
hrep.A = computeNormalFromAngle(houghAngleVector(edgeAngleIndexVector))';
hrep.h = houghDistanceVector(distanceIndexVector(edgeAngleIndexVector))';
% tmp = distanceIndexVector(edgeAngleIndexVector)
% houghDistanceVector(tmp)

% %% adjust facet distance to mean square
% for iFacet = 1:size(hrep.A, 1)
%     distanceVector = hrep.A(iFacet, :) * [bnd.cartX; bnd.cartY] - ...
%                      hrep.h(iFacet);
%     filter = abs(distanceVector) < 0.5*optionStruct.tolDistance;
%     hrep.h = hrep.h + mean(distanceVector(filter));
% end

varargout{1} = hrep;
if nargout > 1
    hough.angleGrid = phiGrid;
    hough.distanceGrid = dGrid;
    hough.angleVector = houghAngleVector;
    hough.distanceVector = houghDistanceVector;
    hough.image = houghMatrix;
    varargout{2} = hough;
end