hwpfilter/source/cspline.cpp

   1 /*************************************************************************
   2  *
   3  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
   4  *
   5  * Copyright 2008 by Sun Microsystems, Inc.
   6  *
   7  * OpenOffice.org - a multi-platform office productivity suite
   8  *
   9  * $RCSfile: cspline.cpp,v $
  10  * $Revision: 1.3 $
  11  *
  12  * This file is part of OpenOffice.org.
  13  *
  14  * OpenOffice.org is free software: you can redistribute it and/or modify
  15  * it under the terms of the GNU Lesser General Public License version 3
  16  * only, as published by the Free Software Foundation.
  17  *
  18  * OpenOffice.org is distributed in the hope that it will be useful,
  19  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  20  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21  * GNU Lesser General Public License version 3 for more details
  22  * (a copy is included in the LICENSE file that accompanied this code).
  23  *
  24  * You should have received a copy of the GNU Lesser General Public License
  25  * version 3 along with OpenOffice.org.  If not, see
  26  * <http://www.openoffice.org/license.html>
  27  * for a copy of the LGPLv3 License.
  28  *
  29  ************************************************************************/
  30
  31 // Natural, Clamped, or Periodic Cubic Splines
  32 //
  33 // Input:  A list of N+1 points (x_i,a_i), 0 <= i <= N, which are sampled
  34 // from a function, a_i = f(x_i).  The function f is unknown.  Boundary
  35 // conditions are
  36 //   (1) Natural splines:  f"(x_0) = f"(x_N) = 0
  37 //   (2) Clamped splines:  f'(x_0) and f'(x_N) are user-specified.
  38 //   (3) Periodic splines:  f(x_0) = f(x_N) [in which case a_N = a_0 is
  39 //       required in the input], f'(x_0) = f'(x_N), and f"(x_0) = f"(x_N).
  40 //
  41 // Output: b_i, c_i, d_i, 0 <= i <= N-1, which are coefficients for the cubic
  42 // spline S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3 for
  43 // x_i <= x < x_{i+1}.
  44 //
  45 // The natural and clamped algorithms were implemented from
  46 //
  47 //    Numerical Analysis, 3rd edition
  48 //    Richard L. Burden and J. Douglas Faires
  49 //    Prindle, Weber & Schmidt
  50 //    Boston, 1985, pp. 122-124.
  51 //
  52 // The algorithm sets up a tridiagonal linear system of equations in the
  53 // c_i values.  This can be solved in O(N) time.
  54 //
  55 // The periodic spline algorithm was implemented from my own derivation.  The
  56 // linear system of equations is not tridiagonal.  For now I use a standard
  57 // linear solver that does not take advantage of the sparseness of the
  58 // matrix.  Therefore for very large N, you may have to worry about memory
  59 // usage.
  60
  61 #include "solver.h"
  62 //-----------------------------------------------------------------------------
  63 void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
  64     double*& d)
  65 {
  66   const double oneThird = 1.0/3.0;
  67
  68   int i;
  69   double* h = new double[N];
  70   double* hdiff = new double[N];
  71   double* alpha = new double[N];
  72
  73   for (i = 0; i < N; i++){
  74     h[i] = x[i+1]-x[i];
  75   }
  76
  77   for (i = 1; i < N; i++)
  78     hdiff[i] = x[i+1]-x[i-1];
  79
  80   for (i = 1; i < N; i++)
  81   {
  82     double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
  83     double denom = h[i-1]*h[i];
  84     alpha[i] = numer/denom;
  85   }
  86
  87   double* ell = new double[N+1];
  88   double* mu = new double[N];
  89   double* z = new double[N+1];
  90   double recip;
  91
  92   ell[0] = 1.0;
  93   mu[0] = 0.0;
  94   z[0] = 0.0;
  95
  96   for (i = 1; i < N; i++)
  97   {
  98     ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
  99     recip = 1.0/ell[i];
 100     mu[i] = recip*h[i];
 101     z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
 102   }
 103   ell[N] = 1.0;
 104   z[N] = 0.0;
 105
 106   b = new double[N];
 107   c = new double[N+1];
 108   d = new double[N];
 109
 110   c[N] = 0.0;
 111
 112   for (i = N-1; i >= 0; i--)
 113   {
 114     c[i] = z[i]-mu[i]*c[i+1];
 115     recip = 1.0/h[i];
 116     b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
 117     d[i] = oneThird*recip*(c[i+1]-c[i]);
 118   }
 119
 120   delete[] h;
 121   delete[] hdiff;
 122   delete[] alpha;
 123   delete[] ell;
 124   delete[] mu;
 125   delete[] z;
 126 }
 127
 128 void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
 129     double*& d)
 130 {
 131   double* h = new double[N];
 132   int i;
 133   for (i = 0; i < N; i++)
 134     h[i] = x[i+1]-x[i];
 135
 136   mgcLinearSystemD sys;
 137   double** mat = sys.NewMatrix(N+1);  // guaranteed to be zeroed memory
 138   c = sys.NewVector(N+1);   // guaranteed to be zeroed memory
 139
 140   // c[0] - c[N] = 0
 141   mat[0][0] = +1.0f;
 142   mat[0][N] = -1.0f;
 143
 144   // h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
 145   //   3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
 146   for (i = 1; i <= N-1; i++)
 147   {
 148     mat[i][i-1] = h[i-1];
 149     mat[i][i  ] = 2.0f*(h[i-1]+h[i]);
 150     mat[i][i+1] = h[i];
 151     c[i] = 3.0f*((a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]);
 152   }
 153
 154   // "wrap around equation" for periodicity
 155   // h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
 156   //   3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
 157   mat[N][N-1] = h[N-1];
 158   mat[N][0  ] = 2.0f*(h[N-1]+h[0]);
 159   mat[N][1  ] = h[0];
 160   c[N] = 3.0f*((a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]);
 161
 162   // solve for c[0] through c[N]
 163   sys.Solve(N+1,mat,c);
 164
 165   const double oneThird = 1.0/3.0;
 166   b = new double[N];
 167   d = new double[N];
 168   for (i = 0; i < N; i++)
 169   {
 170     b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
 171     d[i] = oneThird*(c[i+1]-c[i])/h[i];
 172   }
 173
 174   sys.DeleteMatrix(N+1,mat);
 175 }