Commit [01a21d] default Maximize Restore History

geometry: adding polynomial curve functions

jpicarbajal jpicarbajal 2013-03-03

removed devel/polynomialCurves2d/polynomialCurveFit.m
changed INDEX
changed PKG_ADD
changed NEWS
changed PKG_DEL
copied devel/polygons2d/drawPolyline.m -> inst/polygons2d/drawPolyline.m
copied devel/polynomialCurves2d/polynomialDerivate.m -> inst/polynomialCurves2d/drawPolynomialCurve.m
copied devel/polynomialCurves2d/polynomialCurvePoint.m -> inst/polynomialCurves2d/polynomialCurvePoint.m
copied devel/polynomialCurves2d/polynomialCurveSetFit.m -> inst/polynomialCurves2d/polynomialCurveSetFit.m
copied devel/polynomialCurves2d/drawPolynomialCurve.m -> inst/polynomialCurves2d/polynomialCurveDerivative.m
copied devel/polynomialCurves2d/polynomialCurveDerivative.m -> inst/polynomialCurves2d/polynomialCurveCentroid.m
copied devel/polynomialCurves2d/polynomialCurveCentroid.m -> inst/polynomialCurves2d/polynomialCurveFit.m
INDEX Diff Switch to side-by-side view
Loading...
PKG_ADD Diff Switch to side-by-side view
Loading...
NEWS Diff Switch to side-by-side view
Loading...
PKG_DEL Diff Switch to side-by-side view
Loading...
devel/polygons2d/drawPolyline.m to inst/polygons2d/drawPolyline.m
--- a/devel/polygons2d/drawPolyline.m
+++ b/inst/polygons2d/drawPolyline.m
@@ -26,106 +26,114 @@
 ## those of the authors and should not be interpreted as representing official
 ## policies, either expressed or implied, of the copyright holders.
 
