function [coefs lblBranches] = polynomialCurveSetFit(seg, varargin)
%POLYNOMIALCURVESETFIT Fit a set of polynomial curves to a segmented image
%
%   COEFS = polynomialCurveSetFit(IMG);
%   COEFS = polynomialCurveSetFit(IMG, DEG);
%   Result is a cell array of matrices. Each matrix is DEG+1-by-2, and
%   contains coefficients of polynomial curve for each coordinate.
%   IMG is first binarised, then skeletonized. Each cure
%
%   [COEFS LBL] = polynomialCurveSetFit(...);
%   also returns an image of labels for the segmented curves. The max label
%   is the number of curves, and the length of COEFS.
%
%   Requires the toolboxes:
%   - Optimization
%   - Image Processing
%
%   Example
%   polynomialCurveSetFit
%
%   See also
%   polynomialCurves2d, polynomialCurveFit
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-03-21
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.

% default degree for curves
deg = 2;
if ~isempty(varargin)
    deg = varargin{1};
end


% ajoute un contour
seg([1 end], :) = 1;
seg(:, [1 end]) = 1;

% skeletise le segmentat
seg = bwmorph(seg, 'shrink', Inf);

% compute image of multiple points (intersections between curves)
imgNodes = imfilter(double(seg), ones([3 3])) .* seg > 3;

% compute coordinate of nodes, as c entroids of the multiple points
lblNodes = bwlabel(imgNodes, 4);
struct   = regionprops(lblNodes, 'Centroid');
nodes = zeros(length(struct), 2);
for i=1:length(struct)
    nodes(i, 1:2) = struct(i).Centroid;
end

% enleve les bords de l'image
seg([1 end], :) = 0;
seg(:, [1 end]) = 0;

% Isoles les branches
imgBranches = seg & ~imgNodes;
lblBranches = bwlabel(imgBranches, 8);

% % donne une couleur a chaque branche, et affiche
% map = colorcube(max(lblBranches(:))+1);
% rgbBranches = label2rgb(lblBranches, map, 'w', 'shuffle');
% imshow(rgbBranches);

% number of curves
nBranches = max(lblBranches(:));

% allocate memory
coefs = cell(nBranches, 1);


% For each curve, find interpolated polynomial curve
for i = 1:nBranches
    disp(i);
    
    % extract points corresponding to current curve
    imgBranch = lblBranches == i;
    points = chainPixels(imgBranch);
    
    % check number of points is sufficient
    if size(points, 1) < max(deg+1-2, 2)
        % find labels of nodes
        inds = unique(lblNodes(imdilate(imgBranch, ones(3,3))));
        inds = inds(inds ~= 0);
        
        if length(inds) < 2
            disp(['Could not find extremities of branch number ' num2str(i)]);
            continue;
        end
        
        % consider extremity nodes
        node0 = nodes(inds(1), :);
        node1 = nodes(inds(2), :);
        
        % use only a linear approximation
        xc = zeros(1, deg+1);
        yc = zeros(1, deg+1);
        xc(1) = node0(1);
        yc(1) = node0(2);
        xc(2) = node1(1)-node0(1);
        yc(2) = node1(2)-node0(2);
        
        % assigne au tableau de courbes
        coefs{i} = [xc;yc];
        
        % next branch
        continue;
    end

    % find nodes closest to first and last points of the current curve
    [dist, ind0] = minDistancePoints(points(1, :), nodes); %#ok<*ASGLU>
    [dist, ind1] = minDistancePoints(points(end, :), nodes);
    
    % add nodes to the curve.
    points = [nodes(ind0,:); points; nodes(ind1,:)]; %#ok<AGROW>
    
    % parametrization of the polyline
    t = parametrize(points);
    t = t / max(t);
    
    % fit a polynomial curve to the set of points
    [xc yc] = polynomialCurveFit(...
        t, points, deg, ...
        0, {points(1,1), points(1,2)},...
        1, {points(end,1), points(end,2)});
    
    % stores result
    coefs{i} = [xc;yc];
end




function points = chainPixels(img, varargin)
%CHAINPIXELS return the list of points which constitute a curve on image
%   output = chainPixels(input)
%
%   Example
%   chainPixels
%
%   See also
%
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2007-03-21
% Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.


conn = 8;
if ~isempty(varargin)
    conn = varargin{1};
end

% matrice de voisinage
if conn == 4
    f = [0 1 0;1 1 1;0 1 0];
elseif conn == 8
    f = ones([3 3]);
end

% find extremity points
nb = imfilter(double(img), f) .* img;
imgEnding = nb == 2 | nb == 1;
[yi xi] = find(imgEnding);

% extract coordinates of points
[y x] = find(img);

% index of first point
if isempty(xi)
    % take arbitrary point
    ind = 1;
else
    ind = find(x==xi(1) & y==yi(1));
end

% allocate memory
points  = zeros(length(x), 2);

if conn == 8
    for i = 1:size(points, 1)
        % avoid multiple neighbors (can happen in loops)
        ind = ind(1);
        
        % add current point to chained curve
        points(i,:) = [x(ind) y(ind)];

        % remove processed coordinate
        x(ind) = [];
        y(ind) = [];

        % find next candidate
        ind = find(abs(x-points(i,1))<=1 & abs(y-points(i,2))<=1);
    end
    
else
    for i = 1:size(points, 1)
        % avoid multiple neighbors (can happen in loops)
        ind = ind(1);
        
        % add current point to chained curve
        points(i,:) = [x(ind) y(ind)];

        % remove processed coordinate
        x(ind) = [];
        y(ind) = [];

        % find next candidate
        ind = find(abs(x-points(i,1)) + abs(y-points(i,2)) <=1 );
    end
end    
