Intermediate C/C++ Data Interpolation (curve fitting)

  0, 0.0
  10, .1736481776669303
  20, .3420201433256687
  30, 0.5
  40, .6427876096865393
  50, 0.766044443118978
  60, .8660254037844386
  70, .9396926207859083
  80, 0.984807753012208
  90, 1

  table2graph sines.txt -p 91

Here are the 91 evenly spaced data points output.

0.000000, 0.000000
1.000000, 0.017452
2.000000, 0.034899
3.000000, 0.052336
4.000000, 0.069756
5.000000, 0.087155
6.000000, 0.104528
7.000000, 0.121869
8.000000, 0.139173
9.000000, 0.156434
10.000000, 0.173648
11.000000, 0.190809
12.000000, 0.207911
13.000000, 0.224951
14.000000, 0.241921
15.000000, 0.258819
16.000000, 0.275637
17.000000, 0.292371
18.000000, 0.309017
19.000000, 0.325568
20.000000, 0.342020
21.000000, 0.358368
22.000000, 0.374606
23.000000, 0.390730
24.000000, 0.406735
25.000000, 0.422617
26.000000, 0.438370
27.000000, 0.453989
28.000000, 0.469471
29.000000, 0.484809
30.000000, 0.500000
31.000000, 0.515038
32.000000, 0.529919
33.000000, 0.544639
34.000000, 0.559193
35.000000, 0.573577
36.000000, 0.587786
37.000000, 0.601816
38.000000, 0.615663
39.000000, 0.629321
40.000000, 0.642788
41.000000, 0.656057
42.000000, 0.669127
43.000000, 0.681992
44.000000, 0.694651
45.000000, 0.707098
46.000000, 0.719331
47.000000, 0.731345
48.000000, 0.743138
49.000000, 0.754706
50.000000, 0.766044
51.000000, 0.777151
52.000000, 0.788021
53.000000, 0.798651
54.000000, 0.809038
55.000000, 0.819177
56.000000, 0.829065
57.000000, 0.838697
58.000000, 0.848070
59.000000, 0.857181
60.000000, 0.866025
61.000000, 0.874600
62.000000, 0.882905
63.000000, 0.890940
64.000000, 0.898707
65.000000, 0.906205
66.000000, 0.913436
67.000000, 0.920400
68.000000, 0.927097
69.000000, 0.933527
70.000000, 0.939693
71.000000, 0.945591
72.000000, 0.951212
73.000000, 0.956545
74.000000, 0.961577
75.000000, 0.966298
76.000000, 0.970694
77.000000, 0.974756
78.000000, 0.978469
79.000000, 0.981824
80.000000, 0.984808
81.000000, 0.987417
82.000000, 0.989683
83.000000, 0.991642
84.000000, 0.993333
85.000000, 0.994795
86.000000, 0.996066
87.000000, 0.997183
88.000000, 0.998186
89.000000, 0.999112
90.000000, 1.000000

Create subdirs 
  src                   sources for this app
  usr32/include/LQ      headers for the libLQ library
  usr32/lib/            the libLQ library goes here.
  o                     a place for temporary object files

  makefile-creator c++ table2graph

  PREFIX = ./usr32

  INCLUDE = -I$(PREFIX)/include/LQ -I./$(SRCDIR) -I ./usr32/include

  LDFLAGS = -lLQ
  LIB = -L/usr/lib -L$(PREFIX)/lib

 
#include <stdio.h>  // FILE*
#include <string.h> // strcmp()
#include <stdlib.h> // exit()

#include <slist.h>
#include <cspline.h>

void dbg(){}

/////////////////////////////////////////////////////
// protos

void args_print(int argc, char** argv)
{
  for(int i = 0; i < argc; i++)
  {
    printf("%s\n", argv[i]);
  }
}

// returns 0 on success
int write_interpolate(const char* ofileName, 
                      char** data_in, 
                      int nin, 
                      int nout);

// returns number of args handled, -1 if exit after
int parse_args(int argc, char** argv);

// shows usage note and returns errcode
int usage(int errcode);

// rotates last to first in list starting at base
void args_rol(char** base, int indx);


/////////////////////////////////////////////////////
// Vars
int numPoints = 10; // default
const char *ifileName, *ofileName;
int format_g = 0; // default %lf

