hwpfilter/source/cspline.cxx

   1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
   2 /*
   3  * This file is part of the LibreOffice project.
   4  *
   5  * This Source Code Form is subject to the terms of the Mozilla Public
   6  * License, v. 2.0. If a copy of the MPL was not distributed with this
   7  * file, You can obtain one at http://mozilla.org/MPL/2.0/.
   8  *
   9  * This file incorporates work covered by the following license notice:
  10  *
  11  *   Licensed to the Apache Software Foundation (ASF) under one or more
  12  *   contributor license agreements. See the NOTICE file distributed
  13  *   with this work for additional information regarding copyright
  14  *   ownership. The ASF licenses this file to you under the Apache
  15  *   License, Version 2.0 (the "License"); you may not use this file
  16  *   except in compliance with the License. You may obtain a copy of
  17  *   the License at http://www.apache.org/licenses/LICENSE-2.0 .
  18  */
  19
  20 // Natural, Clamped, or Periodic Cubic Splines
  21 //
  22 // Input:  A list of N+1 points (x_i,a_i), 0 <= i <= N, which are sampled
  23 // from a function, a_i = f(x_i).  The function f is unknown.  Boundary
  24 // conditions are
  25 //   (1) Natural splines:  f"(x_0) = f"(x_N) = 0
  26 //   (2) Clamped splines:  f'(x_0) and f'(x_N) are user-specified.
  27 //   (3) Periodic splines:  f(x_0) = f(x_N) [in which case a_N = a_0 is
  28 //       required in the input], f'(x_0) = f'(x_N), and f"(x_0) = f"(x_N).
  29 //
  30 // Output: b_i, c_i, d_i, 0 <= i <= N-1, which are coefficients for the cubic
  31 // spline S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3 for
  32 // x_i <= x < x_{i+1}.
  33 //
  34 // The natural and clamped algorithms were implemented from
  35 //
  36 //    Numerical Analysis, 3rd edition
  37 //    Richard L. Burden and J. Douglas Faires
  38 //    Prindle, Weber & Schmidt
  39 //    Boston, 1985, pp. 122-124.
  40 //
  41 // The algorithm sets up a tridiagonal linear system of equations in the
  42 // c_i values.  This can be solved in O(N) time.
  43 //
  44 // The periodic spline algorithm was implemented from my own derivation.  The
  45 // linear system of equations is not tridiagonal.  For now I use a standard
  46 // linear solver that does not take advantage of the sparseness of the
  47 // matrix.  Therefore for very large N, you may have to worry about memory
  48 // usage.
  49
  50 #include <sal/config.h>
  51
  52 #include "cspline.h"
  53 #include "solver.h"
  54
  55 void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
  56     double*& d)
  57 {
  58   const double oneThird = 1.0/3.0;
  59
  60   int i;
  61   double* h = new double[N];
  62   double* hdiff = new double[N];
  63   double* alpha = new double[N];
  64
  65   for (i = 0; i < N; i++){
  66     h[i] = x[i+1]-x[i];
  67   }
  68
  69   for (i = 1; i < N; i++)
  70     hdiff[i] = x[i+1]-x[i-1];
  71
  72   for (i = 1; i < N; i++)
  73   {
  74     double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
  75     double denom = h[i-1]*h[i];
  76     alpha[i] = numer/denom;
  77   }
  78
  79   double* ell = new double[N+1];
  80   double* mu = new double[N];
  81   double* z = new double[N+1];
  82   double recip;
  83
  84   ell[0] = 1.0;
  85   mu[0] = 0.0;
  86   z[0] = 0.0;
  87
  88   for (i = 1; i < N; i++)
  89   {
  90     ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
  91     recip = 1.0/ell[i];
  92     mu[i] = recip*h[i];
  93     z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
  94   }
  95   ell[N] = 1.0;
  96   z[N] = 0.0;
  97
  98   b = new double[N];
  99   c = new double[N+1];
 100   d = new double[N];
 101
 102   c[N] = 0.0;
 103
 104   for (i = N-1; i >= 0; i--)
 105   {
 106     c[i] = z[i]-mu[i]*c[i+1];
 107     recip = 1.0/h[i];
 108     b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
 109     d[i] = oneThird*recip*(c[i+1]-c[i]);
 110   }
 111
 112   delete[] h;
 113   delete[] hdiff;
 114   delete[] alpha;
 115   delete[] ell;
 116   delete[] mu;
 117   delete[] z;
 118 }
 119
 120 void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
 121     double*& d)
 122 {
 123   double* h = new double[N];
 124   int i;
 125   for (i = 0; i < N; i++)
 126     h[i] = x[i+1]-x[i];
 127
 128   mgcLinearSystemD sys;
 129   double** mat = mgcLinearSystemD::NewMatrix(N+1);  // guaranteed to be zeroed memory
 130   c = mgcLinearSystemD::NewVector(N+1);   // guaranteed to be zeroed memory
 131
 132   // c[0] - c[N] = 0
 133   mat[0][0] = +1.0f;
 134   mat[0][N] = -1.0f;
 135
 136   // h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
 137   //   3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
 138   for (i = 1; i <= N-1; i++)
 139   {
 140     mat[i][i-1] = h[i-1];
 141     mat[i][i  ] = 2.0f*(h[i-1]+h[i]);
 142     mat[i][i+1] = h[i];
 143     c[i] = 3.0f*((a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]);
 144   }
 145
 146   // "wrap around equation" for periodicity
 147   // h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
 148   //   3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
 149   mat[N][N-1] = h[N-1];
 150   mat[N][0  ] = 2.0f*(h[N-1]+h[0]);
 151   mat[N][1  ] = h[0];
 152   c[N] = 3.0f*((a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]);
 153
 154   // solve for c[0] through c[N]
 155   mgcLinearSystemD::Solve(N+1,mat,c);
 156
 157   const double oneThird = 1.0/3.0;
 158   b = new double[N];
 159   d = new double[N];
 160   for (i = 0; i < N; i++)
 161   {
 162     b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
 163     d[i] = oneThird*(c[i+1]-c[i])/h[i];
 164   }
 165
 166   delete[] h;
 167   mgcLinearSystemD::DeleteMatrix(N+1,mat);
 168 }
 169
 170 /* vim:set shiftwidth=4 softtabstop=4 expandtab: */