"Letters to the Editor"

Modification to entitled "A Cubic Spline Extrema Algorithm" 
by Mike Courtney (DDJ, April 1996).
mike.courtney@mci.com


/*
  File:    cubic_io.c
  Purpose: Input/output program for the FindCubicExtrema() routine.
  Reads an ascii file of x,y data, calls the routine, and lists the returned extrema.
  Author: Mike J. Courtney
*/


#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include "cubic.h"

main()
{
  /* the following are required for main() */
  unsigned int i;            /* array indices */
  unsigned int num_extr;     /* number of extrema found */
  BYTE len;                  /* input data parsing variable */
  char data_file[80];        /* ascii data file name */
  char str[40], line[80];    /* input data parsing variables */
  FILE *fp;                  /* ascii data file pointer */
  struct point *first_extr;  /* pointer to the first extreme structure */

  /* the following are required for FindCubicExtrema() */
  unsigned int num_pnts;     /* number of x,y pairs */
  float *x_in, *y_in;        /* arrays of x and y input data values */
  struct point *extr;        /* pointer to current extreme structure */

  /* open input data file and determine number of input points */
  printf("\n Enter name of input file: ");
  gets(data_file);
  if((fp = fopen(data_file, "r")) != NULL)
    {
    num_pnts = 0;
    while(!feof(fp))
      {
      fgets(line, 80, fp);
      num_pnts++;
      }

    printf("\n The number of input points was %d.",num_pnts);
    
    /* allocate the input data arrays */
    x_in   = malloc(num_pnts * sizeof(float));
    y_in   = malloc(num_pnts * sizeof(float));
  
    /* read in the each data line and parse into x and y arrays */
    rewind(fp);
    for (i = 0; i < num_pnts; i++)
      {
      /* get a line */
      fgets(line, 80, fp);
      len = strcspn(line,",");
      /* get the x value */
      strncpy(str, line, len);
      /* NULL terminate the x value */
      str[len] = '\0';
      x_in[i] = atof(str);
      /* get the y value */
      strcpy(str, line+len+1);
      y_in[i] = atof(str);
      }
    fclose(fp);

    /* allocate first structure of linked list of output extrema */
    extr = (struct point *)malloc(sizeof(struct point));
  
    /* save the address of the first structure in the list */ 
    first_extr = extr;

    /* call the routine that computes the extrema */
    if (FindCubicExtrema(num_pnts, x_in, y_in, extr) == FAILURE)
      printf("\n\7 No extrema found !\n");
    else
      {
      /* print the linked list of extrema */
      printf("\n\n Relative extrema computed:");
      extr = first_extr;
      num_extr = 0;
      while (extr)
        {
        printf("\n  %d.  y = %f at x = %f", num_extr+1, extr->y, extr->x);
        extr = extr->next;
        num_extr++;
        }
      }
    free(x_in);
    free(y_in);
    
    /* FREE THE LINKED LIST OF EXTREME STRUCTURES */
    do
      {
      /* point to first structure */
      extr = first_extr;
      /* save address of next extreme structure so that when we free the
         current one, we still have a pointer to the next one */
      first_extr = extr->next;  
      free(extr);
      }
    while (first_extr != NULL);
    }
  else
    printf("\n Couldn't open %s !", data_file);
}


/*
  File:    cubic_extrema.c
  Purpose: Contains FindCubicExtrema() and supporting routines that implement
           the cubic spline extrema algorithm. Given a set of x,y data points,
           determine in the cubic spline sense the relative extrema of the 
           function describing the data.
  Author:  Mike J. Courtney
*/

#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <assert.h>
#include "cubic.h"


/* 
   Primary routine that implements the cubic spline extrema algorithm. Calls
   ComputeSecDerivs() to compute the second derivatives, computes the 
   quadratic coefficients, calls FindQuadRoots() to solve the quadratic roots,
   determines if the roots are valid abscissa, and calls ComputeY() to compute
   ordinates.
*/

BOOL FindCubicExtrema (
  unsigned int num_pnts,  /* input - number of x & y points */
  float *x,               /* input - array of x values */
  float *y,               /* input - array of y values */
  struct point *extr)     /* output - singly linked list of extrema */
{
  float a, b, c;    /* coefficients of quadratic equation */
  float x1, x2;     /* roots of quadratic equation to be computed */
  unsigned int i;   /* array index */
  float *sec_deriv; /* second derivatives of each data interval */
  BOOL root_stat;   /* computation status returned by FindRoots() */
  BOOL valid_flag;  /* TRUE if at least one valid root found */
  BOOL first_flag;  /* TRUE if current root is the first one */
  
  /* allocate array for second derivatives */
  sec_deriv = malloc(num_pnts * sizeof(float));
  