int main(int argc, char** argv)
{
  dbg();
  int offset = parse_args(argc, argv);
  // args_print(argc, argv); // diags
  if(offset < 0) 
    return 0;
  switch(argc - offset)
  {
    case 3: // both files named
      ifileName = argv[offset + 1];
      ofileName = argv[offset + 2];
      break;
    case 2: // input file
      ifileName = argv[offset +1];
      ofileName = "/dev/stdout";
      break;
    default:
      return usage(1); // error
  }
  char** data_in = slist_new();
  slist_fileRead(&data_in, ifileName);
  
  if(slist_count(&data_in) < 3)
    return usage(1);
  
  write_interpolate(ofileName, data_in, slist_count(&data_in), numPoints);
  slist_delete(&data_in);
  return 0;
}

#define streq(a, b) (strcmp(a, b) == 0)

// returns number of args handled, -1 if exit after
int parse_args(int argc, char** argv)
{
  char** args = argv + 1;
  int nargs = argc - 1;
  int args_offset = 0; // count number of args handled
  
  for(int i = 0; i < nargs; i++)
  {
    if(streq(args[i], "--help"))
      return usage(-1);
    
    if(streq(args[i], "-v") || streq(args[i], "--version"))
       return usage(-1);
   
    if(streq(args[i], "-p"))
    {
      char* bugstr = args[i];
      i++;
      bugstr = args[i];
      if(args[i] == 0)
        exit(usage(1));
      if(!sscanf(args[i], "%d", &numPoints))
        exit(usage(1));
      args_rol(args, i);
      args_rol(args, i);
      args_offset += 2;
    }
    if(streq(args[i], "-g"))
    {
      format_g = 1;
      args_rol(args, i);
      args_offset += 1;
    }
  }
  return args_offset;
}

int usage(int errcode)
{
#include "usage_str.dat"
  printf("%s", usage_str);
  return errcode;
}

// rotates num args left, last to first in list
void args_rol(char** base, int indx)
{
  char* last = base[indx];
  memmove(base+ 1, base, (indx) * sizeof(char*));
  base[0] = last;
}


// returns 0 on success
int load_data(SPLINE_DATA* pts_in, char** table, int ntable)
{
  int chk = 0;
  double x, y;
  for(int i = 0; i < ntable; i++)
  {
    chk = sscanf(table[i], "%lf %lf", &x, &y);
    if(chk != 2)
    {
      chk = sscanf(table[i], "%lf, %lf", &x, &y);
      if(chk != 2)    
        return 1;
    }
    pts_in[i].x = x;
    pts_in[i].y = y;
  }    
  return 0;
}


// returns 0 on success
int write_interpolate(
                      const char* ofileName, 
                      char** table, 
                      int ntable, 
                      int nout)
{
  int err = 0;
  FILE* fp = fopen(ofileName, "w");
  if(!fp)
    return 1;
  
  SPLINE_DATA *pts_in, *pts_out;
  pts_in = (SPLINE_DATA*)malloc(ntable * sizeof(SPLINE_DATA));
  pts_out = (SPLINE_DATA*)malloc(nout * sizeof(SPLINE_DATA));
  
  err = load_data(pts_in, table, ntable);
  err = cspline_interpolate(pts_in, ntable, pts_out, nout);
  
  for(int i = 0; i < nout; i++)
  {
    if(format_g)
      fprintf(fp, "%lg, %lg\n", pts_out[i].x, pts_out[i].y);
    else
      fprintf(fp, "%lf, %lf\n", pts_out[i].x, pts_out[i].y);
  }
  
  free(pts_in);
  free(pts_out);
}

// cspline.h -- cubic spline for interpolation/curve fitting

#ifndef cspline_h
#define cspline_h

/*
  The input X should be monotonic increasing, prescaled relative to 
  0.5 and output will be clamped to 0.0 <= out <= 1.0
  Returns 0 on if successful
*/

typedef struct
{
  double x;
  double y;
}SPLINE_DATA;

// used to set incrementing or to restore descending X for 
// input or output data, which must be a working read/writable copy.
void spline_flipX(SPLINE_DATA* data, int ndata);

// Interpolates data from mid to ending points.  End points are copied
// verbatim so these may be 'guesses' if the data is part of a continuum
// where the actual data beyond the limits is unknown.
int cspline_interpolate(SPLINE_DATA* tab, int ntab,
                       SPLINE_DATA* data_out, int nout);

#endif // cspline_h

/*
  cspline.c
  
  Interpolation/ curve fitting algorithm optimized for mid point data.
  
  (C) 2012, by Rainbow Sally, based on code from XForms by T. C. Zhao
  (center to edges algorithm) and from Maxima (interpol.mac), Released 
  under GPL.
 
  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 at
  http://www.gnu.org/copyleft/gpl.html
*/


