summaryrefslogtreecommitdiff
path: root/gsl-1.9/bspline/bspline.c
diff options
context:
space:
mode:
Diffstat (limited to 'gsl-1.9/bspline/bspline.c')
-rw-r--r--gsl-1.9/bspline/bspline.c497
1 files changed, 497 insertions, 0 deletions
diff --git a/gsl-1.9/bspline/bspline.c b/gsl-1.9/bspline/bspline.c
new file mode 100644
index 0000000..847a498
--- /dev/null
+++ b/gsl-1.9/bspline/bspline.c
@@ -0,0 +1,497 @@
+/* bspline/bspline.c
+ *
+ * Copyright (C) 2006 Patrick Alken
+ *
+ * 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 2 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.
+ */
+
+static int bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
+ gsl_bspline_workspace *w);
+
+static inline size_t bspline_find_interval(const double x, int *flag,
+ gsl_bspline_workspace *w);
+
+/*
+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 *)
+ calloc(1, 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)
+ {
+ gsl_bspline_free(w);
+ GSL_ERROR_NULL("failed to allocate space for knots vector", GSL_ENOMEM);
+ }
+
+ w->deltal = gsl_vector_alloc(k);
+ w->deltar = gsl_vector_alloc(k);
+ if (!w->deltal || !w->deltar)
+ {
+ gsl_bspline_free(w);
+ GSL_ERROR_NULL("failed to allocate space for delta vectors", GSL_ENOMEM);
+ }
+
+ w->B = gsl_vector_alloc(k);
+ if (w->B == 0)
+ {
+ gsl_bspline_free(w);
+ GSL_ERROR_NULL("failed to allocate space for temporary spline vector", GSL_ENOMEM);
+ }
+
+ return (w);
+ }
+} /* gsl_bspline_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;
+}
+
+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 bspline workspace
+
+Inputs: w - workspace to free
+
+Return: none
+*/
+
+void
+gsl_bspline_free(gsl_bspline_workspace *w)
+{
+ if (!w)
+ return;
+
+ if (w->knots)
+ gsl_vector_free(w->knots);
+
+ if (w->deltal)
+ gsl_vector_free(w->deltal);
+
+ if (w->deltar)
+ gsl_vector_free(w->deltar);
+
+ if (w->B)
+ gsl_vector_free(w->B);
+
+ free(w);
+} /* gsl_bspline_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
+basically a wrapper function for bspline_eval_all()
+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("B vector length does not match n", GSL_EBADLEN);
+ }
+ else
+ {
+ size_t i; /* looping */
+ size_t idx = 0;
+ size_t start; /* first non-zero spline */
+
+ /* find all non-zero B_i(x) values */
+ bspline_eval_all(x, w->B, &idx, w);
+
+ /* store values in appropriate part of given vector */
+
+ start = idx - w->k + 1;
+ for (i = 0; i < start; ++i)
+ gsl_vector_set(B, i, 0.0);
+
+ for (i = start; i <= idx; ++i)
+ gsl_vector_set(B, i, gsl_vector_get(w->B, i - start));
+
+ for (i = idx + 1; i < w->n; ++i)
+ gsl_vector_set(B, i, 0.0);
+
+ return GSL_SUCCESS;
+ }
+} /* gsl_bspline_eval() */
+
+/****************************************
+ * INTERNAL ROUTINES *
+ ****************************************/
+
+/*
+bspline_eval_all()
+ Evaluate all non-zero B-splines B_i(x) using algorithm (8)
+of chapter X of [1]
+
+The idea is something like this. Given x \in [t_i, t_{i+1}]
+with t_i < t_{i+1} and the t_i are knots, the values of the
+B-splines not automatically zero fit into a triangular
+array as follows:
+ 0
+ 0
+0 B_{i-2,3}
+ B_{i-1,2}
+B_{i,1} B_{i-1,3} ...
+ B_{i,2}
+0 B_{i,3}
+ 0
+ 0
+
+where B_{i,k} is the ith B-spline of order k. The boundary 0s
+indicate that those B-splines not in the table vanish at x.
+
+To compute the non-zero B-splines of a given order k, we use
+Eqs. (4) and (5) of chapter X of [1]:
+
+(4) B_{i,1}(x) = { 1, t_i <= x < t_{i+1}
+ 0, else }
+
+(5) B_{i,k}(x) = (x - t_i)
+ ----------------- B_{i,k-1}(x)
+ (t_{i+k-1} - t_i)
+
+ t_{i+k} - x
+ + ----------------- B_{i+1,k-1}(x)
+ t_{i+k} - t_{i+1}
+
+So (4) gives us the first column of the table and we can use
+the recurrence relation (5) to get the rest of the columns.
+
+Inputs: x - point at which to evaluate splines
+ B - (output) where to store B-spline values (length k)
+ idx - (output) B-spline function index of last output
+ value (B_{idx}(x) is stored in the last slot of 'B')
+ 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 = [B_{i-k+1,k}, B_{i-k+2,k}, ..., B_{i-1,k}, B_{i,k}]
+
+ where 'i' is stored in 'idx' on output
+
+ 3) based on PPPACK bsplvb
+*/
+
+static int
+bspline_eval_all(const double x, gsl_vector *B, size_t *idx,
+ gsl_bspline_workspace *w)
+{
+ if (B->size != w->k)
+ {
+ GSL_ERROR("B vector not of length k", GSL_EBADLEN);
+ }
+ else
+ {
+ size_t i; /* spline index */
+ size_t j; /* looping */
+ size_t ii; /* looping */
+ int flag = 0; /* error flag */
+ double saved;
+ double term;
+
+ i = bspline_find_interval(x, &flag, w);
+
+ if (flag == -1)
+ {
+ GSL_ERROR("x outside of knot interval", GSL_EINVAL);
+ }
+ else if (flag == 1)
+ {
+ if (x <= gsl_vector_get(w->knots, i) + GSL_DBL_EPSILON)
+ {
+ --i;
+ }
+ else
+ {
+ GSL_ERROR("x outside of knot interval", GSL_EINVAL);
+ }
+ }
+
+ *idx = i;
+
+ gsl_vector_set(B, 0, 1.0);
+
+ for (j = 0; j < w->k - 1; ++j)
+ {
+ gsl_vector_set(w->deltar, j,
+ gsl_vector_get(w->knots, i + j + 1) - x);
+ gsl_vector_set(w->deltal, j,
+ x - gsl_vector_get(w->knots, i - j));
+
+ saved = 0.0;
+
+ for (ii = 0; ii <= j; ++ii)
+ {
+ term = gsl_vector_get(B, ii) /
+ (gsl_vector_get(w->deltar, ii) +
+ gsl_vector_get(w->deltal, j - ii));
+
+ gsl_vector_set(B, ii,
+ saved +
+ gsl_vector_get(w->deltar, ii) * term);
+
+ saved = gsl_vector_get(w->deltal, j - ii) * term;
+ }
+
+ gsl_vector_set(B, j + 1, saved);
+ }
+
+ return GSL_SUCCESS;
+ }
+} /* bspline_eval_all() */
+
+/*
+bspline_find_interval()
+ Find knot interval such that
+
+ t_i <= x < t_{i + 1}
+
+where the t_i are knot values.
+
+Inputs: x - x value
+ flag - (output) error flag
+ w - bspline workspace
+
+Return: i (index in w->knots corresponding to left limit of interval)
+
+Notes: The error conditions are reported as follows:
+
+Condition Return value Flag
+--------- ------------ ----
+x < t_0 0 -1
+t_i <= x < t_{i+1} i 0
+t_{n+k-1} <= x l+k-1 +1
+*/
+
+static inline size_t
+bspline_find_interval(const double x, int *flag,
+ gsl_bspline_workspace *w)
+{
+ size_t i;
+
+ if (x < gsl_vector_get(w->knots, 0))
+ {
+ *flag = -1;
+ return 0;
+ }
+
+ /* find i such that t_i <= x < t_{i+1} */
+ for (i = w->k - 1; i < w->k + w->l - 1; ++i)
+ {
+ const double ti = gsl_vector_get(w->knots, i);
+ const double tip1 = gsl_vector_get(w->knots, i + 1);
+
+ if (tip1 < ti)
+ {
+ GSL_ERROR("knots vector is not increasing", GSL_EINVAL);
+ }
+
+ if (ti <= x && x < tip1)
+ break;
+ }
+
+ if (i == w->k + w->l - 1)
+ *flag = 1;
+ else
+ *flag = 0;
+
+ return i;
+} /* bspline_find_interval() */
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-01 01:22:59 +0000