#include <BSplineBase.h>

Inheritance diagram for BSplineBase< T >:

Public Types
enum	BoundaryConditionTypes { BC_ZERO_ENDPOINTS = 0 , BC_ZERO_FIRST = 1 , BC_ZERO_SECOND = 2 }

typedef T	datum_type

Public Member Functions
	BSplineBase (const T *x, int nx, double wl, int bc_type=BC_ZERO_SECOND, int num_nodes=0)

	BSplineBase (const BSplineBase &bb)
	Copy constructor.

BSplineBase &	operator= (const BSplineBase &right)
	Assignment operator.

bool	setDomain (const T *x, int nx, double wl, int bc_type=BC_ZERO_SECOND, int num_nodes=0)

BSpline< T > *	apply (const T *y)

const T *	nodes (int *nnodes)

int	nNodes ()

int	nX ()

T	Xmin ()
	Minimum x value found.

T	Xmax ()
	Maximum x value found.

double	Alpha (double wavelength)

double	Alpha ()

bool	ok ()

Static Public Member Functions
static const char *	Version ()
	Return a string describing the bspline library version.

static bool	Debug (int on=-1)

Protected Types
typedef BSplineBaseP< T >	Base

Protected Member Functions
bool	Setup (int num_nodes=0)
	Setup number of nodes and deltax for the given domain and cutoff wavelength. More...

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)

double	Ratiod (int &, double &, double &)

Protected Attributes
double	waveLength

int	NX

int	K

int	BC

T	xmax

T	xmin

int	M

double	DX

double	alpha

bool	OK

std::unique_ptr< Base >	base

Static Protected Attributes
static const double	BoundaryConditions [3][4]

static const double	PI = 3.1415927

Detailed Description

template<class T>
class BSplineBase< T >

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:

    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);
    }
    }

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 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)

Copyright (c) 1998-2009
University Corporation for Atmospheric Research, UCAR
All rights reserved.

Member Enumeration Documentation

◆ BoundaryConditionTypes

template<class T >

enum BSplineBase::BoundaryConditionTypes

Boundary condition types.

Enumerator
BC_ZERO_ENDPOINTS	Set the endpoints of the spline to zero.
BC_ZERO_FIRST	Set the first derivative of the spline to zero at the endpoints.
BC_ZERO_SECOND	Set the second derivative to zero.

Constructor & Destructor Documentation

◆ BSplineBase()

template<class T >

BSplineBase< T >::BSplineBase	(	const T *	x,
		int	nx,
		double	wl,
		int	bc_type = `BC_ZERO_SECOND`,
		int	num_nodes = `0`
	)

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.

Member Function Documentation

◆ Alpha() [1/2]

template<class T >

double BSplineBase< T >::Alpha ( )

inline

Return alpha currently in use by this domain.

◆ Alpha() [2/2]

template<class T >

double BSplineBase< T >::Alpha ( double wl )

Return the Alpha value for a given wavelength. Note that this depends on the current node interval length (DX).

Calculate the alpha parameter given a wavelength.

◆ apply()

template<class T >

BSpline< T > * BSplineBase< T >::apply ( const T * y )

Create a BSpline smoothed curve for the given set of NX y values. The returned object will need to be deleted by the caller.

Parameters

y	The array of y values corresponding to each of the nX() x values in the domain.

See also: ok()

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.

◆ Basis()

template<class T >

double BSplineBase< T >::Basis	(	int	m,
		T	x
	)

protected

Evaluate the closed basis function at node m for value x, using the parameters for the current boundary conditions.

◆ Beta()

template<class T >

double BSplineBase< T >::Beta ( int m )

inlineprotected

Return the correct beta value given the node index. The value depends on the node index and the current boundary condition type.

◆ DBasis()

template<class T >

double BSplineBase< T >::DBasis	(	int	m,
		T	x
	)

protected

Evaluate the deriviative of the closed basis function at node m for value x, using the parameters for the current boundary conditions.

◆ Debug()

template<class T >

bool BSplineBase< T >::Debug ( int on = -1 )

inlinestatic

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.

◆ nNodes()

template<class T >

int BSplineBase< T >::nNodes ( )

inline

Return the number of nodes (one more than the number of intervals).

◆ nodes()

template<class T >

const T * BSplineBase< T >::nodes ( int * nnodes )

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.

◆ nX()

template<class T >

int BSplineBase< T >::nX ( )

inline

Number of original x values.

◆ ok()

template<class T >

bool BSplineBase< T >::ok ( )

inline

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.

◆ qDelta()

template<class T >

double BSplineBase< T >::qDelta	(	int	m1,
		int	m2
	)

protected

Return the integral of the product of the basis function derivative restricted to the node domain, 0 to M.

◆ setDomain()

template<class T >

bool BSplineBase< T >::setDomain	(	const T *	x,
		int	nx,
		double	wl,
		int	bc_type = `BC_ZERO_SECOND`,
		int	num_nodes = `0`
	)

Change the domain of this base. [If this is part of a BSpline object, this method {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.

Parameters

x	The array of x values in the domain.
nx	The number of values in the `x` array.
wl	The cutoff wavelength, in the same units as the `x` values. A wavelength of zero disables the derivative constraint.
bc_type	The enumerated boundary condition type. If omitted it defaults to BC_ZERO_SECOND.
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 also: ok().

◆ Setup()

template<class T >

bool BSplineBase< T >::Setup ( int num_nodes = 0 )

protected

Setup number of nodes and 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.

Parameters

num_nodes

Returns: bool

Member Data Documentation

◆ BoundaryConditions

template<class T >

const double BSplineBase< T >::BoundaryConditions

staticprotected

Initial value:

=
{
    
    {  -4,  -1,   -1,  -4 },
    {   0,   1,    1,   0 },
    {   2,  -1,   -1,   2 }
}

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.

The documentation for this class was generated from the following files:

BSpline/BSplineBase.h
BSpline/BSplineBase.cpp

Public Types

Public Member Functions

Static Public Member Functions

Protected Types

Protected Member Functions

Protected Attributes

Static Protected Attributes

Detailed Description

template<class T> class BSplineBase< T >

Member Enumeration Documentation

◆ BoundaryConditionTypes

Constructor & Destructor Documentation

◆ BSplineBase()

Member Function Documentation

◆ Alpha() [1/2]

◆ Alpha() [2/2]

◆ apply()

◆ Basis()

◆ Beta()

◆ DBasis()

◆ Debug()

◆ nNodes()

◆ nodes()

◆ nX()

◆ ok()

◆ qDelta()

◆ setDomain()

◆ Setup()

Member Data Documentation

◆ BoundaryConditions

template<class T>
class BSplineBase< T >