octave_packages/geometry-1.5.0/polygons2d/curvature.m

   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.
   4 %%
   5 %% Redistribution and use in source and binary forms, with or without
   6 %% modification, are permitted provided that the following conditions are met:
   7 %%
   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.
  13 %%
  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.
  24 %%
  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.
  28
  29 %% -*- texinfo -*-
  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.
  36 %%
  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
  39 %% for each point.
  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.
  44 %%
  45 %% If the argument @var{t} is not given, the parametrization is estimated using
  46 %% function @code{parametrize}.
  47 %%
  48 %% If requested, @var{poly} contains the approximating polygon evlauted at the
  49 %% parametrization @var{t}.
  50 %%
  51 %% @seealso{parametrize, polygons2d}
  52 %% @end deftypefn
  53
  54 function [kappa, varargout] = curvature(varargin)
  55
  56   % default values
  57   degree = 5;
  58   t=0;                    % parametrization of curve
  59   tc=0;                   % indices of points wished for curvature
  60
  61
  62   % =================================================================
  63
  64   % Extract method and degree ------------------------------
  65
  66   nargin = length(varargin);
  67   varN = varargin{nargin};
  68   varN2 = varargin{nargin-1};
  69
  70   if ischar(varN2)
  71       % method and degree are specified
  72       method = varN2;
  73       degree = varN;
  74       varargin = varargin(1:nargin-2);
  75   elseif ischar(varN)
  76       % only method is specified, use degree 6 as default
  77       method = varN;
  78       varargin = varargin{1:nargin-1};
  79   else
  80       % method and degree are implicit : use 'polynom' and 6
  81       method = 'polynom';
  82   end
  83
  84   % extract input parametrization and curve. -----------------------
  85   nargin = length(varargin);
  86   if nargin==1
  87       % parameters are just the points -> compute caracterization.
  88       var = varargin{1};
  89       px = var(:,1);
  90       py = var(:,2);
  91   elseif nargin==2
  92       var = varargin{2};
  93       if size(var, 2)==2
  94           % parameters are t and POINTS
  95           px = var(:,1);
  96           py = var(:,2);
  97           t = varargin{1};
  98       else
  99           % parameters are px and py
 100           px = varargin{1};
 101           py = var;
 102       end
 103   elseif nargin==3
 104       var = varargin{2};
 105       if size(var, 2)==2
 106           % parameters are t, POINTS, and tc
 107           px = var(:,1);
 108           py = var(:,2);
 109           t = varargin{1};
 110       else
 111           % parameters are t, px and py
 112           t = varargin{1};
 113           px = var;
 114           py = varargin{3};
 115       end
 116   elseif nargin==4
 117       % parameters are t, px, py and tc
 118       t  = varargin{1};
 119       px = varargin{2};
 120       py = varargin{3};
 121       tc = varargin{4};
 122   end
 123
 124   % compute implicit parameters --------------------------
 125
 126   % if t and/or tc are not computed, use implicit definition
 127   if t==0
 128       t = parametrize(px, py, 'norm');
 129   end
 130
 131   % if tc not defined, compute curvature for all points
 132   if tc==0
 133       tc = t;
 134   else
 135       % else convert from indices to parametrization values
 136       tc = t(tc);
 137   end
 138
 139
 140   % =================================================================
 141   %    compute curvature for each point of the curve
 142
 143   if strcmp(method, 'polynom')
 144       % compute coefficients of interpolation functions
 145       x0 = polyfit(t, px, degree);
 146       y0 = polyfit(t, py, degree);
 147
 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]);
 156
 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);
 162
 163       % compute value of curvature
 164       kappa = (xprime.*ysec - xsec.*yprime)./ ...
 165           power(xprime.*xprime + yprime.*yprime, 3/2);
 166
 167       if nargout > 1
 168         varargout{1} = [polyval(x0,tc(:)) polyval(y0,tc(:))];
 169         if nargout > 2
 170           varargout{2} = tc;
 171         end
 172       end
 173   else
 174       error('unknown method');
 175   end
 176
 177 endfunction