g2_splines.c

ViennaRNA-1.6.1

/*****************************************************************************
**  Copyright (C) 1998-2001  Ljubomir Milanovic & Horst Wagner
**  This file is part of the g2 library
**
**  This library is free software; you can redistribute it and/or
**  modify it under the terms of the GNU Lesser General Public
**  License as published by the Free Software Foundation; either
**  version 2.1 of the License, or (at your option) any later version.
**
**  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
******************************************************************************/
/*
 * g2_splines.c
 * Tijs Michels
 * [e-mail]
 * 06/16/99
 */

#include <math.h>
#include <stdio.h>
#include "g2.h"
#include "g2_util.h"

static void g2_split(int n, const double *points, double *x, double *y);
static void g2_c_spline(int n, const double *points, int m, double *sxy);
static void g2_c_b_spline(int n, const double *points, int m, double *sxy);
static void g2_c_raspln(int n, const double *points, double tn, double *sxy);
static void g2_c_newton(int n, const double *c1, const double *c2, int o, const double *xv, double *yv);
static void g2_c_para_3(int n, const double *points, double *sxy);
static void g2_c_para_5(int n, const double *points, double *sxy);

void g2_split(int n, const double *points, double *x, double *y)
{
   int i;
   for (i = 0; i < n; i++) {
      x[i] = points[i+i];
      y[i] = points[i+i+1];
   }
}

#define eps 1.e-12

void g2_c_spline(int n, const double *points, int m, double *sxy)

/*
 *	FUNCTIONAL DESCRIPTION:
 *
 *	Compute a curve of m points (sx[j],sy[j])
 *	-- j being a positive integer < m --
 *	passing through the n data points (x[i],y[i])
 *	-- i being a positive integer < n --
 *	supplied by the user.
 *	The procedure to determine sy[j] involves
 *	Young's method of successive over-relaxation.
 *
 *	FORMAL ARGUMENTS:
 *
 *	n			number of data points
 *	points			data points (x[i],y[i])
 *	m			number of interpolated points; m = (n-1)*o+1
 *				for o curve points for every data point
 *	sxy			interpolated points (sx[j],sy[j])
 *
 *	IMPLICIT INPUTS:	NONE
 *	IMPLICIT OUTPUTS:	NONE
 *	SIDE EFFECTS:		NONE
 *
 *	REFERENCES:
 *
 *	1. Ralston and Wilf, Mathematical Methods for Digital Computers,
 *	   Vol. II, John Wiley and Sons, New York 1967, pp. 156-158.
 *	2. Greville, T.N.E., Ed., Proceedings of An Advanced Seminar
 *	   Conducted by the Mathematics Research Center, U.S. Army,
 *	   University of Wisconsin, Madison. October 7-9, 1968. Theory
 *	   and Applications of Spline Functions, Academic Press,
 *	   New York / London 1969, pp. 156-167.
 *
 *	AUTHORS:
 *
 *	Josef Heinen	04/06/88	<[e-mail]>
 *	Tijs Michels	06/16/99	<[e-mail]>
 */

{
   int i, j;
   double *x, *y, *g, *h;
   double k, u, delta_g;

   if (n < 3) {
      fputs("\nERROR calling function \"g2_c_spline\":\n"
	    "number of data points input should be at least three\n", stderr);
      return;
   }
   if ((m-1)%(n-1)) {
      fputs("\nWARNING from function \"g2_c_spline\":\n"
	    "number of curve points output for every data point input "
	    "is not an integer\n", stderr);
   }

   x = (double *) g2_malloc(n*4*sizeof(double));
   y = x + n;
   g = y + n;
   h = g + n; /* for the constant copy of g */
   g2_split(n, points, x, y);

   n--; /* last value index */
   k =  x[0]; /* look up once */
   u = (x[n] - k) / (m - 1); /* calculate step outside loop */
   for (j = 0; j < m; j++)	sxy[j+j] = j * u + k; /* x-coordinates */

   for (i = 1; i < n; i++) {
      g[i] = 2. * ((y[i+1] - y[i]) / (x[i+1] - x[i]) -
		   (y[i] - y[i-1]) / (x[i] - x[i-1]))
	/ (x[i+1] - x[i-1]); /* whereas g[i] will later be changed repeatedly */
      h[i] = 1.5 * g[i];     /* copy h[i] of g[i] will remain constant */
   }

   k = 0.;

   do {
      for (u = 0., i = 1; i < n; i++) {
	 delta_g = .5 * (x[i] - x[i-1]) / (x[i+1] - x[i-1]);
	 delta_g = (h[i] -
		    g[i] -
		    g[i-1] * delta_g -      /* 8. - 4 * sqrt(3.) */
		    g[i+1] * (.5 - delta_g)) * 1.0717967697244907832;
	 g[i] += delta_g;

	 if (fabs(delta_g) > u) u = fabs(delta_g);
      }	/* On loop termination u holds the largest delta_g. */

      if (k == 0.)	k = u * eps;
	/* Only executed once, at the end of pass one. So k preserves
	 * the largest delta_g of pass one, multiplied by eps.
	 */
   } while (u > k);

   m += m, i = 1, j = 0;
   do {
      u = sxy[j++]; /* x-coordinate */

      while (x[i] < u)	i++;

      if (--i > n)	i = n;

      k = (u - x[i]) / (x[i+1] - x[i]); /* calculate outside loop */
      sxy[j++] = y[i] +
	(y[i+1] - y[i]) * k +
	(u - x[i]) * (u - x[i+1]) *
	((2. - k) * g[i] +
	 (1. + k) * g[i+1]) / 6.; /* y-coordinate */
   } while (j < m);
   g2_free(x);
}

