bspline.c

来自「math library from gnu」· C语言 代码 · 共 1,008 行 · 第 1/2 页

/* bspline/bspline.c
 *
 * Copyright (C) 2006, 2007, 2008 Patrick Alken
 * Copyright (C) 2008 Rhys Ulerich
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_bspline.h>

/*
 * This module contains routines related to calculating B-splines.
 * The algorithms used are described in
 *
 * [1] Carl de Boor, "A Practical Guide to Splines", Springer
 *     Verlag, 1978.
 *
 * The bspline_pppack_* internal routines contain code adapted from
 *
 * [2] "PPPACK - Piecewise Polynomial Package",
 *     http://www.netlib.org/pppack/
 *
 */

#include "bspline.h"

/*
gsl_bspline_alloc()
  Allocate space for a bspline workspace. The size of the
workspace is O(5k + nbreak)

Inputs: k      - spline order (cubic = 4)
        nbreak - number of breakpoints

Return: pointer to workspace
*/

gsl_bspline_workspace *
gsl_bspline_alloc (const size_t k, const size_t nbreak)
{
  if (k == 0)
    {
      GSL_ERROR_NULL ("k must be at least 1", GSL_EINVAL);
    }
  else if (nbreak < 2)
    {
      GSL_ERROR_NULL ("nbreak must be at least 2", GSL_EINVAL);
    }
  else
    {
      gsl_bspline_workspace *w;

      w = (gsl_bspline_workspace *) malloc (sizeof (gsl_bspline_workspace));

      if (w == 0)
	{
	  GSL_ERROR_NULL ("failed to allocate space for workspace",
			  GSL_ENOMEM);
	}

      w->k = k;
      w->km1 = k - 1;
      w->nbreak = nbreak;
      w->l = nbreak - 1;
      w->n = w->l + k - 1;

      w->knots = gsl_vector_alloc (w->n + k);
      if (w->knots == 0)
	{
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for knots vector",
			  GSL_ENOMEM);
	}

      w->deltal = gsl_vector_alloc (k);
      if (w->deltal == 0)
	{
	  free (w->knots);
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for deltal vector",
			  GSL_ENOMEM);
	}

      w->deltar = gsl_vector_alloc (k);
      if (w->deltar == 0)
	{
	  free (w->deltal);
	  free (w->knots);
	  free (w);
	  GSL_ERROR_NULL ("failed to allocate space for deltar vector",
			  GSL_ENOMEM);
	}


      w->B = gsl_vector_alloc (k);
      if (w->B == 0)
	{
	  free (w->deltar);;
	  free (w->deltal);
	  free (w->knots);
	  free (w);
	  GSL_ERROR_NULL
	    ("failed to allocate space for temporary spline vector",
	     GSL_ENOMEM);
	}

      return w;
    }
}				/* gsl_bspline_alloc() */

/*
gsl_bspline_deriv_alloc()
  Allocate space for a bspline derivative workspace. The size of the
workspace is O(2k^2)

Inputs: k      - spline order (cubic = 4)

Return: pointer to workspace
*/

gsl_bspline_deriv_workspace *
gsl_bspline_deriv_alloc (const size_t k)
{
  if (k == 0)
    {
      GSL_ERROR_NULL ("k must be at least 1", GSL_EINVAL);
    }
  else
    {
      gsl_bspline_deriv_workspace *dw;

      dw =
	(gsl_bspline_deriv_workspace *)
	malloc (sizeof (gsl_bspline_deriv_workspace));

      if (dw == 0)
	{
	  GSL_ERROR_NULL ("failed to allocate space for workspace",
			  GSL_ENOMEM);
	}

      dw->A = gsl_matrix_alloc (k, k);
      if (dw->A == 0)
	{
	  free (dw);
	  GSL_ERROR_NULL
	    ("failed to allocate space for derivative work matrix",
	     GSL_ENOMEM);
	}

      dw->dB = gsl_matrix_alloc (k, k + 1);
      if (dw->dB == 0)
	{
	  free (dw->A);
	  free (dw);
	  GSL_ERROR_NULL
	    ("failed to allocate space for temporary derivative matrix",
	     GSL_ENOMEM);
	}

      dw->k = k;

      return dw;
    }
}				/* gsl_bspline_deriv_alloc() */

/* Return number of coefficients */
size_t
gsl_bspline_ncoeffs (gsl_bspline_workspace * w)
{
  return w->n;
}

