GitHub-Proxy/Slic3r
1
0
Fork 0
mirror of https://git.mirrors.martin98.com/https://github.com/slic3r/Slic3r.git synced 2025-08-06 00:06:02 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity

Include external BSpline solver

Browse Source
This commit is contained in:
Florens Wasserfall 2017-01-26 08:44:34 +01:00
parent 3409ad8ae9
commit c3f7a226a0
4 changed files with 1750 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

879
xs/src/BSpline/BSpline.cpp Normal file
View File

@ -0,0 +1,879 @@
// -*- mode: c++; c-basic-offset: 4; -*-
/************************************************************************
* Copyright 2009 University Corporation for Atmospheric Research.
* All rights reserved.
*
* Use of this code is subject to the standard BSD license:
*
* http://www.opensource.org/licenses/bsd-license.html
*
* See the COPYRIGHT file in the source distribution for the license text,
* or see this web page:
*
* http://www.eol.ucar.edu/homes/granger/bspline/doc/
*
*************************************************************************/
/**
* @file
*
* This file defines the implementation for the BSpline and BSplineBase
* templates.
**/
#include "BSpline.h"
#include "BandedMatrix.h"
#include <vector>
#include <algorithm>
#include <iterator>
#include <iostream>
#include <iomanip>
#include <map>
#include <assert.h>
/*
* Original BSplineBase.cpp start here
*/
/*
* This class simulates a namespace for private symbols used by this template
* implementation which should not pollute the global namespace.
*/
class my
{
public:
template<class T> static inline
T abs(const T t)
{
return (t < 0) ? -t : t;
}
template<class T> static inline
const T& min(const T& a,
const T& b)
{
return (a < b) ? a : b;
}
template<class T> static inline
const T& max(const T& a,
const T& b)
{
return (a > b) ? a : b;
}
};
//////////////////////////////////////////////////////////////////////
template<class T> class Matrix : public BandedMatrix<T>
{
public:
Matrix &operator +=(const Matrix &B)
{
Matrix &A = *this;
typename Matrix::size_type M = A.num_rows();
typename Matrix::size_type N = A.num_cols();
assert(M==B.num_rows());
assert(N==B.num_cols());
typename Matrix::size_type i, j;
for (i=0; i<M; i++)
for (j=0; j<N; j++)
A[i][j] += B[i][j];
return A;
}
inline Matrix & operator=(const Matrix &b)
{
return Copy(*this, b);
}
inline Matrix & operator=(const T &e)
{
BandedMatrix<T>::operator= (e);
return *this;
}
};
//////////////////////////////////////////////////////////////////////
// Our private state structure, which hides our use of some matrix
// template classes.
template<class T> struct BSplineBaseP
{
typedef Matrix<T> MatrixT;
MatrixT Q; // Holds P+Q and its factorization
std::vector<T> X;
std::vector<T> Nodes;
};
//////////////////////////////////////////////////////////////////////
// This array contains the beta parameter for the boundary condition
// constraints. The boundary condition type--0, 1, or 2--is the first
// index into the array, followed by the index of the endpoints. See the
// Beta() method.
template<class T> const double BSplineBase<T>::BoundaryConditions[3][4] =
{
// 0 1 M-1 M
{
-4,
-1,
-1,
-4 },
{
0,
1,
1,
0 },
{
2,
-1,
-1,
2 } };
//////////////////////////////////////////////////////////////////////
template<class T> inline bool BSplineBase<T>::Debug(int on)
{
static bool debug = false;
if (on >= 0)
debug = (on > 0);
return debug;
}
//////////////////////////////////////////////////////////////////////
template<class T> const double BSplineBase<T>::BS_PI = 3.1415927;
//////////////////////////////////////////////////////////////////////
template<class T> const char * BSplineBase<T>::ImplVersion()
{
return ("$Id: BSpline.cpp 6352 2008-05-05 04:40:39Z martinc $");
}
//////////////////////////////////////////////////////////////////////
template<class T> const char * BSplineBase<T>::IfaceVersion()
{
return (_BSPLINEBASE_IFACE_ID);
}
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
template<class T> BSplineBase<T>::~BSplineBase()
{
delete base;
}
// This is a member-wise copy except for replacing our
// private base structure with the source's, rather than just copying
// the pointer. But we use the compiler's default copy constructor for
// constructing our BSplineBaseP.
template<class T> BSplineBase<T>::BSplineBase(const BSplineBase<T> &bb) :
K(bb.K), BC(bb.BC), OK(bb.OK), base(new BSplineBaseP<T>(*bb.base))
{
xmin = bb.xmin;
xmax = bb.xmax;
alpha = bb.alpha;
waveLength = bb.waveLength;
DX = bb.DX;
M = bb.M;
NX = base->X.size();
}
//////////////////////////////////////////////////////////////////////
template<class T> BSplineBase<T>::BSplineBase(const T *x,
int nx,
double wl,
int bc,
int num_nodes) :
NX(0), K(2), OK(false), base(new BSplineBaseP<T>)
{
setDomain(x, nx, wl, bc, num_nodes);
}
//////////////////////////////////////////////////////////////////////
// Methods
template<class T> bool BSplineBase<T>::setDomain(const T *x,
int nx,
double wl,
int bc,
int num_nodes)
{
if ((nx <= 0) || (x == 0) || (wl< 0) || (bc< 0) || (bc> 2)) {
return false;
}
OK = false;
waveLength = wl;
BC = bc;
// Copy the x array into our storage.
base->X.resize(nx);
std::copy(x, x+nx, base->X.begin());
NX = base->X.size();
// The Setup() method determines the number and size of node intervals.
if (Setup(num_nodes)) {
if (Debug()) {
std::cerr << "Using M node intervals: " << M << " of length DX: "
<< DX << std::endl;
std::cerr << "X min: " << xmin << " ; X max: " << xmax << std::endl;
std::cerr << "Data points per interval: " << (float)NX/(float)M
<< std::endl;
std::cerr << "Nodes per wavelength: " << (float)waveLength
/(float)DX << std::endl;
std::cerr << "Derivative constraint degree: " << K << std::endl;
}
// Now we can calculate alpha and our Q matrix
alpha = Alpha(waveLength);
if (Debug()) {
std::cerr << "Cutoff wavelength: " << waveLength << " ; "
<< "Alpha: " << alpha << std::endl;
std::cerr << "Calculating Q..." << std::endl;
}
calculateQ();
if (Debug() && M < 30) {
std::cerr.fill(' ');
std::cerr.precision(2);
std::cerr.width(5);
std::cerr << base->Q << std::endl;
}
if (Debug())
std::cerr << "Calculating P..." << std::endl;
addP();
if (Debug()) {
std::cerr << "Done." << std::endl;
if (M < 30) {
std::cerr << "Array Q after addition of P." << std::endl;
std::cerr << base->Q;
}
}
// Now perform the LU factorization on Q
if (Debug())
std::cerr << "Beginning LU factoring of P+Q..." << std::endl;
if (!factor()) {
if (Debug())
std::cerr << "Factoring failed." << std::endl;
} else {
if (Debug())
std::cerr << "Done." << std::endl;
OK = true;
}
}
return OK;
}
//////////////////////////////////////////////////////////////////////
/*
* Calculate the alpha parameter given a wavelength.
*/
template<class T> double BSplineBase<T>::Alpha(double wl)
{
// K is the degree of the derivative constraint: 1, 2, or 3
double a = (double) (wl / (2 * BS_PI * DX));
a *= a; // a^2
if (K == 2)
a = a * a; // a^4
else if (K == 3)
a = a * a * a; // a^6
return a;
}
//////////////////////////////////////////////////////////////////////
/*
* Return the correct beta value given the node index. The value depends
* on the node index and the current boundary condition type.
*/
template<class T> inline double BSplineBase<T>::Beta(int m)
{
if (m > 1 && m < M-1)
return 0.0;
if (m >= M-1)
m -= M-3;
assert(0 <= BC && BC <= 2);
assert(0 <= m && m <= 3);
return BoundaryConditions[BC][m];
}
//////////////////////////////////////////////////////////////////////
/*
* Given an array of y data points defined over the domain
* of x data points in this BSplineBase, create a BSpline
* object which contains the smoothed curve for the y array.
*/
template<class T> BSpline<T> * BSplineBase<T>::apply(const T *y)
{
return new BSpline<T> (*this, y);
}
//////////////////////////////////////////////////////////////////////
/*
* Evaluate the closed basis function at node m for value x,
* using the parameters for the current boundary conditions.
*/
template<class T> double BSplineBase<T>::Basis(int m,
T x)
{
double y = 0;
double xm = xmin + (m * DX);
double z = my::abs((double)(x - xm) / (double)DX);
if (z < 2.0) {
z = 2 - z;
y = 0.25 * (z*z*z);
z -= 1.0;
if (z > 0)
y -= (z*z*z);
}
// Boundary conditions, if any, are an additional addend.
if (m == 0 || m == 1)
y += Beta(m) * Basis(-1, x);
else if (m == M-1 || m == M)
y += Beta(m) * Basis(M+1, x);
return y;
}
//////////////////////////////////////////////////////////////////////
/*
* Evaluate the deriviative of the closed basis function at node m for
* value x, using the parameters for the current boundary conditions.
*/
template<class T> double BSplineBase<T>::DBasis(int m,
T x)
{
double dy = 0;
double xm = xmin + (m * DX);
double delta = (double)(x - xm) / (double)DX;
double z = my::abs(delta);
if (z < 2.0) {
z = 2.0 - z;
dy = 0.25 * z * z;
z -= 1.0;
if (z > 0) {
dy -= z * z;
}
dy *= ((delta > 0) ? -1.0 : 1.0) * 3.0 / DX;
}
// Boundary conditions, if any, are an additional addend.
if (m == 0 || m == 1)
dy += Beta(m) * DBasis(-1, x);
else if (m == M-1 || m == M)
dy += Beta(m) * DBasis(M+1, x);
return dy;
}
//////////////////////////////////////////////////////////////////////
template<class T> double BSplineBase<T>::qDelta(int m1,
int m2)
/*
* Return the integral of the product of the basis function derivative
* restricted to the node domain, 0 to M.
*/
{
// These are the products of the Kth derivative of the
// normalized basis functions
// given a distance m nodes apart, qparts[K-1][m], 0 <= m <= 3
// Each column is the integral over each unit domain, -2 to 2
static const double qparts[3][4][4] =
{
{
{
0.11250,
0.63750,
0.63750,
0.11250 },
{
0.00000,
0.13125,
-0.54375,
0.13125 },
{
0.00000,
0.00000,
-0.22500,
-0.22500 },
{
0.00000,
0.00000,
0.00000,
-0.01875 } },
{
{
0.75000,
2.25000,
2.25000,
0.75000 },
{
0.00000,
-1.12500,
-1.12500,
-1.12500 },
{
0.00000,
0.00000,
0.00000,
0.00000 },
{
0.00000,
0.00000,
0.00000,
0.37500 } },
{
{
2.25000,
20.25000,
20.25000,
2.25000 },
{
0.00000,
-6.75000,
-20.25000,
-6.75000 },
{
0.00000,
0.00000,
6.75000,
6.75000 },
{
0.00000,
0.00000,
0.00000,
-2.25000 } } };
if (m1 > m2)
std::swap(m1, m2);
if (m2 - m1 > 3)
return 0.0;
double q = 0;
for (int m = my::max(m1-2, 0); m < my::min(m1+2, M); ++m)
q += qparts[K-1][m2-m1][m-m1+2];
return q * alpha;
}
//////////////////////////////////////////////////////////////////////
template<class T> void BSplineBase<T>::calculateQ()
{
Matrix<T> &Q = base->Q;
Q.setup(M+1, 3);
Q = 0;
if (alpha == 0)
return;
// First fill in the q values without the boundary constraints.
int i;
for (i = 0; i <= M; ++i) {
Q[i][i] = qDelta(i, i);
for (int j = 1; j < 4 && i+j <= M; ++j) {
Q[i][i+j] = Q[i+j][i] = qDelta(i, i+j);
}
}
// Now add the boundary constraints:
// First the upper left corner.
float b1, b2, q;
for (i = 0; i <= 1; ++i) {
b1 = Beta(i);
for (int j = i; j < i+4; ++j) {
b2 = Beta(j);
assert(j-i >= 0 && j - i < 4);
q = 0.0;
if (i+1 < 4)
q += b2*qDelta(-1, i);
if (j+1 < 4)
q += b1*qDelta(-1, j);
q += b1*b2*qDelta(-1, -1);
Q[j][i] = (Q[i][j] += q);
}
}
// Then the lower right
for (i = M-1; i <= M; ++i) {
b1 = Beta(i);
for (int j = i - 3; j <= i; ++j) {
b2 = Beta(j);
q = 0.0;
if (M+1-i < 4)
q += b2*qDelta(i, M+1);
if (M+1-j < 4)
q += b1*qDelta(j, M+1);
q += b1*b2*qDelta(M+1, M+1);
Q[j][i] = (Q[i][j] += q);
}
}
}
//////////////////////////////////////////////////////////////////////
template<class T> void BSplineBase<T>::addP()
{
// Add directly to Q's elements
Matrix<T> &P = base->Q;
std::vector<T> &X = base->X;
// For each data point, sum the product of the nearest, non-zero Basis
// nodes
int mx, m, n, i;
for (i = 0; i < NX; ++i) {
// Which node does this put us in?
T &x = X[i];
mx = (int)((x - xmin) / DX);
// Loop over the upper triangle of nonzero basis functions,
// and add in the products on each side of the diagonal.
for (m = my::max(0, mx-1); m <= my::min(M, mx+2); ++m) {
float pn;
float pm = Basis(m, x);
float sum = pm * pm;
P[m][m] += sum;
for (n = m+1; n <= my::min(M, mx+2); ++n) {
pn = Basis(n, x);
sum = pm * pn;
P[m][n] += sum;
P[n][m] += sum;
}
}
}
}
//////////////////////////////////////////////////////////////////////
template<class T> bool BSplineBase<T>::factor()
{
Matrix<T> &LU = base->Q;
if (LU_factor_banded(LU, 3) != 0) {
if (Debug())
std::cerr << "LU_factor_banded() failed." << std::endl;
return false;
}
if (Debug() && M < 30)
std::cerr << "LU decomposition: " << std::endl << LU << std::endl;
return true;
}
//////////////////////////////////////////////////////////////////////
template<class T> inline double BSplineBase<T>::Ratiod(int &ni,
double &deltax,
double &ratiof)
{
deltax = (xmax - xmin) / ni;
ratiof = waveLength / deltax;
double ratiod = (double) NX / (double) (ni + 1);
return ratiod;
}
//////////////////////////////////////////////////////////////////////
// Setup the number of nodes (and hence deltax) for the given domain and
// cutoff wavelength. According to Ooyama, the derivative constraint
// approximates a lo-pass filter if the cutoff wavelength is about 4*deltax
// or more, but it should at least be 2*deltax. We can increase the number
// of nodes to increase the number of nodes per cutoff wavelength.
// However, to get a reasonable representation of the data, the setup
// enforces at least as many nodes as data points in the domain. (This
// constraint assumes reasonably even distribution of data points, since
// its really the density of data points which matters.)
//
// Return zero if the setup fails, non-zero otherwise.
//
// The algorithm in this routine is mostly taken from the FORTRAN
// implementation by James Franklin, NOAA/HRD.
//
template<class T> bool BSplineBase<T>::Setup(int num_nodes)
{
std::vector<T> &X = base->X;
// Find the min and max of the x domain
xmin = X[0];
xmax = X[0];
int i;
for (i = 1; i < NX; ++i) {
if (X[i] < xmin)
xmin = X[i];
else if (X[i] > xmax)
xmax = X[i];
}
if (Debug())
std::cerr << "Xmax=" << xmax << ", Xmin=" << xmin << std::endl;
// Number of node intervals (number of spline nodes - 1).
int ni;
double deltax;
if (num_nodes >= 2) {
// We've been told explicitly the number of nodes to use.
ni = num_nodes - 1;
if (waveLength == 0) {
waveLength = 1.0;
}
if (Debug())
{
std::cerr << "Num nodes explicitly given as " << num_nodes
<< ", wavelength set to " << waveLength << std::endl;
}
} else if (waveLength == 0) {
// Turn off frequency constraint and just set two node intervals per
// data point.
ni = NX * 2;
waveLength = 1;
if (Debug())
{
std::cerr << "Frequency constraint disabled, using 2 intervals "
<< "per node: " << ni << " intervals, wavelength="
<< waveLength << std::endl;
}
} else if (waveLength > xmax - xmin) {
if (Debug())
{
std::cerr << "Wavelength " << waveLength << " exceeds X span: "
<< xmin << " - " << xmax << std::endl;
}
return (false);
} else {
if (Debug())
{
std::cerr << "Searching for a reasonable number of "
<< "intervals for wavelength " << waveLength
<< " while keeping at least 2 intervals per "
<< "wavelength ..." << std::endl;
}
// Minimum acceptable number of node intervals per cutoff wavelength.
static const double fmin = 2.0;
// Start out at a minimum number of intervals, meaning the maximum
// number of points per interval, then work up to the maximum
// number of intervals for which the intervals per wavelength is
// still adequate. I think the minimum must be more than 2 since
// the basis function is evaluated on multiple nodes.
ni = 5;
double ratiof; // Nodes per wavelength for current deltax
double ratiod; // Points per node interval
// Increase the number of node intervals until we reach the minimum
// number of intervals per cutoff wavelength, but only as long as
// we can maintain at least one point per interval.
do {
if (Ratiod(++ni, deltax, ratiof) < 1.0)
{
if (Debug())
{
std::cerr << "At " << ni << " intervals, fewer than "
<< "one point per interval, and "
<< "intervals per wavelength is "
<< ratiof << "." << std::endl;
}
return false;
}
} while (ratiof < fmin);
// Now increase the number of intervals until we have at least 4
// intervals per cutoff wavelength, but only as long as we can
// maintain at least 2 points per node interval. There's also no
// point to increasing the number of intervals if we already have
// 15 or more nodes per cutoff wavelength.
//
do {
if ((ratiod = Ratiod(++ni, deltax, ratiof)) < 1.0 || ratiof > 15.0) {
--ni;
break;
}
} while (ratiof < 4 || ratiod > 2.0);
if (Debug())
{
std::cerr << "Found " << ni << " intervals, "
<< "length " << deltax << ", "
<< ratiof << " nodes per wavelength " << waveLength
<< ", "
<< ratiod << " data points per interval." << std::endl;
}
}
// Store the calculations in our state
M = ni;
DX = (xmax - xmin) / ni;
return (true);
}
//////////////////////////////////////////////////////////////////////
template<class T> const T * BSplineBase<T>::nodes(int *nn)
{
if (base->Nodes.size() == 0) {
base->Nodes.reserve(M+1);
for (int i = 0; i <= M; ++i) {
base->Nodes.push_back(xmin + (i * DX));
}
}
if (nn)
*nn = base->Nodes.size();
assert(base->Nodes.size() == (unsigned)(M+1));
return &base->Nodes[0];
}
//////////////////////////////////////////////////////////////////////
template<class T> std::ostream &operator<<(std::ostream &out,
const std::vector<T> &c)
{
for (typename std::vector<T>::const_iterator it = c.begin(); it < c.end(); ++it)
out << *it << ", ";
out << std::endl;
return out;
}
/*
* Original BSpline.cpp start here
*/
//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////
// BSpline Class
//////////////////////////////////////////////////////////////////////
template<class T> struct BSplineP {
std::vector<T> spline;
std::vector<T> A;
};
//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
/*
* This BSpline constructor constructs and sets up a new base and
* solves for the spline curve coeffiecients all at once.
*/
template<class T> BSpline<T>::BSpline(const T *x,
int nx,
const T *y,
double wl,
int bc_type,
int num_nodes) :
BSplineBase<T>(x, nx, wl, bc_type, num_nodes), s(new BSplineP<T>) {
solve(y);
}
//////////////////////////////////////////////////////////////////////
/*
* Create a new spline given a BSplineBase.
*/
template<class T> BSpline<T>::BSpline(BSplineBase<T> &bb,
const T *y) :
BSplineBase<T>(bb), s(new BSplineP<T>) {
solve(y);
}
//////////////////////////////////////////////////////////////////////
/*
* (Re)calculate the spline for the given set of y values.
*/
template<class T> bool BSpline<T>::solve(const T *y) {
if (!OK)
return false;
// Any previously calculated curve is now invalid.
s->spline.clear();
OK = false;
// Given an array of data points over x and its precalculated
// P+Q matrix, calculate the b vector and solve for the coefficients.
std::vector<T> &B = s->A;
std::vector<T> &A = s->A;
A.clear();
A.resize(M+1);
if (Debug())
std::cerr << "Solving for B..." << std::endl;
// Find the mean of these data
mean = 0.0;
int i;
for (i = 0; i < NX; ++i) {
mean += y[i];
}
mean = mean / (double)NX;
if (Debug())
std::cerr << "Mean for y: " << mean << std::endl;
int mx, m, j;
for (j = 0; j < NX; ++j) {
// Which node does this put us in?
T &xj = base->X[j];
T yj = y[j] - mean;
mx = (int)((xj - xmin) / DX);
for (m = my::max(0, mx-1); m <= my::min(mx+2, M); ++m) {
B[m] += yj * Basis(m, xj);
}
}
if (Debug() && M < 30) {
std::cerr << "Solution a for (P+Q)a = b" << std::endl;
std::cerr << " b: " << B << std::endl;
}
// Now solve for the A vector in place.
if (LU_solve_banded(base->Q, A, 3) != 0) {
if (Debug())
std::cerr << "LU_solve_banded() failed." << std::endl;
} else {
OK = true;
if (Debug())
std::cerr << "Done." << std::endl;
if (Debug() && M < 30) {
std::cerr << " a: " << A << std::endl;
std::cerr << "LU factor of (P+Q) = " << std::endl << base->Q
<< std::endl;
}
}
return (OK);
}
//////////////////////////////////////////////////////////////////////
template<class T> BSpline<T>::~BSpline() {
delete s;
}
//////////////////////////////////////////////////////////////////////
template<class T> T BSpline<T>::coefficient(int n) {
if (OK)
if (0 <= n && n <= M)
return s->A[n];
return 0;
}
//////////////////////////////////////////////////////////////////////
template<class T> T BSpline<T>::evaluate(T x) {
T y = 0;
if (OK) {
int n = (int)((x - xmin)/DX);
for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) {
y += s->A[i] * Basis(i, x);
}
y += mean;
}
return y;
}
//////////////////////////////////////////////////////////////////////
template<class T> T BSpline<T>::slope(T x) {
T dy = 0;
if (OK) {
int n = (int)((x - xmin)/DX);
for (int i = my::max(0, n-1); i <= my::min(M, n+2); ++i) {
dy += s->A[i] * DBasis(i, x);
}
}
return dy;
}
/// Instantiate BSplineBase
template class BSplineBase<double>;
//template class BSplineBase<float>;
/// Instantiate BSpline
template class BSpline<double>;
//template class BSpline<float>;

