curvature

CURVATURE Estimate curvature of a polyline defined by points.

   KAPPA = curvature(T, PX, PY, METHOD, DEGREE)
   First compute an approximation of the curve given by PX and PY, with
   the parametrization T. METHOD used for approximation can be only:
   'polynom', with specified degree
   Further methods will be provided in a future version.
   T, PX, and PY are N*1 array of the same length.
   Then compute the curvature of approximated curve for each point.

   For example:
   KAPPA = curvature(t, px, py, 'polynom', 6)

   KAPPA = curvature(T, POINTS, METHOD, DEGREE)
   specify curve as a suite of points. POINTS is size [N*2].

   KAPPA = curvature(PX, PY, METHOD, DEGREE)
   KAPPA = curvature(POINTS, METHOD, DEGREE)
   compute implicite normalization of the curve, based on euclidian
   distance between 2 consecutive points, and normalized between 0 and 1.


   See also
   polygons2d, parametrize

0001 function kappa = curvature(varargin)
0002 %CURVATURE Estimate curvature of a polyline defined by points.
0003 %
0004 %   KAPPA = curvature(T, PX, PY, METHOD, DEGREE)
0005 %   First compute an approximation of the curve given by PX and PY, with
0006 %   the parametrization T. METHOD used for approximation can be only:
0007 %   'polynom', with specified degree
0008 %   Further methods will be provided in a future version.
0009 %   T, PX, and PY are N*1 array of the same length.
0010 %   Then compute the curvature of approximated curve for each point.
0011 %
0012 %   For example:
0013 %   KAPPA = curvature(t, px, py, 'polynom', 6)
0014 %
0015 %   KAPPA = curvature(T, POINTS, METHOD, DEGREE)
0016 %   specify curve as a suite of points. POINTS is size [N*2].
0017 %
0018 %   KAPPA = curvature(PX, PY, METHOD, DEGREE)
0019 %   KAPPA = curvature(POINTS, METHOD, DEGREE)
0020 %   compute implicite normalization of the curve, based on euclidian
0021 %   distance between 2 consecutive points, and normalized between 0 and 1.
0022 %
0023 %
0024 %   See also
0025 %   polygons2d, parametrize
0026 
0027 % ------
0028 % Author: David Legland
0029 % E-mail: david.legland@inrae.fr
0030 % Created: 2003-04-07
0031 % Copyright 2003-2024 INRA - TPV URPOI - BIA IMASTE
0032 
0033 % default values
0034 degree = 5;
0035 t=0;                    % parametrization of curve
0036 tc=0;                   % indices of points wished for curvature
0037 
0038 % Extract method and degree
0039 
0040 nargin = length(varargin);
0041 varN = varargin{nargin};
0042 varN2 = varargin{nargin-1};
0043 
0044 if ischar(varN2)
0045     % method and degree are specified
0046     method = varN2;
0047     degree = varN;
0048     varargin = varargin(1:nargin-2);
0049 elseif ischar(varN)
0050     % only method is specified, use degree 6 as default
0051     method = varN;
0052     varargin = varargin{1:nargin-1};
0053 else
0054     % method and degree are implicit : use 'polynom' and 6
0055     method = 'polynom';
0056 end
0057 
0058 % extract input parametrization and curve
0059 nargin = length(varargin);
0060 if nargin==1
0061     % parameters are just the points -> compute caracterization.
0062     var = varargin{1};
0063     px = var(:,1);
0064     py = var(:,2);
0065 elseif nargin==2
0066     var = varargin{2};
0067     if size(var, 2)==2
0068         % parameters are t and POINTS
0069         px = var(:,1);
0070         py = var(:,2);
0071         t = varargin{1};
0072     else
0073         % parameters are px and py
0074         px = varargin{1};
0075         py = var;
0076     end
0077 elseif nargin==3
0078     var = varargin{2};
0079     if size(var, 2)==2
0080         % parameters are t, POINTS, and tc
0081         px = var(:,1);
0082         py = var(:,2);
0083         t = varargin{1};
0084     else
0085         % parameters are t, px and py
0086         t = varargin{1};
0087         px = var;
0088         py = varargin{3};
0089     end
0090 elseif nargin==4
0091     % parameters are t, px, py and tc
0092     t  = varargin{1};
0093     px = varargin{2};
0094     py = varargin{3};
0095     tc = varargin{4};
0096 end
0097 
0098 % compute implicit parameters
0099 
0100 % if t and/or tc are not computed, use implicit definition
0101 if t==0
0102     t = parametrize(px, py);
0103     t = t/t(length(t)); % normalize between 0 and 1
0104 end
0105 
0106 % if tc not defined, compute curvature for all points
0107 if tc==0
0108     tc = t;
0109 else
0110     % else convert from indices to parametrization values
0111     tc = t(tc);
0112 end
0113  
0114 %% compute curvature for each point of the curve
0115 
0116 if strcmp(method, 'polynom')
0117     % compute coefficients of interpolation functions
0118     x0 = polyfit(t, px, degree);
0119     y0 = polyfit(t, py, degree);
0120 
0121     % compute coefficients of first and second derivatives. In the case of a
0122     % polynom, it is possible to compute coefficient of derivative by
0123     % multiplying with a matrix.
0124     derive = diag(degree:-1:0);
0125     xp = circshift(x0*derive, [0 1]);
0126     yp = circshift(y0*derive, [0 1]);
0127     xs = circshift(xp*derive, [0 1]);
0128     ys = circshift(yp*derive, [0 1]);
0129 
0130     % compute values of first and second derivatives for needed points
0131     xprime = polyval(xp, tc);
0132     yprime = polyval(yp, tc);
0133     xsec = polyval(xs, tc);
0134     ysec = polyval(ys, tc);
0135 
0136     % compute value of curvature
0137     kappa = (xprime.*ysec - xsec.*yprime)./ ...
0138         power(xprime.*xprime + yprime.*yprime, 3/2);
0139 else
0140     error('unknown method');
0141 end