#include <malloc.h>
#include "cspline.h"

int cspline_interpolate(SPLINE_DATA *tab, int ntab,
                       SPLINE_DATA* data_out, int nout)
{
  int i, k, jo, ilo, ihi, im;
  double sigma, Un, p, Qn, aux, a, b;
  static int nwork;
  static double *Y2, *U;

  // need 3 points minimum
  if (ntab <= 3)
    return 1;

  if (nwork < ntab)
  {
    if (Y2)
    {
      Y2 = (double*)realloc(Y2, sizeof(*Y2) * ntab);
      U = (double*)realloc(U, sizeof(*U) * ntab);
    }
    else
    {
      Y2 = (double*)malloc(sizeof(*Y2) * ntab);
      U = (double*)malloc(sizeof(*U) * ntab);
    }
    nwork = ntab;
  }

  // controlling the lower boundary condition
  Y2[0] = U[0] = 0.0;

  // decomposition loop of the triangular algorithm
  for (i = 1; i < ntab - 1; i++)
  {
    sigma = (tab[i].x - tab[i - 1].x) / (tab[i + 1].x - tab[i - 1].x);
    p = sigma * Y2[i - 1] + 2.0;
    Y2[i] = (sigma - 1.0) / p;
    U[i] = ((tab[i + 1].y - tab[i].y) / (tab[i + 1].x - tab[i].x)) -
        ((tab[i].y - tab[i - 1].y) / (tab[i].x - tab[i - 1].x));
    U[i] = (6.0 * U[i] / (tab[i + 1].x - tab[i - 1].x) - sigma * U[i - 1]) / p;
  }

  // controlling the upper boundary condition
  Qn = Un = 0.0;
  Y2[ntab - 1] = (Un - Qn * U[ntab - 2]) / (Qn * Y2[ntab - 2] + 1.0);
  
  // backsubstitution loop of the tridiagonal algorithm
  for (k = ntab - 2; k >= 0; k--)
    Y2[k] = Y2[k] * Y2[k + 1] + U[k];

  // setup
  double xrange = tab[ntab - 1].x - tab[0].x;
  double stepsize = xrange / (nout - 1);

  data_out[0].x = tab[0].x;
  data_out[0].y = tab[0].y;

  /* constructing the cubic splines */
  for (jo = 0, i = 1; i < nout; i++)
  {
    data_out[i].x = data_out[0].x + (i * stepsize);

    /* center */
    ilo = jo;   // lo
    ihi = ntab;  // hi
    while ((ihi - ilo) > 1)
    {
      im = (ihi + ilo) / 2;
      if (data_out[i].x > tab[im].x)
        ilo = im;
      else
        ihi = im;
    }

    jo = ilo;

    /* interpolate from center outward */
    aux = tab[ihi].x - tab[ilo].x;
    a = (tab[ihi].x - data_out[i].x) / aux;
    b = (data_out[i].x - tab[ilo].x) / aux;

    data_out[i].y = (a * tab[ilo].y + b * tab[ihi].y +
        ((a * a * a - a) * Y2[ilo] + (b * b * b - b) * Y2[ihi]) *
        (aux * aux) / 6.0);
  }

  data_out[nout - 1] = tab[ntab - 1];
  return 0; // no error
}

// used to reverse or restore data order for data with descending X
void spline_flipX(SPLINE_DATA* data_in, int ndata)
{
  int start = 0;
  int end = ndata -1;
  SPLINE_DATA tmp;
  while(start < end)
  {
    tmp = data_in[start];
    data_in[start] = data_in[end];
    data_in[end] = tmp;
    start++;
    end--;
  }
}

/* usage_str.txt converted with txt2cstr */
const char* usage_str =
    "Usage table2graph inputFile [outputFile] [-p nPointsOut ]\n"
    "\n"
    "  Takes input data from a file and interpolates over regular \n"
    "  intervals suitable for conversion to a graph.\n"
    "\n"
    "  If outputFile is not specified the output is sent to stdout.\n"
    "\n"
    "  Switches:\n"
    "\n"
    "  -p N    Default is to output 10 points, to change the number \n"
    "          of output points, use the -p switch followed by the \n"
    "          number of points required.\n"
    "\n"
    "  -g      Default is to output in %lf format.  To use %lg \n"
    "          invoke the -g switch.\n"
    ;

TODO: 
  impliment and document the version switch.
  output the data directly to a plotting window -- kdialog has a msgbox but not a graphbox, oh, what do do? :)