octave_packages/splines-1.0.7/catmullrom.m

   1 ## Copyright (C) 2008  Carlo de Falco
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; if not, see <http://www.fsf.org/licenses>
  15
  16 ## -*- texinfo -*-
  17 ##  @deftypefn {Function File} {@var{pp}} = catmullrom( @var{x},@
  18 ##  @var{f}, @var{v})
  19 ##
  20 ## Returns the piecewise polynomial form of the Catmull-Rom cubic
  21 ## spline interpolating @var{f} at the points @var{x}.
  22 ## If the input @var{v} is supplied it will be interpreted as the
  23 ## values of the tangents at the extremals, if it is
  24 ## missing, the values will be computed from the data via one-sided
  25 ## finite difference formulas. See the wikipedia page for "Cubic
  26 ## Hermite spline" for a description of the algorithm.
  27 ##
  28 ##  @seealso{ppval}
  29 ##  @end deftypefn
  30
  31 function pp = catmullrom(x,f,v)
  32
  33   if ( nargin < 2 )
  34     print_usage();
  35   endif
  36
  37   h00 = [2 -3 0 1];
  38   h10 = [1 -2 1 0];
  39   h01 = [-2 3 0 0];
  40   h11 = [1 -1 0 0];
  41
  42   h  = diff(x(:)');
  43   p0 = f(:)'(1:end-1);
  44   p1 = f(:)'(2:end);
  45
  46   if (nargin < 3)
  47     v(1) = (p1(1)-p0(1))./h(1);
  48     v(2) = (p1(end)-p0(end))./h(end);
  49   endif
  50   m  = (p1(2:end)-p0(1:end-1))./(h(2:end)+h(1:end-1));
  51   m0 = [v(1) m];
  52   m1 = [m v(2)];
  53
  54   for ii = 1:4
  55     coeff(:,ii) =  ((h00(ii)*p0 + h10(ii)*h.*m0 +...
  56                      h01(ii)*p1 + h11(ii)*h.*m1 )./h.^(4-ii))' ;
  57   end
  58
  59   pp = mkpp (x, coeff);
  60
  61 endfunction