LinuxQuestions.org
Help answer threads with 0 replies.
Go Back   LinuxQuestions.org > Blogs > rainbowsally
User Name
Password

Notices

Rate this Entry

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

Posted 03-25-2012 at 02:16 PM by rainbowsally
Updated 04-19-2012 at 02:02 PM by rainbowsally (Note about obsolete makefile-creator and replacement)

[Note: makefile-creator was far too limited and has been replaced by mc2, makefile-creator the 2nd generation. All the sources can be downloaded at once here:
http://www.linuxquestions.org/questi...eration-34648/
There may be some differences in the libLQ functions as well, but they should be easily spotted (need '&' prefix, and readPipe() has been renamed pipeRead() for purposes of alphabetically listing by object first then function.
]

Intermediate C/C++ Table2graph utility to interpolate data and output data suitable for generating graphs.

Example: Table2graph using the cspline_interpolate() function.

Features:
  • A cubic spline interpolator for converting table data to graphical point data
  • Usage of our libLQ static lib.
  • Usage and modification of makefile-creator output files
  • Using a homeless lib in a sandbox.

Requires:
[Requires libLQ (our own lib). See prerequisites and links to other blog entries at the bottom.]

First, if all you want is an interplator, it's here. But don't get in a big hurry to snatch the files and run off. There are a few dependencies that are also posted here and if you can save yourself 20 or 30 hours of coding by starting with working code, why not work this up with the dependencies before you start modifying it to your own needs. Besides, there are a few other things you could learn or perspectives you may gain in the process which could save you much more than the 20 - 30 hours of coding, in the long run.

Second, these are antique codes. You have to shop around for them, and to the shame of us Open Sourcers we are letting the "big guys" in KDE and Gnome cut us off from these resources by allowing them to "protect" our machines from "us". If you can't compile on your own system, you CAN'T be part of Open Source development. Nuff said on that, I suppose. If you know what I'm talking about you know what I mean. Make some noise about it.

Old code, circa 1970s is the brains of my air foil models. Other code circa 1980s is used for the spline interpolation. If you think some Windows app or openFOAM will get you there quicker, go for it... but then you are reading this for a reason, and friends, you will not be disappointed.

Now let's look at the test output before we dive in.

Below 45 degrees it's quite accurate, it's accurate to about 10 millionths at around 45 deg (midway between inputs of 40 and 50 degrees) and gets worse from there up to near 90 degrees where the mid points are off by about a little over 1 thousandth.

A sine function is very non-linear but is a useful test case for the interpolation.

Here's the table data which I put in a file named sines.txt

Code:
  0, 0.0
  10, .1736481776669303
  20, .3420201433256687
  30, 0.5
  40, .6427876096865393
  50, 0.766044443118978
  60, .8660254037844386
  70, .9396926207859083
  80, 0.984807753012208
  90, 1
The input data doesn't need to be in evenly spaced steps but it is here for simplicity's sake.

Then using the commandline to get 91 points out (0 - 90 degrees inclusive) from the 10 data points input from the table above:
Code:
  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
Here's the setup for the build.

Code:
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
Place files named as follows under the subdir named 'src'. At this point they can be empty files or the real things, we just need the file names for generating the Makefile at present.

cspline.h
cspline.cpp
main.cpp
usage_str.dat

If you start with your makefile-creator type
Code:
  makefile-creator c++ table2graph
Then edit the makefile to add near the top
Code:
  PREFIX = ./usr32
And change the following vars further down in the Makefile
Code:
  INCLUDE = -I$(PREFIX)/include/LQ -I./$(SRCDIR) -I ./usr32/include

  LDFLAGS = -lLQ
  LIB = -L/usr/lib -L$(PREFIX)/lib
If you have some other makefile generator, good for you. Have at it. We'll rework the makefile-creator eventually to automate changing these things so that files can be added to the build and the Makefile will track the changes without any need to re-edit anything.

Now here's the code for the application/demo and the sources for the interpolator.

file: main.cpp
purpose: create a spline interpolator for tabular data input using cspline funcs.
Code:
 
#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);
}
file: cspline.h
purpose: header for functions
Code:
// 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
file: cspline.c
purpose: curve fitting, interpolation of table data to data points suitable for graphing
Code:
/*
  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--;
  }
}
file: usage_str.dat
purpose: usage info and easy modification of usage info using a text to string utility
Code:
/* 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"
    ;
Code:
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? :)
Here's the 'prerequisite' stuff previously posted on this blog.

Adds objlist and slist to libLQ
http://www.linuxquestions.org/questi...-part-2-34477/
grab objlist.c, objlist.h, slist.c and slist.h

Take a look at the "after testing" stuff about halfway down the page
and grab the makefile-creator here.
http://www.linuxquestions.org/questi...-part-2-34421/

Adds temp strings to libLQ
http://www.linuxquestions.org/questi...-memory-34482/
get tmpstr.c, tmpstr.h and the makefile data if you need it.
Posted in Uncategorized
Views 2097 Comments 0
« Prev     Main     Next »
Total Comments 0

Comments

 

  



All times are GMT -5. The time now is 12:05 AM.

Main Menu
Advertisement

My LQ
Write for LQ
LinuxQuestions.org is looking for people interested in writing Editorials, Articles, Reviews, and more. If you'd like to contribute content, let us know.
Main Menu
Syndicate
RSS1  Latest Threads
RSS1  LQ News
Twitter: @linuxquestions
identi.ca: @linuxquestions
Facebook: linuxquestions Google+: linuxquestions
Open Source Consulting | Domain Registration