void g2_spline(int id, int n, double *points, int o)

/*
 *	FORMAL ARGUMENTS:
 *
 *	id			device id
 *	n			number of data points
 *	points			data points (x[i],y[i])
 *	o			number of interpolated points per data point
 *
 *	Given an array of n data points {x[1], y[1], ... x[n], y[n]} plot a
 *	spline curve on device id with o interpolated points per data point.
 *	So the larger o, the more fluent the curve.
 */

{
   int m;
   double *sxy;

   m = (n-1)*o+1;
   sxy = (double*)g2_malloc(m*2*sizeof(double));

   g2_c_spline(n, points, m, sxy);
   g2_poly_line(id, m, sxy);

   g2_free(sxy);
}

void g2_filled_spline(int id, int n, double *points, int o)

/*
 *	FORMAL ARGUMENTS:
 *
 *	id			device id
 *	n			number of data points
 *	points			data points (x[i],y[i])
 *	o			number of interpolated points per data point
 */

{
   int m;
   double *sxy;

   m = (n-1)*o+1;
   sxy = (double*)g2_malloc((m+1)*2*sizeof(double));

   g2_c_spline(n, points, m, sxy);
   sxy[m+m] = points[n+n-2];
   sxy[m+m+1] = points[1];
   g2_filled_polygon(id, m+1, sxy);
   g2_free(sxy);
}

void g2_c_b_spline(int n, const double *points, int m, double *sxy)

/*
 * g2_c_b_spline takes n input points. It uses parameter t
 * to compute sx(t) and sy(t) respectively
 */

{
   int i, j;
   double *x, *y;
   double t, bl1, bl2, bl3, bl4;
   double interval, xi_3, yi_3, xi, yi;

   if (n < 3) {
      fputs("\nERROR calling function \"g2_c_b_spline\":\n"
	    "number of data points input should be at least three\n", stderr);
      return;
   }
   x = (double *) g2_malloc(n*2*sizeof(double));
   y = x + n;
   g2_split(n, points, x, y);

   m--; /* last value index */
   n--; /* last value index */
   interval = (double)n / m;

   for (m += m, i = 2, j = 0; i <= n+1; i++) {
      if (i == 2) {
	 xi_3 = 2 * x[0] - x[1];
	 yi_3 = 2 * y[0] - y[1];
      } else {
	 xi_3 = x[i-3];
	 yi_3 = y[i-3];
      }
      if (i == n+1) {
	 xi = 2 * x[n] - x[n-1];
	 yi = 2 * y[n] - y[n-1];
      } else {
	 xi = x[i];
	 yi = y[i];
      }

      t = fmod(j * interval, 1.);

      while (t < 1. && j < m) {
	 bl1 = (1. - t);
	 bl2 = t * t;	/* t^2 */
	 bl4 = t * bl2;	/* t^3 */
	 bl3 = bl4 - bl2;

	 bl1 = bl1 * bl1 * bl1;
	 bl2 = 3. * (bl3 - bl2) + 4.;
	 bl3 = 3. * (  t - bl3) + 1.;

	 sxy[j++] = (bl1 * xi_3 + bl2 * x[i-2] + bl3 * x[i-1] + bl4 * xi) / 6.; /* x-coordinate */
	 sxy[j++] = (bl1 * yi_3 + bl2 * y[i-2] + bl3 * y[i-1] + bl4 * yi) / 6.; /* y-coordinate */

	 t += interval;
      }
   }
   sxy[m]   = x[n];
   sxy[m+1] = y[n];
   g2_free(x);
}

void g2_b_spline(int id, int n, double *points, int o)

/*
 *	FORMAL ARGUMENTS:
 *
 *	id			device id
 *	n			number of data points
 *	points			data points (x[i],y[i])
 *	o			number of interpolated points per data point
 */

{
   int m;
   double *sxy;

   m = (n-1)*o+1;
   sxy = (double*)g2_malloc(m*2*sizeof(double));

   g2_c_b_spline(n, points, m, sxy);
   g2_poly_line(id, m, sxy);

   g2_free(sxy);
}

void g2_filled_b_spline(int id, int n, double *points, int o)

/*
 *	FORMAL ARGUMENTS:
 *
 *	id			device id
 *	n			number of data points
 *	points			data points (x[i],y[i])
 *	o			number of interpolated points per data point
 */

