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 --------- author : David Legland INRA - TPV URPOI - BIA IMASTE created the 07/04/2003.
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 % INRA - TPV URPOI - BIA IMASTE 0030 % created the 07/04/2003. 0031 % 0032 0033 0034 % default values 0035 degree = 5; 0036 t=0; % parametrization of curve 0037 tc=0; % indices of points wished for curvature 0038 0039 0040 % ================================================================= 0041 0042 % Extract method and degree ------------------------------ 0043 0044 nargin = length(varargin); 0045 varN = varargin{nargin}; 0046 varN2 = varargin{nargin-1}; 0047 0048 if ischar(varN2) 0049 % method and degree are specified 0050 method = varN2; 0051 degree = varN; 0052 varargin = varargin(1:nargin-2); 0053 elseif ischar(varN) 0054 % only method is specified, use degree 6 as default 0055 method = varN; 0056 varargin = varargin{1:nargin-1}; 0057 else 0058 % method and degree are implicit : use 'polynom' and 6 0059 method = 'polynom'; 0060 end 0061 0062 % extract input parametrization and curve. ----------------------- 0063 nargin = length(varargin); 0064 if nargin==1 0065 % parameters are just the points -> compute caracterization. 0066 var = varargin{1}; 0067 px = var(:,1); 0068 py = var(:,2); 0069 elseif nargin==2 0070 var = varargin{2}; 0071 if size(var, 2)==2 0072 % parameters are t and POINTS 0073 px = var(:,1); 0074 py = var(:,2); 0075 t = varargin{1}; 0076 else 0077 % parameters are px and py 0078 px = varargin{1}; 0079 py = var; 0080 end 0081 elseif nargin==3 0082 var = varargin{2}; 0083 if size(var, 2)==2 0084 % parameters are t, POINTS, and tc 0085 px = var(:,1); 0086 py = var(:,2); 0087 t = varargin{1}; 0088 else 0089 % parameters are t, px and py 0090 t = varargin{1}; 0091 px = var; 0092 py = varargin{3}; 0093 end 0094 elseif nargin==4 0095 % parameters are t, px, py and tc 0096 t = varargin{1}; 0097 px = varargin{2}; 0098 py = varargin{3}; 0099 tc = varargin{4}; 0100 end 0101 0102 % compute implicit parameters -------------------------- 0103 0104 % if t and/or tc are not computed, use implicit definition 0105 if t==0 0106 t = parametrize(px, py); 0107 t = t/t(length(t)); % normalize between 0 and 1 0108 end 0109 0110 % if tc not defined, compute curvature for all points 0111 if tc==0 0112 tc = t; 0113 else 0114 % else convert from indices to parametrization values 0115 tc = t(tc); 0116 end 0117 0118 0119 % ================================================================= 0120 % compute curvature for each point of the curve 0121 0122 if strcmp(method, 'polynom') 0123 % compute coefficients of interpolation functions 0124 x0 = polyfit(t, px, degree); 0125 y0 = polyfit(t, py, degree); 0126 0127 % compute coefficients of first and second derivatives. In the case of a 0128 % polynom, it is possible to compute coefficient of derivative by 0129 % multiplying with a matrix. 0130 derive = diag(degree:-1:0); 0131 xp = circshift(x0*derive, [0 1]); 0132 yp = circshift(y0*derive, [0 1]); 0133 xs = circshift(xp*derive, [0 1]); 0134 ys = circshift(yp*derive, [0 1]); 0135 0136 % compute values of first and second derivatives for needed points 0137 xprime = polyval(xp, tc); 0138 yprime = polyval(yp, tc); 0139 xsec = polyval(xs, tc); 0140 ysec = polyval(ys, tc); 0141 0142 % compute value of curvature 0143 kappa = (xprime.*ysec - xsec.*yprime)./ ... 0144 power(xprime.*xprime + yprime.*yprime, 3/2); 0145 else 0146 error('unknown method'); 0147 end