function bnd = obtainBnd(dataMatrix, dataType, image, ...
    sizeFraction, threshold, restoreFlag, resampleFlag)
% create boundary data structure from raw data
%
% Syntax: bnd = obtainBnd(dataMatrix)
%
% A structure is created that contains the measurement data in dataMatrix
%   (one data point each column) including some basic measures. The 
%   structure contains the fields cartX, cartY, polarAngle (radian), 
%   polarDistance that represent the input data. Additionally, the
%   following fields are returned:
%     polarBinSize: range of angle (radian) that is covered by each
%       distance in polarDistance
%     centerPointOriginal: coordinate origin as for the input points
%     centerPointMean: coordinate origin as centered to the mean
%     polarBoundaryArea: shape area estimated from circle sections
%     convexHullArea: area of convex hull
%     vrepHull: convex hull in V-rep
%     hrepHull: convex hull in H-rep
%     meanWidthHull: mean width of the convex hull
%     convexity: convexity ratio as [shape area]/[hull area]
%     largestAngleGap: largest gap in polarAngle vector
%     sphericity1: 1 - (maxPolar - minPolar)/meanPolar
%     sphericity2: 1 - stdPolar / meanPolar
%     sphericity3: pi*[diameter area eq. circle] / [circumfence of hull]
%     brightFraction: fraction of bright pixels (only with provided image)
%
% obtainBnd(.., type) sets the input data type. For type='cart', cartesian
%   coordinates are exprected. For type='polar', the input data is 
%   interpreted as polar coordinates with the distance in the first row and
%   the angle (radian) in the second row.
%
% obtainBnd(.., type, image, sizeFraction, threshold) Allows to provide the
%   corresponding image data to evaluate a brightness fraction of some 
%   area. The area is a circle with a radius of 0.5*sizeFraction*[minimum 
%   polar distance] that is centered in the bounding box. Bright pixels are
%   any pixels above the threshold. This measure is used to filter for air 
%   bubbles.
%
% See also: fitShapeBnd, viewBoundary

% The elk-library: convex geometry applied to crystallization modeling.
%   Copyright (C) 2014 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('dataType', 'var')
    dataType = 'cart';
end

%% input matrix check
% dimensions of data matrix
if size(dataMatrix, 1) == 2
    % everything OK
elseif size(dataMatrix, 2) == 2
    dataMatrix = dataMatrix';
else
    error('elk:boundaryCurve:wrongInput', ...
          'input data matrix must contains two rows, one data point in each column');
end

%% restore points
if exist('restoreFlag', 'var') && ~isempty(restoreFlag) && ...
        restoreFlag == 1 && strcmpi(dataType, 'cart')
    dataMatrix = restoreBoundary(dataMatrix);
end

%% evaluate brightness
if ~exist('sizeFraction', 'var'),sizeFraction = 0.5;end
if ~exist('threshold', 'var'),threshold=90;end
% if image is provided
if exist('image', 'var') && ~isempty(image) && strcmpi(dataType, 'cart')
    brightFraction = computeBrightFraction(image, ...
        dataMatrix(1,:), dataMatrix(2,:), ...
        sizeFraction, threshold);
end

%% build struct
% extract cartesian and polar data
if strcmpi(dataType, 'cart')
    % get data
    bnd.cartX = dataMatrix(1,:);
    bnd.cartY = dataMatrix(2,:);
    % center data
    [bnd.cartX, bnd.cartY bnd.centerPointOriginal] = ...
        centerCart(bnd.cartX, bnd.cartY);
    % create polar data
    [bnd.polarAngle bnd.polarDistance] = ...
        convertCartToPolar(bnd.cartX, bnd.cartY);
    % sort data
    [bnd.polarAngle bnd.polarDistance, permutationVector] = ...
        standardizePolar(bnd.polarAngle, bnd.polarDistance);
    bnd.cartX = bnd.cartX(permutationVector);
    bnd.cartY = bnd.cartY(permutationVector);
    