/* Return order */
size_t
gsl_bspline_order (gsl_bspline_workspace * w)
{
  return w->k;
}

/* Return number of breakpoints */
size_t
gsl_bspline_nbreak (gsl_bspline_workspace * w)
{
  return w->nbreak;
}

/* Return the location of the i-th breakpoint*/
double
gsl_bspline_breakpoint (size_t i, gsl_bspline_workspace * w)
{
  size_t j = i + w->k - 1;
  return gsl_vector_get (w->knots, j);
}

/*
gsl_bspline_free()
  Free a gsl_bspline_workspace.

Inputs: w - workspace to free

Return: none
*/

void
gsl_bspline_free (gsl_bspline_workspace * w)
{
  gsl_vector_free (w->knots);
  gsl_vector_free (w->deltal);
  gsl_vector_free (w->deltar);
  gsl_vector_free (w->B);
  free (w);
}				/* gsl_bspline_free() */

/*
gsl_bspline_deriv_free()
  Free a gsl_bspline_deriv_workspace.

Inputs: dw - workspace to free

Return: none
*/

void
gsl_bspline_deriv_free (gsl_bspline_deriv_workspace * dw)
{
  gsl_matrix_free (dw->A);
  gsl_matrix_free (dw->dB);
  free (dw);
}				/* gsl_bspline_deriv_free() */

/*
gsl_bspline_knots()
  Compute the knots from the given breakpoints:

   knots(1:k) = breakpts(1)
   knots(k+1:k+l-1) = breakpts(i), i = 2 .. l
   knots(n+1:n+k) = breakpts(l + 1)

where l is the number of polynomial pieces (l = nbreak - 1) and
   n = k + l - 1
(using matlab syntax for the arrays)

The repeated knots at the beginning and end of the interval
correspond to the continuity condition there. See pg. 119
of [1].

Inputs: breakpts - breakpoints
        w        - bspline workspace

Return: success or error
*/