  /* compute the second derivatives */
  ComputeSecDerivs(num_pnts, x, y, sec_deriv);
  
  /* initialize extrema flags */
  valid_flag = FALSE;
  first_flag = TRUE;
  
  /* loop through all the input points and find the extrema */
  for (i = 0; i < num_pnts - 1; i++)
    {
    /* compute the quadratic coefficients */
    a = 3.0 * (sec_deriv[i+1] - sec_deriv[i]);
    b = 6.0 * (x[i+1] * sec_deriv[i] - x[i] * sec_deriv[i+1]);
    c = 6.0 * (y[i+1] - y[i]) - (2.0 * x[i+1] * x[i+1] - x[i] * x[i] + 2 *
        x[i] * x[i+1]) * sec_deriv[i];
    c -= (x[i+1] * x[i+1] - 2.0 * x[i] * x[i+1] - 2.0 * x[i] * x[i]) *
        sec_deriv[i+1];
    /* determine the roots of the cubic extrema quadratic equation */
    root_stat = FindQuadRoots(a, b, c, &x1, &x2);
    if (root_stat != FAILURE)
      {
      /* if root x1 was calculated successfully */
      if (root_stat != FAILURE1)
        {
        /* Determine if root is within the interval */
        if ((x1 > x[i]) && (x1 < x[i+1]))
          {
          /* first root (extremum) */
          if (first_flag == TRUE)
            first_flag = FALSE;
          /* beyond first valid root so allocate next extremum structure */
          else
            {
            extr->next = (struct point *)malloc(sizeof(struct point));
            assert(extr->next != NULL);
            extr = extr->next;
            }
          extr->next = NULL;
          extr->x = x1;
          /* compute the corresponding value of y at the extreme x value */
          ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
          valid_flag = TRUE;
          }
        }
     /* if root x2 was calculated successfully */
     if (root_stat != FAILURE2)
        {
        /* Determine if root is within the current interval */
        if ((x2 > x[i]) && (x2 < x[i+1]))
          {
          /* first root (extremum) */
          if (first_flag == TRUE)
            first_flag = FALSE;
          /* beyond first valid root so allocate next extremum structure */
          else
            {
            extr->next = (struct point *)malloc(sizeof(struct point));
            extr = extr->next;
            }
          extr->next = NULL;
          extr->x = x2;
          /* compute the corresponding value of y at the extreme x value
*/
          ComputeY(i, x, extr->x, y, sec_deriv, &extr->y);
          valid_flag = TRUE;
          }
        }
      } /* end of if(root_stat != FAILURE) */
    } /* end of for(i) */
    
  free(sec_deriv);
  if (valid_flag == TRUE)
    return SUCCESS;
  else
    {
    /* set next to NULL just in case it was not set in the loop - this
is
       so that free loop will operate properly upon return */
    extr->next = NULL;
    return FAILURE;
    }
  }


/* 
   Use input x,y data to form tridiagonal matrix and compute second 
   derivatives of function in the cubic spline sense. 
*/

void ComputeSecDerivs (
  unsigned int num_pnts,   /* input  - number of x & y points */
  float *x,                /* input  - x values */
  float *y,                /* input  - y values */
  float *sec_deriv)        /* output - 2nd derivatives of intervals */
{
  unsigned int i;          /* index */
  float        ftemp;      /* temporary float */
  float        *main_diag; /* matrix main diagonal array */
  float        *diag;      /* matrix diagonal array */
  float        *right;     /* equation right hand side values */

  main_diag = malloc((num_pnts - 2) * sizeof(float));
  diag      = malloc((num_pnts - 1) * sizeof(float));
  right     = malloc((num_pnts - 2) * sizeof(float));
  
  /* compute the matrix main and off-diagonal values */
  /* even though the calling program is suppose to have guaranteed that
the
     x values are increasing, assert that neither of the diagonal 
     differences are zero to avoid a divide by zero condition */
  for (i = 1; i < num_pnts - 1; i++)
    {
    main_diag[i-1] = 2.0 * (x[i+1] - x[i-1]);
    assert(main_diag[i-1] > 0);
    }
    
  for (i = 0; i < num_pnts - 1; i++)
    {
    diag[i-1] = x[i+1] - x[i];
    assert(diag[i-1] > 0);
    }
       
  /* compute right hand side of equation */
  for (i = 1; i < num_pnts - 1; i++)
    right[i-1] = 6.0 * ((y[i+1] - y[i]) / diag[i-1]
                    - (y[i] - y[i-1]) / diag[i-2]);
  /* forward eliminate tridiagonal */
  sec_deriv[0] = 0.0;
  sec_deriv[num_pnts - 1] = 0.0;
  for (i = 1; i < num_pnts - 2; i++)
    {
    ftemp = diag[i] / main_diag[i];
    right[i] -= (right[i-1] * ftemp);
    main_diag[i] -= (diag[i-1] * ftemp);
    }
  /* backward substitution to solve for second derivative at each knot
*/
  for (i = num_pnts - 2; i > 0; i--)
    sec_deriv[i] = (right[i-1] - diag[i-1] * sec_deriv[i+1]) /
main_diag[i-1];
  free(main_diag);
  free(diag);
  free(right);

  return;
}