{
   int m;
   double *sxy;

   m = (n-1)*o+1;
   sxy = (double*)g2_malloc((m+1)*2*sizeof(double));

   g2_c_b_spline(n, points, m, sxy);
   sxy[m+m] = points[n+n-2];
   sxy[m+m+1] = points[1];
   g2_filled_polygon(id, m+1, sxy);

   g2_free(sxy);
}

/*
 *	FUNCTION g2_c_raspln
 *
 *	FUNCTIONAL DESCRIPTION:
 *
 *	This function draws a piecewise cubic polynomial through
 *	the specified data points. The (n-1) cubic polynomials are
 *	basically parametric cubic Hermite polynomials through the
 *	n specified data points with tangent values at the data
 *	points determined by a weighted average of the slopes of
 *	the secant lines. A tension parameter "tn" is provided to
 *	adjust the length of the tangent vector at the data points.
 *	This allows the "roundness" of the curve to be adjusted.
 *	For further information and references on this technique see:
 *
 *	D. Kochanek and R. Bartels, Interpolating Splines With Local
 *	Tension, Continuity and Bias Control, Computer Graphics,
 *	18(1984)3.
 *
 *	AUTHORS:
 *
 *	Dennis Mikkelson	distributed in GPLOT	Jan 7, 1988	F77
 *	Tijs Michels		[e-mail]	Jun 7, 1999	C
 *
 *	FORMAL ARGUMENTS:
 *
 *	n	number of data points, n > 2
 *	points	double array holding the x and y-coords of the data points
 *	tn	double parameter in [0.0, 2.0]. When tn = 0.0,
 *		the curve through the data points is very rounded.
 *		As tn increases the curve is gradually pulled tighter.
 *		When tn = 2.0, the curve is essentially a polyline
 *		through the given data points.
 *	sxy	double array holding the coords of the spline curve
 *
 *	IMPLICIT INPUTS:	NONE
 *	IMPLICIT OUTPUTS:	NONE
 *	SIDE EFFECTS:		NONE
 */

#define nb 40
/*
 * Number of straight connecting lines of which each polynomial consists.
 * So between one data point and the next, (nb-1) points are placed.
 */

void g2_c_raspln(int n, const double *points, double tn, double *sxy)
{
   int i, j;
   double *x, *y;
   double bias, tnFactor, tangentL1, tangentL2;
   double D1x, D1y, D2x, D2y, t1x, t1y, t2x, t2y;
   double h1[nb+1];	/* Values of the Hermite basis functions */
   double h2[nb+1];	/* at nb+1 evenly spaced points in [0,1] */
   double h3[nb+1];
   double h4[nb+1];

   x = (double *) g2_malloc(n*2*sizeof(double));
   y = x + n;
   g2_split(n, points, x, y);

/*
 * First, store the values of the Hermite basis functions in a table h[ ]
 * so no time is wasted recalculating them
 */
   for (i = 0; i < nb+1; i++) {
      double t, tt, ttt;
      t = (double) i / nb;
      tt  = t * t;
      ttt = t * tt;
      h1[i] =  2. * ttt - 3. * tt + 1.;
      h2[i] = -2. * ttt + 3. * tt;
      h3[i] =       ttt - 2. * tt + t;
      h4[i] =       ttt -      tt;
   }

/*
 * Set local tnFactor based on input parameter tn
 */
   if (tn <= 0.) {
      tnFactor = 2.;
      fputs("g2_c_raspln: Using Tension Factor 0.0: very rounded", stderr);
   }
   else if (tn >= 2.) {
      tnFactor = 0.;
      fputs("g2_c_raspln: Using Tension Factor 2.0: not rounded at all", stderr);
   }
   else			tnFactor = 2. - tn;

   D1x = D1y = 0.; /* first point has no preceding point */
   for (j = 0; j < n - 2; j++) {
      t1x = x[j+1] - x[j];
      t1y = y[j+1] - y[j];
      t2x = x[j+2] - x[j+1];
      t2y = y[j+2] - y[j+1];
      tangentL1 = t1x * t1x + t1y * t1y;
      tangentL2 = t2x * t2x + t2y * t2y;
      if (tangentL1 + tangentL2 == 0) bias = .5;
      else bias = tangentL2 / (tangentL1 + tangentL2);
      D2x = tnFactor * (bias  * t1x + (1 - bias) * t2x);
      D2y = tnFactor * (bias  * t1y + (1 - bias) * t2y);
      for (i = 0; i < nb; i++) {
	sxy[2 * nb * j + i + i] =
	   h1[i] * x[j] + h2[i] * x[j+1] + h3[i] * D1x + h4[i] * D2x;
	sxy[2 * nb * j + i + i + 1] =
	   h1[i] * y[j] + h2[i] * y[j+1] + h3[i] * D1y + h4[i] * D2y;
      }
      D1x = D2x; /* store as preceding point in */
      D1y = D2y; /* the next pass */
   }