452
xs/src/BSpline/BSpline.h Normal file
View File

@ -0,0 +1,452 @@
/* -*- mode: c++; c-basic-offset: 4; -*- */
/************************************************************************
* Copyright 2009 University Corporation for Atmospheric Research.
* All rights reserved.
*
* Use of this code is subject to the standard BSD license:
*
* http://www.opensource.org/licenses/bsd-license.html
*
* See the COPYRIGHT file in the source distribution for the license text,
* or see this web page:
*
* http://www.eol.ucar.edu/homes/granger/bspline/doc/
*
*************************************************************************/
#ifndef BSPLINE_H
#define BSPLINE_H
#include <vector>
/*
* Warning: original BSplineBase.h/cpp merged into BSpline.h/cpp to avoid dependency issues caused by Build::WithXSpp which tries to compile all .cpp files in /src
* Original BSplineBase.h starts here!!!
*/
#ifndef _BSPLINEBASE_IFACE_ID
#define _BSPLINEBASE_IFACE_ID "$Id: BSpline.h 6353 2008-05-05 19:30:48Z martinc $"
#endif
/**
* @file
*
* This file defines the BSpline library interface.
*
*/
template <class T> class BSpline;
/*
* Opaque member structure to hide the matrix implementation.
*/
template <class T> struct BSplineBaseP;
/**
* @class BSplineBase
*
* The base class for a spline object containing the nodes for a given
* domain, cutoff wavelength, and boundary condition.
* To smooth a single curve, the BSpline interface contains a constructor
* which both sets up the domain and solves for the spline. Subsequent
* curves over the same domain can be created by apply()ing them to the
* BSpline object, where apply() is a BSplineBase method. [See apply().]
* New curves can also be smoothed within the same BSpline object by
* calling solve() with the new set of y values. [See BSpline.] A
* BSplineBase can be created on its own, in which case all of the
* computations dependent on the x values, boundary conditions, and cutoff
* wavelength have already been completed.
*
* The solution of the cubic b-spline is divided into two parts. The first
* is the setup of the domain given the x values, boundary conditions, and
* wavelength. The second is the solution of the spline for a set of y
* values corresponding to the x values in the domain. The first part is
* done in the creation of the BSplineBase object (or when calling the
* setDomain method). The second part is done when creating a BSpline
* object (or calling solve() on a BSpline object).
*
* A BSpline object can be created with either one of its constructors, or
* by calling apply() on an existing BSplineBase object. Once a spline has
* been solved, it can be evaluated at any x value. The following example
* creates a spline curve and evaluates it over the domain:
*
@verbatim
vector<float> x;
vector<float> y;
{ ... }
int bc = BSplineBase<float>::BC_ZERO_SECOND;
BSpline<float>::Debug = true;
BSpline<float> spline (x.begin(), x.size(), y.begin(), wl, bc);
if (spline.ok())
{
ostream_iterator<float> of(cout, "\t ");
float xi = spline.Xmin();
float xs = (spline.Xmax() - xi) / 2000.0;
for (; xi <= spline.Xmax(); xi += xs)
{
*of++ = spline.evaluate (xi);
}
}
@endverbatim
*
* In the usual usage, the BSplineBase can compute a reasonable number of
* nodes for the spline, balancing between a few desirable factors. There
* needs to be at least 2 nodes per cutoff wavelength (preferably 4 or
* more) for the derivative constraint to reliably approximate a lo-pass
* filter. There should be at least 1 and preferably about 2 data points
* per node (measured just by their number and not by any check of the
* density of points across the domain). Lastly, of course, the fewer the
* nodes then the faster the computation of the spline. The computation of
* the number of nodes happens in the Setup() method during BSplineBase
* construction and when setDomain() is called. If the setup fails to find
* a desirable number of nodes, then the BSplineBase object will return
* false from ok().
*
* The ok() method returns false when a BSplineBase or BSpline could not
* complete any operation successfully. In particular, as mentioned above,
* ok() will return false if some problem was detected with the domain
* values or if no reasonable number of nodes could be found for the given
* cutoff wavelength. Also, ok() on a BSpline object will return false if
* the matrix equation could not be solved, such as after BSpline
* construction or after a call to apply().
*
* If letting Setup() determine the number of nodes is not acceptable, the
* constructors and setDomain() accept the parameter num_nodes. By
* default, num_nodes is passed as zero, forcing Setup() to calculate the
* number of nodes. However, if num_nodes is passed as 2 or greater, then
* Setup() will bypass its own algorithm and accept the given number of
* nodes instead. Obviously, it's up to the programmer to understand the
* affects of the number of nodes on the representation of the data and on
* the solution (or non-solution) of the spline. Remember to check the
* ok() method to detect when the spline solution has failed.
*
* The interface for the BSplineBase and BSpline templates is defined in
* the header file BSpline.h. The implementation is defined in BSpline.cpp.
* Source files which will instantiate the template should include the
* implementation file and @em not the interface. If the implementation
* for a specific type will be linked from elsewhere, such as a
* static library or Windows DLL, source files should only include the
* interface file. On Windows, applications should link with the import
* library BSpline.lib and make sure BSpline.dll is on the path. The DLL
* contains an implementation for BSpline<float> and BSpline<double>.
* For debugging, an application can include the implementation to get its
* own instantiation.
*
* The algorithm is based on the cubic spline described by Katsuyuki Ooyama
* in Montly Weather Review, Vol 115, October 1987. This implementation
* has benefited from comparisons with a previous FORTRAN implementation by
* James L. Franklin, NOAA/Hurricane Research Division. In particular, the
* algorithm in the Setup() method is based mostly on his implementation
* (VICSETUP). The Setup() method finds a suitable default for the number
* of nodes given a domain and cutoff frequency. This implementation
* adopts most of the same constraints, including a constraint that the
* cutoff wavelength not be greater than the span of the domain values: wl
* < max(x) - min(x). If this is not an acceptable constraint, then use the
* num_nodes parameter to specify the number of nodes explicitly.
*
* The cubic b-spline is formulated as the sum of some multiple of the
* basis function centered at each node in the domain. The number of nodes
* is determined by the desired cutoff wavelength and a desirable number of
* x values per node. The basis function is continuous and differentiable
* up to the second degree. A derivative constraint is included in the
* solution to achieve the effect of a low-pass frequency filter with the
* given cutoff wavelength. The derivative constraint can be disabled by
* specifying a wavelength value of zero, which reduces the analysis to a
* least squares fit to a cubic b-spline. The domain nodes, boundary
* constraints, and wavelength determine a linear system of equations,
* Qa=b, where a is the vector of basis function coefficients at each node.
* The coefficient vector is solved by first LU factoring along the
* diagonally banded matrix Q in BSplineBase. The BSpline object then
* computes the B vector for a set of y values and solves for the
* coefficient vector with the LU matrix. Only the diagonal bands are
* stored in memory and calculated during LU factoring and back
* substitution, and the basis function is evaluated as few times as
* possible in computing the diagonal matrix and B vector.
*
* @author Gary Granger (http://www.eol.ucar.edu/homes/granger)
*
@verbatim
Copyright (c) 1998-2009
University Corporation for Atmospheric Research, UCAR
All rights reserved.
@endverbatim
**/
template <class T>
class BSplineBase
{
public:
// Datum type
typedef T datum_type;
/// Return a string describing the implementation version.
static const char *ImplVersion();
/// Return a string describing the interface version.
static const char *IfaceVersion();
/**
* Call this class method with a value greater than zero to enable
* debug messages, or with zero to disable messages. Calling with
* no arguments returns true if debugging enabled, else false.
*/
static bool Debug (int on = -1);
/**
* Boundary condition types.
*/
enum BoundaryConditionTypes
{
/// Set the endpoints of the spline to zero.
BC_ZERO_ENDPOINTS = 0,
/// Set the first derivative of the spline to zero at the endpoints.
BC_ZERO_FIRST = 1,
/// Set the second derivative to zero.
BC_ZERO_SECOND = 2
};
public:
/**
* Construct a spline domain for the given set of x values, cutoff
* wavelength, and boundary condition type. The parameters are the
* same as for setDomain(). Call ok() to check whether domain
* setup succeeded after construction.
*/
BSplineBase (const T *x, int nx,
double wl, int bc_type = BC_ZERO_SECOND,
int num_nodes = 0);
/// Copy constructor
BSplineBase (const BSplineBase &);
/**
* Change the domain of this base. [If this is part of a BSpline
* object, this method {\em will not} change the existing curve or
* re-apply the smoothing to any set of y values.]
*
* The x values can be in any order, but they must be of sufficient
* density to support the requested cutoff wavelength. The setup of
* the domain may fail because of either inconsistency between the x
* density and the cutoff wavelength, or because the resulting matrix
* could not be factored. If setup fails, the method returns false.
*
* @param x The array of x values in the domain.
* @param nx The number of values in the @p x array.
* @param wl The cutoff wavelength, in the same units as the
* @p x values. A wavelength of zero disables
* the derivative constraint.
* @param bc_type The enumerated boundary condition type. If
* omitted it defaults to BC_ZERO_SECOND.
* @param num_nodes The number of nodes to use for the cubic b-spline.
* If less than 2 a reasonable number will be
* calculated automatically, if possible, taking
* into account the given cutoff wavelength.
*
* @see ok().
*/
bool setDomain (const T *x, int nx, double wl,
int bc_type = BC_ZERO_SECOND,
int num_nodes = 0);
/**
* Create a BSpline smoothed curve for the given set of NX y values.
* The returned object will need to be deleted by the caller.
* @param y The array of y values corresponding to each of the nX()
* x values in the domain.
* @see ok()
*/
BSpline<T> *apply (const T *y);
/**
* Return array of the node coordinates. Returns 0 if not ok(). The
* array of nodes returned by nodes() belongs to the object and should
* not be deleted; it will also be invalid if the object is destroyed.
*/
const T *nodes (int *nnodes);
/**
* Return the number of nodes (one more than the number of intervals).
*/
int nNodes () { return M+1; }
/**
* Number of original x values.
*/
int nX () { return NX; }
/// Minimum x value found.
T Xmin () { return xmin; }
/// Maximum x value found.
T Xmax () { return xmin + (M * DX); }
/**
* Return the Alpha value for a given wavelength. Note that this
* depends on the current node interval length (DX).
*/
double Alpha (double wavelength);
/**
* Return alpha currently in use by this domain.
*/
double Alpha () { return alpha; }
/**
* Return the current state of the object, either ok or not ok.
* Use this method to test for valid state after construction or after
* a call to setDomain(). ok() will return false if either fail, such
* as when an appropriate number of nodes and node interval cannot be
* found for a given wavelength, or when the linear equation for the
* coefficients cannot be solved.
*/
bool ok () { return OK; }
virtual ~BSplineBase();
protected:
typedef BSplineBaseP<T> Base;
// Provided
double waveLength; // Cutoff wavelength (l sub c)
int NX;
int K; // Degree of derivative constraint (currently fixed at 2)
int BC; // Boundary conditions type (0,1,2)
// Derived
T xmax;
T xmin;
int M; // Number of intervals (M+1 nodes)
double DX; // Interval length in same units as X
double alpha;
bool OK;
Base *base; // Hide more complicated state members
// from the public interface.
bool Setup (int num_nodes = 0);
void calculateQ ();
double qDelta (int m1, int m2);
double Beta (int m);
void addP ();
bool factor ();
double Basis (int m, T x);
double DBasis (int m, T x);
static const double BoundaryConditions[3][4];
static const double BS_PI;
double Ratiod (int&, double &, double &);
};
/*
* Original BSpline.h start here
*/
template <class T> struct BSplineP;
/**
* Used to evaluate a BSpline.
* Inherits the BSplineBase domain information and interface and adds
* smoothing. See the BSplineBase documentation for a summary of the
* BSpline interface.
*/
template <class T>
class BSpline : public BSplineBase<T>
{
public:
/**
* Create a single spline with the parameters required to set up
* the domain and subsequently smooth the given set of y values.
* The y values must correspond to each of the values in the x array.
* If either the domain setup fails or the spline cannot be solved,
* the state will be set to not ok.
*
* @see ok().
*
* @param x The array of x values in the domain.
* @param nx The number of values in the @p x array.
* @param y The array of y values corresponding to each of the
* nX() x values in the domain.
* @param wl The cutoff wavelength, in the same units as the
* @p x values. A wavelength of zero disables
* the derivative constraint.
* @param bc_type The enumerated boundary condition type. If
* omitted it defaults to BC_ZERO_SECOND.
* @param num_nodes The number of nodes to use for the cubic b-spline.
* If less than 2 a "reasonable" number will be
* calculated automatically, taking into account
* the given cutoff wavelength.
*/
BSpline (const T *x, int nx, /* independent variable */
const T *y, /* dependent values @ ea X */
double wl, /* cutoff wavelength */
int bc_type = BSplineBase<T>::BC_ZERO_SECOND,
int num_nodes = 0);
/**
* A BSpline curve can be derived from a separate @p base and a set
* of data points @p y over that base.
*/
BSpline (BSplineBase<T> &base, const T *y);
/**
* Solve the spline curve for a new set of y values. Returns false
* if the solution fails.
*
* @param y The array of y values corresponding to each of the nX()
* x values in the domain.
*/
bool solve (const T *y);
/**
* Return the evaluation of the smoothed curve
* at a particular @p x value. If current state is not ok(), returns 0.
*/
T evaluate (T x);
/**
* Return the first derivative of the spline curve at the given @p x.
* Returns zero if the current state is not ok().
*/
T slope (T x);
/**
* Return the @p n-th basis coefficient, from 0 to M. If the current
* state is not ok(), or @p n is out of range, the method returns zero.
*/
T coefficient (int n);
virtual ~BSpline();
using BSplineBase<T>::Debug;
using BSplineBase<T>::Basis;
using BSplineBase<T>::DBasis;
protected:
using BSplineBase<T>::OK;
using BSplineBase<T>::M;
using BSplineBase<T>::NX;
using BSplineBase<T>::DX;
using BSplineBase<T>::base;
using BSplineBase<T>::xmin;
using BSplineBase<T>::xmax;
// Our hidden state structure
BSplineP<T> *s;
T mean; // Fit without mean and add it in later
};
#endif