elseif strcmpi(dataType, 'polar')
    % get data
    bnd.polarAngle = dataMatrix(1,:);
    bnd.polarDistance = dataMatrix(2,:);
    % check origin in interior
    if min(bnd.polarDistance) < 0
        error('elk:boundary:wrongInput', ['The polar input data must have ' ...
            'positive distances only']);
    end
    % sort data
    [bnd.polarAngle bnd.polarDistance, ~] = ...
        standardizePolar(bnd.polarAngle, bnd.polarDistance);
    % create cart data
    [bnd.cartX bnd.cartY] = convertPolarToCart(bnd.polarAngle, ...
        bnd.polarDistance);
    bnd.centerPointOriginal = [0; 0];
else
    error('elk:boundaryCurve:wrongInput', ...
          'the specified data type must be ''cart'' or ''polar''');
end

if exist('resampleFlag', 'var') && resampleFlag > 0
    % resample data to an accuracy of resampleFlag
    angleVector = 0:resampleFlag:2*pi;
    distanceVector = 0*angleVector;
    for iAngle = 1:length(angleVector)
        thisFilter = (bnd.polarAngle >= (angleVector(iAngle) - resampleFlag/2)) | ...
                     (bnd.polarAngle >= (angleVector(iAngle) - resampleFlag/2 + 2*pi));
        thisFilter = thisFilter & ...
                     (bnd.polarAngle <= (angleVector(iAngle) + resampleFlag/2)) | ...
                     (bnd.polarAngle <= (angleVector(iAngle) + resampleFlag/2 - 2*pi));
        distanceVector(iAngle) = mean(bnd.polarDistance(thisFilter));
    end
    angleVector(isnan(distanceVector)) = [];
    distanceVector(isnan(distanceVector)) = [];
    bnd.polarAngle = angleVector;
    bnd.polarDistance = distanceVector;
    [bnd.cartX bnd.cartY] = convertPolarToCart(bnd.polarAngle, ...
        bnd.polarDistance);
end


%% add angle bin sizes
leftPolarAngle = 0.5*[bnd.polarAngle(end)-2*pi bnd.polarAngle(1:(end-1))] + ...
    0.5*bnd.polarAngle;
rightPolarAngle = 0.5*[bnd.polarAngle(2:end) bnd.polarAngle(1)+2*pi] + ...
    0.5*bnd.polarAngle;
bnd.polarBinSize = rightPolarAngle - leftPolarAngle;

%% correct bnd
bnd = moveBndToCenterOfMass(bnd);

%% add convexity information
bnd.polarBoundaryArea = sum(bnd.polarDistance.^2.*bnd.polarBinSize)/2;
[K bnd.convexHullArea] = convhull(bnd.cartX, bnd.cartY);
bnd.vrepHull.V = [bnd.cartX(K)' bnd.cartY(K)'];
bnd.hrepHull = convertVrepToHrep(bnd.vrepHull);
bnd.meanWidthHull = sum(computeEdgeLengthHrep(bnd.hrepHull, ...
    bnd.hrepHull.A))/pi;
bnd.convexity = bnd.polarBoundaryArea/bnd.convexHullArea;
bnd.largestAngleGap = max(diff([bnd.polarAngle bnd.polarAngle(1)+2*pi]));
bnd.sphericity = 1 - (max(bnd.polarDistance) - min(bnd.polarDistance)) ...
            / mean(bnd.polarDistance);
bnd.sphericity2 = 1 - std(bnd.polarDistance) ...
            / mean(bnd.polarDistance);
bnd.sphericity3 = 2*sqrt(bnd.convexHullArea/pi)/bnd.meanWidthHull;
if exist('brightFraction', 'var')
    bnd.brightFraction = brightFraction;
end

%% add circle fit
bnd.lambda = 0.5*max(bnd.polarDistance) + 0.5*min(bnd.polarDistance);

end