-function varargout = drawPolyline(varargin)
-#DRAWPOLYLINE Draw a polyline specified by a list of points
-#
-#   drawPolyline(COORD);
-#   packs coordinates in a single [N*2] array.
-#
-#   drawPolyline(PX, PY);
-#   specifies coordinates in separate arrays. PX and PY must be column
-#   vectors with the same length.
-#
-#   drawPolyline(..., TYPE);
-#   where TYPE is either 'closed' or 'open', specifies if last point must
-#   be connected to the first one ('closed') or not ('open').
-#   Default is 'open'.
-#
-#   drawPolyline(..., PARAM, VALUE);
-#   specify plot options as described for plot command.
-#
-#   H = drawPolyline(...) also return a handle to the list of line objects.
-#
-#   Example:
-#   # Draw a curve representing an ellipse
-#   t = linspace(0, 2*pi, 100)';
-#   px = 10*cos(t); py = 5*sin(t);
-#   drawPolyline([px py], 'closed');
-#   axis equal;
-#
-#   # The same, with different drawing options
-#   drawPolyline([px py], 'closed', 'lineWidth', 2, 'lineStyle', '--');
-#
-#   See Also:
-#   polygons2d, drawPolygon
-#
-#   ---------
-#   author : David Legland 
-#   INRA - TPV URPOI - BIA IMASTE
-#   created the 06/04/2004.
-#
-
-#   HISTORY
-#   03/01/2007: better processing of input, and update doc (drawing
-#       options and CLOSE option)
-#   30/04/2009 rename as drawPolyline.
-
-
-# default values
-closed = false;
-
-# If first argument is a cell array, draw each curve individually,
-# and eventually returns handle of each plot.
-var = varargin{1};
-if iscell(var)
-    h = [];
-    for i=1:length(var(:))
-        h = [h ; drawPolyline(var{i}, varargin{2:end})];
-    end
-    if nargout>0
-        varargout{1}=h;
-    end
-    return;
-end
-
-# extract curve coordinate
-if size(var, 2)==1
-    # first argument contains x coord, second argument contains y coord
-    px = var;
-    if length(varargin)==1
-        error('Wrong number of arguments in drawPolyline');
-    end
-    py = varargin{2};
-    varargin = varargin(3:end);
-else
-    # first argument contains both coordinate
-    px = var(:, 1);
-    py = var(:, 2);
-    varargin = varargin(2:end);
-end
-
-# check if curve is closed or open
-if ~isempty(varargin)
-    var = varargin{1};
-    if strncmpi(var, 'close', 5)
-        closed = true;
-        varargin = varargin(2:end);
-    elseif strncmpi(var, 'open', 4)
-        closed = false;
-        varargin = varargin(2:end);
-    end
-end
-
-# if curve is closed, add first point at the end of the list
-if closed
-    px = [px; px(1)];
-    py = [py; py(1)];
-end
-
-# plot the curve, with eventually optional parameters
-h = plot(px, py, varargin{:});
-
-# format output arguments
-if nargout>0
-    varargout{1}=h;
-end
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{} =} drawPolyline (@var{}, @var{})
+##DRAWPOLYLINE Draw a polyline specified by a list of points
+##
+##   drawPolyline(COORD);
+##   packs coordinates in a single [N*2] array.
+##
+##   drawPolyline(PX, PY);
+##   specifies coordinates in separate arrays. PX and PY must be column
+##   vectors with the same length.
+##
+##   drawPolyline(..., TYPE);
+##   where TYPE is either 'closed' or 'open', specifies if last point must
+##   be connected to the first one ('closed') or not ('open').
+##   Default is 'open'.
+##
+##   drawPolyline(..., PARAM, VALUE);
+##   specify plot options as described for plot command.
+##
+##   H = drawPolyline(...) also return a handle to the list of line objects.
+##
+##   Example:
+##   @example
+##   # Draw a curve representing an ellipse
+##   t = linspace(0, 2*pi, 100)';
+##   px = 10*cos(t); py = 5*sin(t);
+##   drawPolyline([px py], 'closed');
+##   axis equal;
+##
+##   # The same, with different drawing options
+##   drawPolyline([px py], 'closed', 'lineWidth', 2, 'lineStyle', '--');
+##   @end example
+##
+## @seealso{polygons2d, drawPolygon}
+## @end deftypefn
+
+function varargout = drawPolyline(varargin)
+  # default values
+  closed = false;
+
+  # If first argument is a cell array, draw each curve individually,
+  # and eventually returns handle of each plot.
+  var = varargin{1};
+  if iscell(var)
+      h = [];
+      for i=1:length(var(:))
+          h = [h ; drawPolyline(var{i}, varargin{2:end})];
+      end
+      if nargout>0
+          varargout{1}=h;
+      end
+      return;
+  end
+
+  # extract curve coordinate
+  if size(var, 2)==1
+      # first argument contains x coord, second argument contains y coord
+      px = var;
+      if length(varargin)==1
+          error('Wrong number of arguments in drawPolyline');
+      end
+      py = varargin{2};
+      varargin = varargin(3:end);
+  else
+      # first argument contains both coordinate
+      px = var(:, 1);
+      py = var(:, 2);
+      varargin = varargin(2:end);
+  end
+
+  # check if curve is closed or open
+  if ~isempty(varargin)
+      var = varargin{1};
+      if strncmpi(var, 'close', 5)
+          closed = true;
+          varargin = varargin(2:end);
+      elseif strncmpi(var, 'open', 4)
+          closed = false;
+          varargin = varargin(2:end);
+      end
+  end
+
+  # if curve is closed, add first point at the end of the list
+  if closed
+      px = [px; px(1)];
+      py = [py; py(1)];
+  end
+
+  # plot the curve, with eventually optional parameters
+  h = plot(px, py, varargin{:});
+
+  # format output arguments
+  if nargout>0
+      varargout{1}=h;
+  end
+
+endfunction
+
+%!demo
+%! t  = linspace (0, 2*pi, 100)';
+%! px = 10*cos (t); py = 5*sin (t);
+%! 
+%! subplot (1,2,1)
+%! drawPolyline ([px py], 'closed');
+%! axis equal;
+%!
+%! subplot (1,2,2)
+%! drawPolyline([px py], 'closed', 'lineWidth', 2, 'lineStyle', '--');
+%! axis equal;
+%! # -------------------------------------------------
+%! # Draw a curve representing an ellipse with two drawing options
devel/polynomialCurves2d/polynomialDerivate.m to inst/polynomialCurves2d/drawPolynomialCurve.m
--- a/devel/polynomialCurves2d/polynomialDerivate.m
+++ b/inst/polynomialCurves2d/drawPolynomialCurve.m
@@ -23,35 +23,32 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function deriv = polynomialDerivate(poly)
-#POLYNOMIALDERIVATE Derivate a polynomial
-#
-#   DERIV = polynomialDERIVATE(POLY)
-#   POLY is a row vector of [n+1] coefficients, in the form:
-#       [a0 a1 a2 ... an]
-#   DERIV has the same format, with length n:
-#       [a1 a2*2 ... an*n]
-#
-#
-#   Example
-#   T = polynomialDerivate([2 3 4])
-#   returns:
-#   T = [3 8]
-#
-#   See also
-#   polynomialCurves2d
-#
-# ------
-# Author: David Legland
-# e-mail: david.legland@grignon.inra.fr
-# Created: 2007-02-23
-# Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
-
-
-# create the derivation matrices
-matrix = diag(0:length(poly)-1);
-
-# compute coefficients of first derivative polynomials
-deriv = circshift(poly*matrix, [0 -1]);
-
-
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{h} =} drawPolynomialCurve (@var{bnd}, @var{xcoefs}, @var{ycoesf})
+## @deftypefnx {Function File} {@var{h} =} drawPolynomialCurve (@dots{}, @var{npts})
+## Draw a polynomial curve approximation
+## @end deftypefn
+
+function varargout = drawPolynomialCurve(tBounds, xCoef, yCoef, nPts=120)
+  ## Extract input parameters
+
+  # parametrization bounds
+  t0 = tBounds(1);
+  t1 = tBounds(end);
+
+  ## Drawing the polyline approximation
+
+  # generate vector of absissa
+  t = linspace (t0, t1, nPts+1)';
+
+  # compute corresponding positions
+  pts = polynomialCurvePoint (t, xCoef, yCoef);
+
+  # draw the resulting curve
+  h = drawPolyline (pts);
+
+  if nargout > 0
+      varargout{1} = h;
+  end
+
+endfunction
devel/polynomialCurves2d/polynomialCurvePoint.m to inst/polynomialCurves2d/polynomialCurvePoint.m
--- a/devel/polynomialCurves2d/polynomialCurvePoint.m
+++ b/inst/polynomialCurves2d/polynomialCurvePoint.m
@@ -23,56 +23,47 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function point = polynomialCurvePoint(t, varargin)
-#POLYNOMIALCURVEPOINT Compute point corresponding to a position
-#
-#   POINT = polynomialCurvePoint(T, XCOEF, YCOEF)
-#   XCOEF and YCOEF are row vectors of coefficients, in the form:
-#       [a0 a1 a2 ... an]
-#   T is a either a scalar, or a column vector, containing values of the
-#   parametrization variable.
-#   POINT is a 1x2 array containing coordinate of point corresponding to
-#   position given by T. If T is a vector, POINT has as many rows as T.
-#
-#   POINT = polynomialCurvePoint(T, COEFS)
-#   COEFS is either a 2xN matrix (one row for the coefficients of each
-#   coordinate), or a cell array.
-#
-#   Example
-#   polynomialCurvePoint
-#
-#   See also
-#   polynomialCurves2d, polynomialCurveLength
-#
-#
-# ------
-# Author: David Legland
-# e-mail: david.legland@grignon.inra.fr
-# Created: 2007-02-23
-# Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
-
-## Extract input parameters
-
-# polynomial coefficients for each coordinate
-var = varargin{1};
-if iscell(var)
-    xCoef = var{1};
-    yCoef = var{2};
-elseif size(var, 1)==1
-    xCoef = varargin{1};
-    yCoef = varargin{2};
-else
-    xCoef = var(1,:);
-    yCoef = var(2,:);
-end
-    
-
-## compute length by numerical integration
-
-# convert polynomial coefficients to polyval convention
-cx = xCoef(end:-1:1);
-cy = yCoef(end:-1:1);
-
-# numerical integration of the Jacobian of parametrized curve
-point = [polyval(cx, t) polyval(cy, t)];
-
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{point} =} polynomialCurvePoint (@var{t}, @var{xcoef},@var{ycoef})
+## @deftypefnx {Function File} {@var{point} =} polynomialCurvePoint (@var{t}, @var{coefs})
+## Compute point corresponding to a position
+##
+##   @var{xcoef} and @var{ycoef} are row vectors of coefficients, in the form:
+##       [a0 a1 a2 ... an]
+##   @var{t} is a either a scalar, or a column vector, containing values of the
+##   parametrization variable.
+##   @var{point} is a 1x2 array containing coordinate of point corresponding to
+##   position given by @var{t}. If @var{t} is a vector, @var{point} has as many rows as @var{t}.
+##
+##   @var{coefs} is either a 2xN matrix (one row for the coefficients of each
+##   coordinate), or a cell array.
+## 
+## @seealso{polynomialCurves2d, polynomialCurveLength}
+## @end deftypefn
+
+function point = polynomialCurvePoint (t, varargin)
+  ## Extract input parameters
+
+  # polynomial coefficients for each coordinate
+  var = varargin{1};
+  if iscell (var)
+      xCoef = var{1};
+      yCoef = var{2};
+  elseif size (var, 1)==1
+      xCoef = varargin{1};
+      yCoef = varargin{2};
+  else
+      xCoef = var(1,:);
+      yCoef = var(2,:);
+  end
+
+  ## compute length by numerical integration
+
+  # convert polynomial coefficients to polyval convention
+  cx = xCoef(end:-1:1);
+  cy = yCoef(end:-1:1);
+
+  # numerical integration of the Jacobian of parametrized curve
+  point = [polyval(cx, t) polyval(cy, t)];
+
+endfunction
devel/polynomialCurves2d/polynomialCurveSetFit.m to inst/polynomialCurves2d/polynomialCurveSetFit.m
--- a/devel/polynomialCurves2d/polynomialCurveSetFit.m
+++ b/inst/polynomialCurves2d/polynomialCurveSetFit.m
@@ -23,220 +23,220 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-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    
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{coefs} @var{bnds}]=} polynomialCurveSetFit (@var{img})
+## @deftypefnx {Function File} {@dots{} =} polynomalCurveSetFit (@var{img}, @var{deg})
+## @deftypefnx {Function File} {[@dots{} @var{lbl}] =} polynomalCurveSetFit (@dots{})
+## Fit a set of polynomial curves to a segmented image
+##
+##   Result is a cell array of matrices. Each matrix is @var{deg}+1-by-2, and
+##   contains coefficients of polynomial curve for each coordinate. 
+##   @var{bnds} contains the boundary of the parametrizations.
+##   @var{img} is first binarised, then skeletonized.
+##
+##   Also returns an image of labels @var{lbl} for the segmented curves. The max label
+##   is the number of curves, and the length of @var{coefs}.
+##
+##   Requires the toolboxes:
+##   - Optimization
+##   - Image Processing
+##
+## @seealso{polynomialCurves2d, polynomialCurveFit}
+## @end deftypefn
+
+function [coefs bnds lblBranches] = polynomialCurveSetFit(seg, varargin)
+  # 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);
+  bnds  = cell(nBranches, 1);
+
+  # For each curve, find interpolated polynomial curve
+  for i = 1:nBranches
+      # 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, true (3,3))));
+          inds = inds(inds ~= 0);
+          
+          if length(inds) < 2
+              warning ("geometry:poylnomialCurveSetFit", ...
+                   ['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];
+          bnds{i}  = [0 1];
+          # 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];
+      bnds{i}  = t([1 end]);
+  end
+
+endfunction
+
+
+function points = chainPixels(img, varargin)
+#CHAINPIXELS return the list of points which constitute a curve on image
+#   output = chainPixels(input)
+
+  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
+
+endfunction
+
+%!demo
+%! [m, cmap] = imread ("default.img");
+%! m         = ind2gray (m, cmap);
+%! mbw       = im2bw(m, graythresh(m)*1.3);
+%!
+%! [c t] = polynomialCurveSetFit (mbw);
+%!
+%! figure(1)
+%! clf;
+%! imshow (m)
+%! hold on
+%! for i=1:numel(c)
+%!   if !isempty (c{i})
+%!     drawPolynomialCurve (t{i}, c{i}(1,:),c{i}(2,:));
+%!   endif
+%! endfor
+%!
+%! hold off
devel/polynomialCurves2d/drawPolynomialCurve.m to inst/polynomialCurves2d/polynomialCurveDerivative.m
--- a/devel/polynomialCurves2d/drawPolynomialCurve.m
+++ b/inst/polynomialCurves2d/polynomialCurveDerivative.m
@@ -23,67 +23,51 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function varargout = drawPolynomialCurve(tBounds, varargin)
-#DRAWPOLYNOMIALCURVE Draw a polynomial curve approximation
-#
-#   Usage
-#   drawPolynomialCurve(BND, XCOEFS, YCOEFS)
-#   drawPolynomialCurve(BND, XCOEFS, YCOEFS, NPTS)
-#
-#   Example
-#   drawPolynomialCurve
-#
-#   See also
-#
-#
-# ------
-# Author: David Legland
-# e-mail: david.legland@grignon.inra.fr
-# Created: 2011-03-21,    using Matlab 7.9.0.529 (R2009b)
-# Copyright 2011 INRA - Cepia Software Platform.
-
-
-## Extract input parameters
-
-# parametrization bounds
-t0 = tBounds(1);
-t1 = tBounds(end);
-
-# polynomial coefficients for each coordinate
-var = varargin{1};
-if iscell(var)
-    xCoef = var{1};
-    yCoef = var{2};
-    varargin(1) = [];
-    
-elseif size(var, 1)==1
-    xCoef = varargin{1};
-    yCoef = varargin{2};
-    varargin(1:2) = [];
-    
-else
-    xCoef = var(1,:);
-    yCoef = var(2,:);
-    varargin(1) = [];
-end
-
-nPts = 120;
-if ~isempty(varargin)
-    nPts = varargin{1};
-end
-
-
-## Drawing the polyline approximation
-
-# generate vector of absissa
-t = linspace(t0, t1, nPts+1)';
-
-# compute corresponding positions
-pts = polynomialCurvePoint(t, xCoef, yCoef);
-
-# draw the resulting curve
-h = drawPolyline(pts);
-
-if nargout > 0
-    varargout = {h};
-end
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{v} =} polynomialCurveDerivative (@var{t}, @var{xcoef},@var{ycoef})
+## @deftypefnx {Function File} {@var{v} =} polynomialCurveDerivative (@var{t}, @var{coefs})
+## Compute derivative vector of a polynomial curve
+##
+##   @var{xcoef} and @var{ycoef} are row vectors of coefficients, in the form:
+##       [a0 a1 a2 ... an]
+##   @var{v} is a 1x2 array containing direction of derivative of polynomial
+##   curve, computed for position @var{t}. If @var{t} is a vector, @var{v} has as many rows
+##   as the length of @var{t}.
+##
+##   @var{coefs} is either a 2xN matrix (one row for the coefficients of each
+##   coordinate), or a cell array.
+## 
+## @seealso{polynomialCurves2d, polynomialCurveNormal, polynomialCurvePoint,
+##   polynomialCurveCurvature}
+## @end deftypefn
+
+
+function v = polynomialCurveDerivative(t, varargin)
+  ## Extract input parameters
+
+  # polynomial coefficients for each coordinate
+  var = varargin{1};
+  if iscell(var)
+      xCoef = var{1};
+      yCoef = var{2};
+  elseif size(var, 1)==1
+      xCoef = varargin{1};
+      yCoef = varargin{2};
+  else
+      xCoef = var(1,:);
+      yCoef = var(2,:);
+  end
+  # convert to Octave polynomial convention
+  xCoef = xCoef(end:-1:1);
+  yCoef = yCoef(end:-1:1);
+      
+
+  # compute derivative of the polynomial
+  dx = polyder (xCoef);
+  dy = polyder (yCoef);
+
+  # numerical integration of the Jacobian of parametrized curve
+  v = [polyval(dx, t) polyval(dy, t)];
+
+endfunction
+
devel/polynomialCurves2d/polynomialCurveDerivative.m to inst/polynomialCurves2d/polynomialCurveCentroid.m
--- a/devel/polynomialCurves2d/polynomialCurveDerivative.m
+++ b/inst/polynomialCurves2d/polynomialCurveCentroid.m
@@ -23,60 +23,86 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function v = polynomialCurveDerivative(t, varargin)
-#POLYNOMIALCURVEDERIVATIVE Compute derivative vector of a polynomial curve
-#
-#   VECT = polynomialCurveLength(T, XCOEF, YCOEF);
-#   XCOEF and YCOEF are row vectors of coefficients, in the form:
-#       [a0 a1 a2 ... an]
-#   VECT is a 1x2 array containing direction of derivative of polynomial
-#   curve, computed for position T. If T is a vector, VECT has as many rows
-#   as the length of T.
-#
-#   VECT = polynomialCurveLength(T, COEFS);
-#   COEFS is either a 2xN matrix (one row for the coefficients of each
-#   coordinate), or a cell array.
-#
-#   Example
-#   polynomialCurveDerivative
-#
-#   See also
-#   polynomialCurves2d, polynomialCurveNormal, polynomialCurvePoint,
-#   polynomialCurveCurvature 
-#
-#
-# ------
-# Author: David Legland
-# e-mail: david.legland@grignon.inra.fr
-# Created: 2007-02-23
-# Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
-
-## Extract input parameters
-
-# polynomial coefficients for each coordinate
-var = varargin{1};
-if iscell(var)
-    xCoef = var{1};
-    yCoef = var{2};
-elseif size(var, 1)==1
-    xCoef = varargin{1};
-    yCoef = varargin{2};
-else
-    xCoef = var(1,:);
-    yCoef = var(2,:);
-end
-    
-
-## compute derivative
-
-# compute derivative of the polynomial
-dx = polynomialDerivate(xCoef);
-dy = polynomialDerivate(yCoef);
-
-# convert to polyval convention
-dx = dx(end:-1:1);
-dy = dy(end:-1:1);
-
-# numerical integration of the Jacobian of parametrized curve
-v = [polyval(dx, t) polyval(dy, t)];
-
+## -*- texinfo -*-
+## @deftypefn {Function File} {@var{c} =} polynomialCurveCentroid (@var{t}, @var{xcoef}, @var{ycoef})
+## @deftypefnx {Function File} {@var{c} =} polynomialCurveCentroid (@var{t}, @var{coefs})
+## @deftypefnx {Function File} {@var{c} =} polynomialCurveCentroid (@dots{}, @var{tol})
+## Compute the centroid of a polynomial curve
+##
+##   @var{xcoef} and @var{ycoef} are row vectors of coefficients, in the form:
+##       [a0 a1 a2 ... an]
+##   @var{t} is a 1x2 row vector, containing the bounds of the parametrization
+##   variable: @var{t} = [T0 T1], with @var{t} taking all values between T0 and T1.
+##   @var{c} contains coordinate of the polynomial curve centroid.
+##
+##   @var{coefs} is either a 2xN matrix (one row for the coefficients of each
+##   coordinate), or a cell array.
+##
+##   @var{tol} is the tolerance fo computation (absolute).
+## 
+## @seealso{polynomialCurves2d, polynomialCurveLength}
+## @end deftypefn
+
+function centroid = polynomialCurveCentroid(tBounds, varargin)
+  ## Extract input parameters
+
+  # parametrization bounds
+  t0 = tBounds(1);
+  t1 = tBounds(end);
+
+  # polynomial coefficients for each coordinate
+  var = varargin{1};
+  if iscell(var)
+      cx = var{1};
+      cy = var{2};
+      varargin(1) = [];
+  elseif size(var, 1)==1
+      cx = varargin{1};
+      cy = varargin{2};
+      varargin(1:2)=[];
+  else
+      cx = var(1,:);
+      cy = var(2,:);
+      varargin(1)=[];
+  end
+  # convert to Octave polyval format
+  cx = cx(end:-1:1);
+  cy = cy(end:-1:1);
+
+  # tolerance
+  tol = 1e-6;
+  if ~isempty(varargin)
+      tol = varargin{1};
+  end
+
+  ## compute length by numerical integration
+
+  # compute derivative of the polynomial
+  dx = polyder (cx);
+  dy = polyder (cy);
+
+  # compute curve length by integrating the Jacobian
+  L = quad(@(t)sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
+
+  # compute first coordinate of centroid
+  xc = quad(@(t)polyval(cx, t).*sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
+
+  # compute first coordinate of centroid
+  yc = quad(@(t)polyval(cy, t).*sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
+
+  # divide result of integration by total length of the curve
+  centroid = [xc yc]/L;
+
+endfunction
+
+%!demo
+%! bounds = [-1 1];
+%! coefs = [0 1 1; 0 -1 2];
+%! c = polynomialCurveCentroid (bounds, coefs);
+%!
+%! drawPolynomialCurve (bounds, coefs(1,:), coefs(2,:));
+%! hold on
+%! plot (c(1),c(2),'sr')
+%! hold off
+%! # -------------------------------------------------
+%! # Centriod of a polynomial curve
devel/polynomialCurves2d/polynomialCurveCentroid.m to inst/polynomialCurves2d/polynomialCurveFit.m
--- a/devel/polynomialCurves2d/polynomialCurveCentroid.m
+++ b/inst/polynomialCurves2d/polynomialCurveFit.m
@@ -23,86 +23,256 @@
 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
-function centroid = polynomialCurveCentroid(tBounds, varargin)
-#POLYNOMIALCURVECENTROID Compute the centroid of a polynomial curve
-#
-#   C = polynomialCurveCentroid(T, XCOEF, YCOEF)
-#   XCOEF and YCOEF are row vectors of coefficients, in the form:
-#       [a0 a1 a2 ... an]
-#   T is a 1x2 row vector, containing the bounds of the parametrization
-#   variable: T = [T0 T1], with T taking all values between T0 and T1.
-#   C contains coordinate of the polynomial curve centroid.
-#
-#   C = polynomialCurveCentroid(T, COEFS)
-#   COEFS is either a 2xN matrix (one row for the coefficients of each
-#   coordinate), or a cell array.
-#
-#   C = polynomialCurveCentroid(..., TOL)
-#   TOL is the tolerance fo computation (absolute).
-#
-#   Example
-#   polynomialCurveCentroid
-#
-#   See also
-#   polynomialCurves2d, polynomialCurveLength
-#
-#
-# ------
-# Author: David Legland
-# e-mail: david.legland@gignon.inra.fr
-# Created: 2007-02-23
-# Copyright 2007 INRA - BIA PV Nantes - MIAJ Jouy-en-Josas.
-
-
-## Extract input parameters
-
-# parametrization bounds
-t0 = tBounds(1);
-t1 = tBounds(end);
-
-# polynomial coefficients for each coordinate
-var = varargin{1};
-if iscell(var)
-    xCoef = var{1};
-    yCoef = var{2};
-    varargin(1) = [];
-elseif size(var, 1)==1
-    xCoef = varargin{1};
-    yCoef = varargin{2};
-    varargin(1:2)=[];
-else
-    xCoef = var(1,:);
-    yCoef = var(2,:);
-    varargin(1)=[];
-end
-
-# tolerance
-tol = 1e-6;
-if ~isempty(varargin)
-    tol = varargin{1};
-end
-
-## compute length by numerical integration
-
-# compute derivative of the polynomial
-dx = polynomialDerivate(xCoef);
-dy = polynomialDerivate(yCoef);
-
-# convert to polyval format
-dx = dx(end:-1:1);
-dy = dy(end:-1:1);
-
-cx = xCoef(end:-1:1);
-cy = yCoef(end:-1:1);
-
-# compute curve length by integrating the Jacobian
-L = quad(@(t)sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
-
-# compute first coordinate of centroid
-xc = quad(@(t)polyval(cx, t).*sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
-
-# compute first coordinate of centroid
-yc = quad(@(t)polyval(cy, t).*sqrt(polyval(dx, t).^2+polyval(dy, t).^2), t0, t1, tol);
-
-# divide result of integration by total length of the curve
-centroid = [xc yc]/L;
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{xc}, @var{yc}] =} polynomialCurveFit (@var{t}, @var{xt}, @var{yt}, @var{order})
+## @deftypefnx {Function File} {[@var{xc}, @var{yc}] =} polynomialCurveFit (@var{t}, @var{points}, @var{order})
+## @deftypefnx {Function File} {[@var{xc}, @var{yc}] =} polynomialCurveFit (@dots{}, @var{ti}, @var{condi})
+## Fit a polynomial curve to a series of points
+##
+##   @var{t} is a Nx1 vector
+##   @var{xt} and @var{yt} are coordinate for each parameter value (column vectors)
+##   @var{order} is the degree of the polynomial used for interpolation
+##   @var{xc} and @var{yc} are polynomial coefficients, given in @var{order}+1 row vectors,
+##   starting from degree 0 and up to degree @var{order}.
+##
+##   @var{points} specifies coordinate of points in a Nx2 array.
+##
+##   Impose some specific conditions using @var{ti} and @var{condi}. 
+##   @var{ti} is a value of the parametrization
+##   variable. @var{condi} is a cell array, with 2 columns, and as many rows as
+##   the derivatives specified for the given @var{ti}. Format for @var{condi} is:
+##   @var{condi} = @@{X_I, Y_I; X_I', Y_I'; X_I", Y_I"; ...@@};
+##   with X_I and Y_I being the imposed coordinate at position @var{ti}, X_I' and
+##   Y_I' being the imposed first derivatives, X_I" and Y_I" the imposed
+##   second derivatives, and so on...
+##   To specify a derivative without specifying derivative with lower
+##   degree, value of lower derivative can be let empty, using '[]'
+##
+##
+##   Requires the optimization Toolbox.
+##
+##   Run @command{demo polynomialCurveFit} to see exaples of use.
+## 
+## @seealso{polynomialCurves2d}
+## @end deftypefn
+
+function varargout = polynomialCurveFit(t, varargin)
+
+  ## extract input arguments
+
+  # extract curve coordinate
+  var = varargin{1};
+  if min(size(var))==1
+      # curve given as separate arguments
+      xt = varargin{1};
+      yt = varargin{2};
+      varargin(1:2)=[];
+  else
+      # curve coordinate bundled in a matrix
+      if size(var, 1)<size(var, 2)
+          var = var';
+      end
+      xt = var(:,1);
+      yt = var(:,2);
+      varargin(1)=[];
+  end
+
+  # order of the polynom
+  var = varargin{1};
+  if length(var)>1
+      Dx = var(1);
+      Dy = var(2);
+  else
+      Dx = var;
+      Dy = var;
+  end
+  varargin(1)=[];
+
+
+  ## Initialize local conditions
+
+  # For a solution vector 'x', the following relation must hold:
+  #   Aeq * x == beq,
+  # with:
+  #   Aeq   Matrix M*N 
+  #   beq   column vector with length M
+  # The coefficients of the Aeq matrix are initialized as follow:
+  # First point and last point are considered successively. For each point,
+  # k-th condition is the value of the (k-1)th derivative. This value is
+  # computed using relation of the form:
+  #   value = sum_i ( fact(i) * t_j^pow(i) )
+  # with:
+  #   i     indice of the (i-1) derivative. 
+  #   fact  row vector containing coefficient of each power of t, initialized
+  #       with a row vector equals to [1 1 ... 1], and updated for each
+  #       derivative by multiplying by corresponding power minus 1.
+  #   pow   row vector of the powers of each monome. It is represented by a
+  #       row vector containing an increasing series of power, eventually
+  #       completed with zeros for lower degrees (for the k-th derivative,
+  #       the coefficients with power lower than k are not relevant).
+
+  # Example for degree 5 polynom:
+  #   iter deriv  pow                 fact
+  #   1    0      [0 1 2 3 4 5]       [1 1 1 1 1 1]
+  #   2    1      [0 0 1 2 3 4]       [0 1 2 3 4 5]
+  #   3    2      [0 0 0 1 2 3]       [0 0 1 2 3 4]
+  #   4    3      [0 0 0 0 1 2]       [0 0 0 1 2 3]
+  #   ...
+  #   The process is repeated for coordinate x and for coordinate y.
+
+  # Initialize empty matrices
+  Aeqx = zeros(0, Dx+1);
+  beqx = zeros(0, 1);
+  Aeqy = zeros(0, Dy+1);
+  beqy = zeros(0, 1);
+
+  # Process local conditions
+  while ~isempty(varargin)
+      if length(varargin)==1
+          warning('MatGeom:PolynomialCurveFit:ArgumentNumber', ...
+              'Wrong number of arguments in polynomialCurvefit');
+      end
+
+      # extract parameter t, and cell array of local conditions
+      ti = varargin{1};
+      cond = varargin{2};
+
+      # factors for coefficients of each polynomial. At the beginning, they
+      # all equal 1. With successive derivatives, their value increase by the
+      # corresponding powers.
+      factX = ones(1, Dx+1);
+      factY = ones(1, Dy+1);
+
+      # start condition initialisations
+      for i = 1:size(cond, 1)
+          # degrees of each polynomial
+          powX = [zeros(1, i) 1:Dx+1-i];
+          powY = [zeros(1, i) 1:Dy+1-i];
+          
+          # update conditions for x coordinate
+          if ~isempty(cond{i,1})
+              Aeqx = [Aeqx ; factY.*power(ti, powX)]; ##ok<AGROW>
+              beqx = [beqx; cond{i,1}]; ##ok<AGROW>
+          end
+
+          # update conditions for y coordinate
+          if ~isempty(cond{i,2})
+              Aeqy = [Aeqy ; factY.*power(ti, powY)]; ##ok<AGROW>
+              beqy = [beqy; cond{i,2}]; ##ok<AGROW>
+          end
+          
+          # update polynomial degrees for next derivative
+          factX = factX.*powX;
+          factY = factY.*powY;
+      end
+      
+      varargin(1:2)=[];
+  end
+
+
+  ## Initialisations
+
+  # ensure column vectors
+  t  = t(:);
+  xt = xt(:);
+  yt = yt(:);
+
+  # number of points to fit
+  L = length(t);
+
+
+  ## Compute coefficients of each polynomial
+
+  # main matrix for x coordinate, size L*(degX+1)
+  T = ones(L, Dx+1);
+  for i = 1:Dx
+      T(:, i+1) = power(t, i);
+  end
+
+  # compute interpolation
+  # Octave compatibility - JPi 2013
+  xc = lsqlin (T, xt, zeros(1, Dx+1), 1, Aeqx, beqx)';
+  
+  # main matrix for y coordinate, size L*(degY+1)
+  T = ones(L, Dy+1);
+  for i = 1:Dy
+      T(:, i+1) = power(t, i);
+  end
+
+  # compute interpolation
+  # Octave compatibility - JPi 2013
+  yc = lsqlin (T, yt, zeros(1, Dy+1), 1, Aeqy, beqy)';
+
+
+  ## Format output arguments
+  if nargout <= 1
+      varargout{1} = {xc, yc};
+  else
+      varargout{1} = xc;
+      varargout{2} = yc;
+  end
+
+endfunction
+
+function x = lsqlin (C, d, A, b, Aeq, beq)
+  H = C'*C;
+  q = -C'*d;
+  x0 = zeros (size(C,2),size(d,2));
+  
+  x = qp (x0, H, q, Aeq, beq, [], [],[], A, b);
+endfunction
+
+%!demo
+%!   # defines a curve (circle arc) with small perturbations
+%!   N  = 50;
+%!   t  = linspace(0, 3*pi/4, N)';
+%!   xp = cos(t) + 5e-2*randn(size(t));
+%!   yp = sin(t) + 5e-2*randn(size(t));
+%!
+%!   [xc yc] = polynomialCurveFit(t, xp, yp, 3);
+%!
+%!   figure(1);
+%!   clf;
+%!   drawPolynomialCurve(t([1 end]), xc, yc);
+%!   hold on
+%!   plot(xp,yp,'.g');
+%!   hold off
+%!   axis tight
+%!   axis equal
+
+%!demo
+%!   # defines a curve (circle arc) with small perturbations
+%!   N  = 100;
+%!   t  = linspace(0, 3*pi/4, N)';
+%!   xp = cos(t) + 7e-2*randn(size(t)); 
+%!   yp = sin(t) + 7e-2*randn(size(t));
+%!
+%!   # plot the points
+%!   figure (1); clf; hold on;
+%!   axis ([-1.2 1.2 -.2 1.2]); axis equal;
+%!   drawPoint (xp, yp, ".g");
+%!
+%!   # fit without knowledge on bounds
+%!   [xc0 yc0] = polynomialCurveFit (t, xp, yp, 5);
+%!   h = drawPolynomialCurve (t([1 end]), xc0, yc0);
+%!   set(h, "color", "b")
+%!
+%!   # fit by imposing coordinate on first point
+%!   [xc1 yc1] = polynomialCurveFit (t, xp, yp, 5, 0, {1, 0});
+%!   h = drawPolynomialCurve (t([1 end]), xc1, yc1);
+%!   set(h, "color", "r")
+%!
+%!   # fit by imposing coordinate (1,0) and derivative (0,1) on first point
+%!   [xc2 yc2] = polynomialCurveFit (t, xp, yp, 5, 0, {1, 0;0 1});
+%!   h = drawPolynomialCurve (t([1 end]), xc2, yc2);
+%!   set(h, "color", "m")
+%!
+%!   # fit by imposing several conditions on various points
+%!   [xc3 yc3] = polynomialCurveFit (t, xp, yp, 5, ...
+%!       0, {1, 0;0 1}, ...      # coord and first derivative of first point
+%!       3*pi/4, {-sqrt(2)/2, sqrt(2)/2}, ...    # coord of last point
+%!       pi/2, {[], [];-1, 0});      # derivative of point on the top of arc
+%!   h = drawPolynomialCurve (t([1 end]), xc3, yc3);
+%!   set(h, "color", "k")
+%!   axis tight
+%!   axis equal