office-gobmx/hwpfilter/source/cspline.cpp
2003-10-15 13:44:37 +00:00

208 lines
6.3 KiB
C++

/*************************************************************************
*
* $RCSfile: cspline.cpp,v $
*
* $Revision: 1.1 $
*
* last change: $Author: dvo $ $Date: 2003-10-15 14:38:47 $
*
* The Contents of this file are made available subject to the terms of
* either of the following licenses
*
* - GNU Lesser General Public License Version 2.1
* - Sun Industry Standards Source License Version 1.1
*
* Sun Microsystems Inc., October, 2000
*
* GNU Lesser General Public License Version 2.1
* =============================================
* Copyright 2001 by Mizi Research Inc.
* Copyright 2003 by Sun Microsystems, Inc.
* 901 San Antonio Road, Palo Alto, CA 94303, USA
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License version 2.1, as published by the Free Software Foundation.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*
*
* Sun Industry Standards Source License Version 1.1
* =================================================
* The contents of this file are subject to the Sun Industry Standards
* Source License Version 1.1 (the "License"); You may not use this file
* except in compliance with the License. You may obtain a copy of the
* License at http://www.openoffice.org/license.html.
*
* Software provided under this License is provided on an "AS IS" basis,
* WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING,
* WITHOUT LIMITATION, WARRANTIES THAT THE SOFTWARE IS FREE OF DEFECTS,
* MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE, OR NON-INFRINGING.
* See the License for the specific provisions governing your rights and
* obligations concerning the Software.
*
* The Initial Developer of the Original Code is: Mizi Research Inc.
*
* Copyright: 2001 by Mizi Research Inc.
* Copyright: 2003 by Sun Microsystems, Inc.
*
* All Rights Reserved.
*
* Contributor(s): _______________________________________
*
*
************************************************************************/
// Natural, Clamped, or Periodic Cubic Splines
//
// Input: A list of N+1 points (x_i,a_i), 0 <= i <= N, which are sampled
// from a function, a_i = f(x_i). The function f is unknown. Boundary
// conditions are
// (1) Natural splines: f"(x_0) = f"(x_N) = 0
// (2) Clamped splines: f'(x_0) and f'(x_N) are user-specified.
// (3) Periodic splines: f(x_0) = f(x_N) [in which case a_N = a_0 is
// required in the input], f'(x_0) = f'(x_N), and f"(x_0) = f"(x_N).
//
// Output: b_i, c_i, d_i, 0 <= i <= N-1, which are coefficients for the cubic
// spline S_i(x) = a_i + b_i(x-x_i) + c_i(x-x_i)^2 + d_i(x-x_i)^3 for
// x_i <= x < x_{i+1}.
//
// The natural and clamped algorithms were implemented from
//
// Numerical Analysis, 3rd edition
// Richard L. Burden and J. Douglas Faires
// Prindle, Weber & Schmidt
// Boston, 1985, pp. 122-124.
//
// The algorithm sets up a tridiagonal linear system of equations in the
// c_i values. This can be solved in O(N) time.
//
// The periodic spline algorithm was implemented from my own derivation. The
// linear system of equations is not tridiagonal. For now I use a standard
// linear solver that does not take advantage of the sparseness of the
// matrix. Therefore for very large N, you may have to worry about memory
// usage.
#include "solver.h"
//-----------------------------------------------------------------------------
void NaturalSpline (int N, double* x, double* a, double*& b, double*& c,
double*& d)
{
const double oneThird = 1.0/3.0;
int i;
double* h = new double[N];
double* hdiff = new double[N];
double* alpha = new double[N];
for (i = 0; i < N; i++){
h[i] = x[i+1]-x[i];
}
for (i = 1; i < N; i++)
hdiff[i] = x[i+1]-x[i-1];
for (i = 1; i < N; i++)
{
double numer = 3.0*(a[i+1]*h[i-1]-a[i]*hdiff[i]+a[i-1]*h[i]);
double denom = h[i-1]*h[i];
alpha[i] = numer/denom;
}
double* ell = new double[N+1];
double* mu = new double[N];
double* z = new double[N+1];
double recip;
ell[0] = 1.0;
mu[0] = 0.0;
z[0] = 0.0;
for (i = 1; i < N; i++)
{
ell[i] = 2.0*hdiff[i]-h[i-1]*mu[i-1];
recip = 1.0/ell[i];
mu[i] = recip*h[i];
z[i] = recip*(alpha[i]-h[i-1]*z[i-1]);
}
ell[N] = 1.0;
z[N] = 0.0;
b = new double[N];
c = new double[N+1];
d = new double[N];
c[N] = 0.0;
for (i = N-1; i >= 0; i--)
{
c[i] = z[i]-mu[i]*c[i+1];
recip = 1.0/h[i];
b[i] = recip*(a[i+1]-a[i])-h[i]*(c[i+1]+2.0*c[i])*oneThird;
d[i] = oneThird*recip*(c[i+1]-c[i]);
}
delete[] h;
delete[] hdiff;
delete[] alpha;
delete[] ell;
delete[] mu;
delete[] z;
}
void PeriodicSpline (int N, double* x, double* a, double*& b, double*& c,
double*& d)
{
double* h = new double[N];
int i;
for (i = 0; i < N; i++)
h[i] = x[i+1]-x[i];
mgcLinearSystemD sys;
double** mat = sys.NewMatrix(N+1); // guaranteed to be zeroed memory
c = sys.NewVector(N+1); // guaranteed to be zeroed memory
// c[0] - c[N] = 0
mat[0][0] = +1.0f;
mat[0][N] = -1.0f;
// h[i-1]*c[i-1]+2*(h[i-1]+h[i])*c[i]+h[i]*c[i+1] =
// 3*{(a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]}
for (i = 1; i <= N-1; i++)
{
mat[i][i-1] = h[i-1];
mat[i][i ] = 2.0f*(h[i-1]+h[i]);
mat[i][i+1] = h[i];
c[i] = 3.0f*((a[i+1]-a[i])/h[i] - (a[i]-a[i-1])/h[i-1]);
}
// "wrap around equation" for periodicity
// h[N-1]*c[N-1]+2*(h[N-1]+h[0])*c[0]+h[0]*c[1] =
// 3*{(a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]}
mat[N][N-1] = h[N-1];
mat[N][0 ] = 2.0f*(h[N-1]+h[0]);
mat[N][1 ] = h[0];
c[N] = 3.0f*((a[1]-a[0])/h[0] - (a[0]-a[N-1])/h[N-1]);
// solve for c[0] through c[N]
sys.Solve(N+1,mat,c);
const double oneThird = 1.0/3.0;
b = new double[N];
d = new double[N];
for (i = 0; i < N; i++)
{
b[i] = (a[i+1]-a[i])/h[i] - oneThird*(c[i+1]+2.0f*c[i])*h[i];
d[i] = oneThird*(c[i+1]-c[i])/h[i];
}
sys.DeleteMatrix(N+1,mat);
}