375
xs/src/BSpline/BandedMatrix.h Normal file
View File

@ -0,0 +1,375 @@
/* -*- mode: c++; c-basic-offset: 4; -*- */
/************************************************************************
* Copyright 2009 University Corporation for Atmospheric Research.
* All rights reserved.
*
* Use of this code is subject to the standard BSD license:
*
* http://www.opensource.org/licenses/bsd-license.html
*
* See the COPYRIGHT file in the source distribution for the license text,
* or see this web page:
*
* http://www.eol.ucar.edu/homes/granger/bspline/doc/
*
*************************************************************************/
/**
* Template for a diagonally banded matrix.
**/
#ifndef _BANDEDMATRIX_ID
#define _BANDEDMATRIX_ID "$Id$"
#include <vector>
#include <algorithm>
#include <iostream>
template <class T> class BandedMatrixRow;
template <class T> class BandedMatrix
{
public:
typedef unsigned int size_type;
typedef T element_type;
// Create a banded matrix with the same number of bands above and below
// the diagonal.
BandedMatrix (int N_ = 1, int nbands_off_diagonal = 0) : bands(0)
{
if (! setup (N_, nbands_off_diagonal))
setup ();
}
// Create a banded matrix by naming the first and last non-zero bands,
// where the diagonal is at zero, and bands below the diagonal are
// negative, bands above the diagonal are positive.
BandedMatrix (int N_, int first, int last) : bands(0)
{
if (! setup (N_, first, last))
setup ();
}
// Copy constructor
BandedMatrix (const BandedMatrix &b) : bands(0)
{
Copy (*this, b);
}
inline bool setup (int N_ = 1, int noff = 0)
{
return setup (N_, -noff, noff);
}
bool setup (int N_, int first, int last)
{
// Check our limits first and make sure they make sense.
// Don't change anything until we know it will work.
if (first > last || N_ <= 0)
return false;
// Need at least as many N_ as columns and as rows in the bands.
if (N_ < abs(first) || N_ < abs(last))
return false;
top = last;
bot = first;
N = N_;
out_of_bounds = T();
// Finally setup the diagonal vectors
nbands = last - first + 1;
if (bands) delete[] bands;
bands = new std::vector<T>[nbands];
int i;
for (i = 0; i < nbands; ++i)
{
// The length of each array varies with its distance from the
// diagonal
int len = N - (abs(bot + i));
bands[i].clear ();
bands[i].resize (len);
}
return true;
}
BandedMatrix<T> & operator= (const BandedMatrix<T> &b)
{
return Copy (*this, b);
}
BandedMatrix<T> & operator= (const T &e)
{
int i;
for (i = 0; i < nbands; ++i)
{
std::fill_n (bands[i].begin(), bands[i].size(), e);
}
out_of_bounds = e;
return (*this);
}
~BandedMatrix ()
{
if (bands)
delete[] bands;
}
private:
// Return false if coordinates are out of bounds
inline bool check_bounds (int i, int j, int &v, int &e) const
{
v = (j - i) - bot;
e = (i >= j) ? j : i;
return !(v < 0 || v >= nbands ||
e < 0 || (unsigned int)e >= bands[v].size());
}
static BandedMatrix & Copy (BandedMatrix &a, const BandedMatrix &b)
{
if (a.bands) delete[] a.bands;
a.top = b.top;
a.bot = b.bot;
a.N = b.N;
a.out_of_bounds = b.out_of_bounds;
a.nbands = a.top - a.bot + 1;
a.bands = new std::vector<T>[a.nbands];
int i;
for (i = 0; i < a.nbands; ++i)
{
a.bands[i] = b.bands[i];
}
return a;
}
public:
T &element (int i, int j)
{
int v, e;
if (check_bounds(i, j, v, e))
return (bands[v][e]);
else
return out_of_bounds;
}
const T &element (int i, int j) const
{
int v, e;
if (check_bounds(i, j, v, e))
return (bands[v][e]);
else
return out_of_bounds;
}
inline T & operator() (int i, int j)
{
return element (i-1,j-1);
}
inline const T & operator() (int i, int j) const
{
return element (i-1,j-1);
}
size_type num_rows() const { return N; }
size_type num_cols() const { return N; }
const BandedMatrixRow<T> operator[] (int row) const
{
return BandedMatrixRow<T>(*this, row);
}
BandedMatrixRow<T> operator[] (int row)
{
return BandedMatrixRow<T>(*this, row);
}
private:
int top;
int bot;
int nbands;
std::vector<T> *bands;
int N;
T out_of_bounds;
};
template <class T>
std::ostream &operator<< (std::ostream &out, const BandedMatrix<T> &m)
{
unsigned int i, j;
for (i = 0; i < m.num_rows(); ++i)
{
for (j = 0; j < m.num_cols(); ++j)
{
out << m.element (i, j) << " ";
}
out << std::endl;
}
return out;
}
/*
* Helper class for the intermediate in the [][] operation.
*/
template <class T> class BandedMatrixRow
{
public:
BandedMatrixRow (const BandedMatrix<T> &_m, int _row) : bm(_m), i(_row)
{ }
BandedMatrixRow (BandedMatrix<T> &_m, int _row) : bm(_m), i(_row)
{ }
~BandedMatrixRow () {}
typename BandedMatrix<T>::element_type & operator[] (int j)
{
return const_cast<BandedMatrix<T> &>(bm).element (i, j);
}
const typename BandedMatrix<T>::element_type & operator[] (int j) const
{
return bm.element (i, j);
}
private:
const BandedMatrix<T> &bm;
int i;
};
/*
* Vector multiplication
*/
template <class Vector, class Matrix>
Vector operator* (const Matrix &m, const Vector &v)
{
typename Matrix::size_type M = m.num_rows();
typename Matrix::size_type N = m.num_cols();
assert (N <= v.size());
//if (M > v.size())
// return Vector();
Vector r(N);
for (unsigned int i = 0; i < M; ++i)
{
typename Matrix::element_type sum = 0;
for (unsigned int j = 0; j < N; ++j)
{
sum += m[i][j] * v[j];
}
r[i] = sum;
}
return r;
}
/*
* LU factor a diagonally banded matrix using Crout's algorithm, but
* limiting the trailing sub-matrix multiplication to the non-zero
* elements in the diagonal bands. Return nonzero if a problem occurs.
*/
template <class MT>
int LU_factor_banded (MT &A, unsigned int bands)
{
typename MT::size_type M = A.num_rows();
typename MT::size_type N = A.num_cols();
if (M != N)
return 1;
typename MT::size_type i,j,k;
typename MT::element_type sum;
for (j = 1; j <= N; ++j)
{
// Check for zero pivot
if ( A(j,j) == 0 )
return 1;
// Calculate rows above and on diagonal. A(1,j) remains as A(1,j).
for (i = (j > bands) ? j-bands : 1; i <= j; ++i)
{
sum = 0;
for (k = (j > bands) ? j-bands : 1; k < i; ++k)
{
sum += A(i,k)*A(k,j);
}
A(i,j) -= sum;
}
// Calculate rows below the diagonal.
for (i = j+1; (i <= M) && (i <= j+bands); ++i)
{
sum = 0;
for (k = (i > bands) ? i-bands : 1; k < j; ++k)
{
sum += A(i,k)*A(k,j);
}
A(i,j) = (A(i,j) - sum) / A(j,j);
}
}
return 0;
}
/*
* Solving (LU)x = B. First forward substitute to solve for y, Ly = B.
* Then backwards substitute to find x, Ux = y. Return nonzero if a
* problem occurs. Limit the substitution sums to the elements on the
* bands above and below the diagonal.
*/
template <class MT, class Vector>
int LU_solve_banded(const MT &A, Vector &b, unsigned int bands)
{
typename MT::size_type i,j;
typename MT::size_type M = A.num_rows();
typename MT::size_type N = A.num_cols();
typename MT::element_type sum;
if (M != N || M == 0)
return 1;
// Forward substitution to find y. The diagonals of the lower
// triangular matrix are taken to be 1.
for (i = 2; i <= M; ++i)
{
sum = b[i-1];
for (j = (i > bands) ? i-bands : 1; j < i; ++j)
{
sum -= A(i,j)*b[j-1];
}
b[i-1] = sum;
}
// Now for the backward substitution
b[M-1] /= A(M,M);
for (i = M-1; i >= 1; --i)
{
if (A(i,i) == 0) // oops!
return 1;
sum = b[i-1];
for (j = i+1; (j <= N) && (j <= i+bands); ++j)
{
sum -= A(i,j)*b[j-1];
}
b[i-1] = sum / A(i,i);
}
return 0;
}
#endif /* _BANDEDMATRIX_ID */

44
xs/src/BSpline/COPYRIGHT Normal file
View File

@ -0,0 +1,44 @@
Copyright (c) 1998-2009,2015
University Corporation for Atmospheric Research, UCAR
All rights reserved.
This software is licensed with the standard BSD license:
http://www.opensource.org/licenses/bsd-license.html
When citing this software, here is a suggested reference:
This software is written by Gary Granger of the National Center for
Atmospheric Research (NCAR), sponsored by the National Science Foundation
(NSF). The algorithm is based on the cubic spline described by Katsuyuki
Ooyama in Montly Weather Review, Vol 115, October 1987. This
implementation has benefited from comparisons with a previous FORTRAN
implementation by James L. Franklin, NOAA/Hurricane Research Division.
The text of the license is reproduced below:
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the distribution.
* Neither the name of the UCAR nor the names of its contributors may be
used to endorse or promote products derived from this software
without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.