LCOV - code coverage report
Current view: top level - hwpfilter/source - cspline.cxx (source / functions) Hit Total Coverage
Test: commit 10e77ab3ff6f4314137acd6e2702a6e5c1ce1fae Lines: 0 70 0.0 %
Date: 2014-11-03 Functions: 0 2 0.0 %
Legend: Lines: hit not hit

          Line data    Source code
       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           0 : void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
      56             :     double*& d)
      57             : {
      58           0 :   const double oneThird = 1.0/3.0;
      59             : 
      60             :   int i;
      61           0 :   double* h = new double[N];
      62           0 :   double* hdiff = new double[N];
      63           0 :   double* alpha = new double[N];
      64             : 
      65           0 :   for (i = 0; i < N; i++){
      66           0 :     h[i] = x[i+1]-x[i];
      67             :   }
      68             : 
      69           0 :   for (i = 1; i < N; i++)
      70           0 :     hdiff[i] = x[i+1]-x[i-1];
      71             : 
      72           0 :   for (i = 1; i < N; i++)
      73             :   {
      74           0 :     double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
      75           0 :     double denom = h[i-1]*h[i];
      76           0 :     alpha[i] = numer/denom;
      77             :   }
      78             : 
      79           0 :   double* ell = new double[N+1];
      80           0 :   double* mu = new double[N];
      81           0 :   double* z = new double[N+1];
      82             :   double recip;
      83             : 
      84           0 :   ell[0] = 1.0;
      85           0 :   mu[0] = 0.0;
      86           0 :   z[0] = 0.0;
      87             : 
      88           0 :   for (i = 1; i < N; i++)
      89             :   {
      90           0 :     ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
      91           0 :     recip = 1.0/ell[i];
      92           0 :     mu[i] = recip*h[i];
      93           0 :     z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
      94             :   }
      95           0 :   ell[N] = 1.0;
      96           0 :   z[N] = 0.0;
      97             : 
      98           0 :   b = new double[N];
      99           0 :   c = new double[N+1];
     100           0 :   d = new double[N];
     101             : 
     102           0 :   c[N] = 0.0;
     103             : 
     104           0 :   for (i = N-1; i >= 0; i--)
     105             :   {
     106           0 :     c[i] = z[i]-mu[i]*c[i+1];
     107           0 :     recip = 1.0/h[i];
     108           0 :     b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
     109           0 :     d[i] = oneThird*recip*(c[i+1]-c[i]);
     110             :   }
     111             : 
     112           0 :   delete[] h;
     113           0 :   delete[] hdiff;
     114           0 :   delete[] alpha;
     115           0 :   delete[] ell;
     116           0 :   delete[] mu;
     117           0 :   delete[] z;
     118           0 : }
     119             : 
     120           0 : void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
     121             :     double*& d)
     122             : {
     123           0 :   double* h = new double[N];
     124             :   int i;
     125           0 :   for (i = 0; i < N; i++)
     126           0 :     h[i] = x[i+1]-x[i];
     127             : 
     128           0 :   mgcLinearSystemD sys;
     129           0 :   double** mat = sys.NewMatrix(N+1);  // guaranteed to be zeroed memory
     130           0 :   c = sys.NewVector(N+1);   // guaranteed to be zeroed memory
     131             : 
     132             :   // c[0] - c[N] = 0
     133           0 :   mat[0][0] = +1.0f;
     134           0 :   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           0 :   for (i = 1; i <= N-1; i++)
     139             :   {
     140           0 :     mat[i][i-1] = h[i-1];
     141           0 :     mat[i][i  ] = 2.0f*(h[i-1]+h[i]);
     142           0 :     mat[i][i+1] = h[i];
     143           0 :     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           0 :   mat[N][N-1] = h[N-1];
     150           0 :   mat[N][0  ] = 2.0f*(h[N-1]+h[0]);
     151           0 :   mat[N][1  ] = h[0];
     152           0 :   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           0 :   sys.Solve(N+1,mat,c);
     156             : 
     157           0 :   const double oneThird = 1.0/3.0;
     158           0 :   b = new double[N];
     159           0 :   d = new double[N];
     160           0 :   for (i = 0; i < N; i++)
     161             :   {
     162           0 :     b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
     163           0 :     d[i] = oneThird*(c[i+1]-c[i])/h[i];
     164             :   }
     165             : 
     166           0 :   delete[] h;
     167           0 :   sys.DeleteMatrix(N+1,mat);
     168           0 : }
     169             : 
     170             : /* vim:set shiftwidth=4 softtabstop=4 expandtab: */

Generated by: LCOV version 1.10