int
gsl_bspline_knots (const gsl_vector * breakpts, gsl_bspline_workspace * w)
{
  if (breakpts->size != w->nbreak)
    {
      GSL_ERROR ("breakpts vector has wrong size", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */

      for (i = 0; i < w->k; i++)
	gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, 0));

      for (i = 1; i < w->l; i++)
	{
	  gsl_vector_set (w->knots, w->k - 1 + i,
			  gsl_vector_get (breakpts, i));
	}

      for (i = w->n; i < w->n + w->k; i++)
	gsl_vector_set (w->knots, i, gsl_vector_get (breakpts, w->l));

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_knots() */

/*
gsl_bspline_knots_uniform()
  Construct uniformly spaced knots on the interval [a,b] using
the previously specified number of breakpoints. 'a' is the position
of the first breakpoint and 'b' is the position of the last
breakpoint.

Inputs: a - left side of interval
        b - right side of interval
        w - bspline workspace

Return: success or error

Notes: 1) w->knots is modified to contain the uniformly spaced
          knots

       2) The knots vector is set up as follows (using octave syntax):

          knots(1:k) = a
          knots(k+1:k+l-1) = a + i*delta, i = 1 .. l - 1
          knots(n+1:n+k) = b
*/

int
gsl_bspline_knots_uniform (const double a, const double b,
			   gsl_bspline_workspace * w)
{
  size_t i;			/* looping */
  double delta;			/* interval spacing */
  double x;

  delta = (b - a) / (double) w->l;

  for (i = 0; i < w->k; i++)
    gsl_vector_set (w->knots, i, a);

  x = a + delta;
  for (i = 0; i < w->l - 1; i++)
    {
      gsl_vector_set (w->knots, w->k + i, x);
      x += delta;
    }

  for (i = w->n; i < w->n + w->k; i++)
    gsl_vector_set (w->knots, i, b);

  return GSL_SUCCESS;
}				/* gsl_bspline_knots_uniform() */

/*
gsl_bspline_eval()
  Evaluate the basis functions B_i(x) for all i. This is
a wrapper function for gsl_bspline_eval_nonzero() which
formats the output in a nice way.

Inputs: x - point for evaluation
        B - (output) where to store B_i(x) values
            the length of this vector is
            n = nbreak + k - 2 = l + k - 1 = w->n
        w - bspline workspace

Return: success or error

Notes: The w->knots vector must be initialized prior to calling
       this function (see gsl_bspline_knots())
*/

int
gsl_bspline_eval (const double x, gsl_vector * B, gsl_bspline_workspace * w)
{
  if (B->size != w->n)
    {
      GSL_ERROR ("vector B not of length n", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */
      size_t istart;		/* first non-zero spline for x */
      size_t iend;		/* last non-zero spline for x, knot for x */
      int error;		/* error handling */

      /* find all non-zero B_i(x) values */
      error = gsl_bspline_eval_nonzero (x, w->B, &istart, &iend, w);
      if (error)
	{
	  return error;
	}

      /* store values in appropriate part of given vector */
      for (i = 0; i < istart; i++)
	gsl_vector_set (B, i, 0.0);

      for (i = istart; i <= iend; i++)
	gsl_vector_set (B, i, gsl_vector_get (w->B, i - istart));

      for (i = iend + 1; i < w->n; i++)
	gsl_vector_set (B, i, 0.0);

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_eval() */

/*
gsl_bspline_eval_nonzero()
  Evaluate all non-zero B-spline functions at point x.
These are the B_i(x) for i in [istart, iend].
Always B_i(x) = 0 for i < istart and for i > iend.

Inputs: x      - point at which to evaluate splines
        Bk     - (output) where to store B-spline values (length k)
        istart - (output) B-spline function index of
                 first non-zero basis for given x
        iend   - (output) B-spline function index of
                 last non-zero basis for given x.
                 This is also the knot index corresponding to x.
        w      - bspline workspace

Return: success or error

Notes: 1) the w->knots vector must be initialized before calling
          this function

       2) On output, B contains

             [B_{istart,k}, B_{istart+1,k},
             ..., B_{iend-1,k}, B_{iend,k}]

          evaluated at the given x.
*/

int
gsl_bspline_eval_nonzero (const double x, gsl_vector * Bk, size_t * istart,
			  size_t * iend, gsl_bspline_workspace * w)
{
  if (Bk->size != w->k)
    {
      GSL_ERROR ("Bk vector length does not match order k", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* spline index */
      size_t j;			/* looping */
      int flag = 0;		/* interval search flag */
      int error = 0;		/* error flag */

      i = bspline_find_interval (x, &flag, w);
      error = bspline_process_interval_for_eval (x, &i, flag, w);
      if (error)
	{
	  return error;
	}

      *istart = i - w->k + 1;
      *iend = i;

      bspline_pppack_bsplvb (w->knots, w->k, 1, x, *iend, &j, w->deltal,
			     w->deltar, Bk);

      return GSL_SUCCESS;
    }
}				/* gsl_bspline_eval_nonzero() */

/*
gsl_bspline_deriv_eval()
  Evaluate d^j/dx^j B_i(x) for all i, 0 <= j <= nderiv.
This is a wrapper function for gsl_bspline_deriv_eval_nonzero()
which formats the output in a nice way.

Inputs: x      - point for evaluation
        nderiv - number of derivatives to compute, inclusive.
        dB     - (output) where to store d^j/dx^j B_i(x)
                 values. the size of this matrix is
                 (n = nbreak + k - 2 = l + k - 1 = w->n)
                 by (nderiv + 1)
        w      - bspline derivative workspace

Return: success or error

Notes: 1) The w->knots vector must be initialized prior to calling
          this function (see gsl_bspline_knots())

       2) based on PPPACK's bsplvd
*/

int
gsl_bspline_deriv_eval (const double x, const size_t nderiv, gsl_matrix * dB,
			gsl_bspline_workspace * w,
			gsl_bspline_deriv_workspace * dw)
{
  if (dB->size1 != w->n)
    {
      GSL_ERROR ("dB matrix first dimension not of length n", GSL_EBADLEN);
    }
  else if (dB->size2 < nderiv + 1)
    {
      GSL_ERROR
	("dB matrix second dimension must be at least length nderiv+1",
	 GSL_EBADLEN);
    }
  else if (dw->k < w->k) 
    {
      GSL_ERROR ("derivative workspace is too small", GSL_EBADLEN);
    }
  else
    {
      size_t i;			/* looping */
      size_t j;			/* looping */
      size_t istart;		/* first non-zero spline for x */
      size_t iend;		/* last non-zero spline for x, knot for x */
      int error;		/* error handling */