1 %% Copyright (C) 2003-2011 David Legland <david.legland@grignon.inra.fr>
2 %% Copyright (C) 2012 Adapted to Octave by Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
3 %% All rights reserved.
5 %% Redistribution and use in source and binary forms, with or without
6 %% modification, are permitted provided that the following conditions are met:
8 %% 1 Redistributions of source code must retain the above copyright notice,
9 %% this list of conditions and the following disclaimer.
10 %% 2 Redistributions in binary form must reproduce the above copyright
11 %% notice, this list of conditions and the following disclaimer in the
12 %% documentation and/or other materials provided with the distribution.
14 %% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
15 %% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 %% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 %% ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
18 %% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 %% DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
20 %% SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
21 %% CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
22 %% OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
23 %% OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 %% The views and conclusions contained in the software and documentation are
26 %% those of the authors and should not be interpreted as representing official
27 %% policies, either expressed or implied, of the copyright holders.
30 %% @deftypefn {Function File} {@var{kappa} = } curvature (@var{t}, @var{px}, @var{py},@var{method},@var{degree})
31 %% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{t}, @var{poly},@var{method},@var{degree})
32 %% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{px}, @var{py},@var{method},@var{degree})
33 %% @deftypefnx {Function File} {@var{kappa} = } curvature (@var{points},@var{method},@var{degree})
34 %% @deftypefnx {Function File} {[@var{kappa} @var{poly} @var{t}] = } curvature (@dots{})
35 %% Estimate curvature of a polyline defined by points.
37 %% First compute an approximation of the curve given by PX and PY, with
38 %% the parametrization @var{t}. Then compute the curvature of approximated curve
40 %% @var{method} used for approximation can be only: 'polynom', with specified degree.
41 %% Further methods will be provided in a future version.
42 %% @var{t}, @var{px}, and @var{py} are N-by-1 array of the same length. The points
43 %% can be specified as a single N-by-2 array.
45 %% If the argument @var{t} is not given, the parametrization is estimated using
46 %% function @code{parametrize}.
48 %% If requested, @var{poly} contains the approximating polygon evlauted at the
49 %% parametrization @var{t}.
51 %% @seealso{parametrize, polygons2d}
54 function [kappa, varargout] = curvature(varargin)
58 t=0; % parametrization of curve
59 tc=0; % indices of points wished for curvature
62 % =================================================================
64 % Extract method and degree ------------------------------
66 nargin = length(varargin);
67 varN = varargin{nargin};
68 varN2 = varargin{nargin-1};
71 % method and degree are specified
74 varargin = varargin(1:nargin-2);
76 % only method is specified, use degree 6 as default
78 varargin = varargin{1:nargin-1};
80 % method and degree are implicit : use 'polynom' and 6
84 % extract input parametrization and curve. -----------------------
85 nargin = length(varargin);
87 % parameters are just the points -> compute caracterization.
94 % parameters are t and POINTS
99 % parameters are px and py
106 % parameters are t, POINTS, and tc
111 % parameters are t, px and py
117 % parameters are t, px, py and tc
124 % compute implicit parameters --------------------------
126 % if t and/or tc are not computed, use implicit definition
128 t = parametrize(px, py, 'norm');
131 % if tc not defined, compute curvature for all points
135 % else convert from indices to parametrization values
140 % =================================================================
141 % compute curvature for each point of the curve
143 if strcmp(method, 'polynom')
144 % compute coefficients of interpolation functions
145 x0 = polyfit(t, px, degree);
146 y0 = polyfit(t, py, degree);
148 % compute coefficients of first and second derivatives. In the case of a
149 % polynom, it is possible to compute coefficient of derivative by
150 % multiplying with a matrix.
151 derive = diag(degree:-1:0);
152 xp = circshift(x0*derive, [0 1]);
153 yp = circshift(y0*derive, [0 1]);
154 xs = circshift(xp*derive, [0 1]);
155 ys = circshift(yp*derive, [0 1]);
157 % compute values of first and second derivatives for needed points
158 xprime = polyval(xp, tc);
159 yprime = polyval(yp, tc);
160 xsec = polyval(xs, tc);
161 ysec = polyval(ys, tc);
163 % compute value of curvature
164 kappa = (xprime.*ysec - xsec.*yprime)./ ...
165 power(xprime.*xprime + yprime.*yprime, 3/2);
168 varargout{1} = [polyval(x0,tc(:)) polyval(y0,tc(:))];
174 error('unknown method');