/* 
   Solve for roots x1 and x2 of a quadratic equation of the form 
        a * (x * x) + b * x + c = 0
   using
        d = -0.5 * [b + sgn(b) * sqrt(b*b - 4ac)]
   where
      x1 = d / a  and
      x2 = c / d.
      
   This root-finder yields particulary accurate results when a and/or c 
   are small values.
*/

BOOL FindQuadRoots (
  float a,   /* input - coefficient a of quadratic equation */
  float b,   /* input - coefficient b of quadratic equation */
  float c,   /* input - coefficient c of quadratic equation */
  float *x1, /* output - first root computed */
  float *x2) /* output - second root computed */
{
  float d;         /* intermediate value */
  BOOL  root_stat; /* status of root computations */
  
  d = b * b - 4 * a * c;
  if (d < 0)
    return FAILURE;
  else
    {
    d = (float)sqrt((double)d);
    /* make the result of sqrt the sign of b */
    if (b < 0 )
      d = -d;
    d = -0.5 * (b + d);
    
    /* solve for the roots of the quadratic */
    /* if both root computations will yield divide by zero ... forget it! */
    if ( (a == 0) && (d == 0) )
      return FAILURE;
      
    root_stat = SUCCESS;
    /* compute first root */
    if (a == 0)
      root_stat = FAILURE1;
    else
      *x1 = d / a;
    /* compute second root */
    if (d == 0)
      root_stat = FAILURE2;
    else
      *x2 = c / d;
    return root_stat;
    }
  }


/* 
   Given an abscissa (x) location, computes the corresponding cubic
spline
   ordinate (y) value.
*/

void ComputeY (
  unsigned int i,          /* input - array index */
  float        *x,         /* input - array of x values */ 
  float        x_value,    /* input - x value at which to solve for y */
  float        *y,         /* input - array of y values */
  float        *sec_deriv, /* input - array of second derivatives of
each data interval */
  float        *y_value)   /* output - y extreme value at x */
{
  float A, B, C, D; /* cubic spline coefficients */
  float ftemp;      /* temporary float value */

  /* compute the standard cubic spline coefficients */
  A = (x[i + 1] - x_value) / (x[i + 1] - x[i]);
  B = 1 - A;
  ftemp = (float) pow((double)(x[i + 1] - x[i]), 2.0) / 6.0;
  C = (A * A * A - A) * ftemp;
  D = (B * B * B - B) * ftemp;
  /* compute the ordinate value at the abscissa location */
  *y_value = A * y[i] + B * y[i + 1] + C * sec_deriv[i] + D * 
            sec_deriv[i + 1];
  return;
}


/*
  File:    cubic.h
  Purpose: #defines, typedefs, and prototypes for cubic_io.c &
cubic_extrema.c
  Author:  Mike J. Courtney 
*/


/* status defines */

#define SUCCESS     1
#define FAILURE     0
#define FAILURE1   -1
#define FAILURE2   -2


/* flag defines */

#define TRUE        1
#define FALSE       0


/* typedefs */

typedef int           BOOL;
typedef unsigned char BYTE;


/* structures */

struct point {
  struct point *next;
  float x;
  float y;
};


/* prototypes */

void ComputeSecDerivs (unsigned int, float *, float *, float *);

BOOL FindCubicExtrema (unsigned int, float *,float *, struct point *);

BOOL FindQuadRoots (float, float, float, float *, float *);

void ComputeY (unsigned int, float *, float, float *, float *, float *);


Four point (simple max) data
-1.5,0
-0.5,1.0
0.5,1.0
1.5,0

Trajectory data
0,      0
10606.6,        8996.6
21213.2,        14773.2
31819.8,        17329.8
42426.4,        16666.4
53033.0,        12783.0
63639.6,        5679.6

System response (underdamped) data
0,      0
0.25,   0.208
0.5,    0.701
0.75,   1.401
1,      1.519
1.25,   1.193
1.5,    0.848
1.75,   0.756
2,      0.886
2.25,   1.053
2.5,    1.114
2.75,   1.064
3,      0.985
3.25,   0.949
3.5,    0.963

Two-root data
-0.3,   0.547
-0.1,   0.879
0.1,    1.081
0.3,    1.147
1.1,    1.011
1.3,    1.117
1.5,    1.375
1.7,    1.833


