bug-gsl
[Top][All Lists]

## [Bug-gsl] bug in akima interpolation

 From: ABX Subject: [Bug-gsl] bug in akima interpolation Date: Fri, 18 Apr 2003 15:58:01 +0200

```I tried to implement akima interpolation in splines of POV-Ray. First I compared
how it looks with code from other sources. When I linked my application with
modules from GSL library I recived shape with bug in the most left point:

http://news.povray.org/povray.binaries.images/31227/221922/akima_spline.png

but when I used code from http://www.magic-software.com/Interpolation1D.html
I received almost similiar shape except all control points are "smooth"

http://news.povray.org/povray.binaries.images/31227/222064/akima_spline.png

Since both libraries are builded different way (temporary and helpers variables
are different) I not tried to compare details of implementations. But I suppose
the difference is located here:

GSL:

m[-2] = 3.0 * m[0] - 2.0 * m[1];
m[-1] = 2.0 * m[0] - m[1];
m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];

Magic Software:

afSlope[1] = 2.0f*afSlope[2] - afSlope[3];
afSlope[0] = 2.0f*afSlope[1] - afSlope[2];
afSlope[iQuantity+1] = 2.0f*afSlope[iQuantity] - afSlope[iQuantity-1];
afSlope[iQuantity+2] = 2.0f*afSlope[iQuantity+1] - afSlope[iQuantity];

In my example splines were open, thought first and last point are equal. The
images were created with following set of data:

t,<x(t),y(t)>

0,<0,0>
1,<0,2>
2,<2,1>
3,<0,0>
4,<-1,2>
5,<-2,1>
6,<-1,-1>
7,<0,0>

I do not have original code I tested this cases and even then I could not
deliver it since I linked it with not released code of next POV. But I imagine
to recreate this bug with current gsl code you can take gsl/interpolation/test.c
file and replace two functions with something like:

static int
test_akima (void)
{
int s;

double data_t[8] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
double data_x[8] = { 0.0, 0.0, 2.0, 0.0,-1.0,-2.0,-1.0, 0.0 };
double data_y[8] = { 0.0, 2.0, 1.0, 0.0, 2.0, 1.0,-1.0, 0.0 };

xy_table data_table_x = make_xy_table(data_t, data_x, 8);
xy_table data_table_y = make_xy_table(data_t, data_y, 8);

s = test_interp (&data_table_x, &data_table_y, gsl_interp_akima);
gsl_test (s, "akima interpolation");
return s;
}

static int
test_interp (
const xy_table * data_table_x,
const xy_table * data_table_y,
const gsl_interp_type * T
)
{
int status = 0;
size_t i;

gsl_interp_accel *a_x = gsl_interp_accel_alloc ();
gsl_interp_accel *a_y = gsl_interp_accel_alloc ();
gsl_interp *interp_x = gsl_interp_alloc (T, data_table_x->n);
gsl_interp *interp_y = gsl_interp_alloc (T, data_table_y->n);
gsl_interp_init (interp_x, data_table_x->x, data_table_x->y, data_table_x->n);
gsl_interp_init (interp_y, data_table_y->x, data_table_y->y, data_table_y->n);

for (DBL i; i<data_table_x->n; i=i+step)  // step has to be defined earlier
{
// loop to get values
double x;
double y;
gsl_interp_eval_e (interp_x, data_table_x->x, data_table_x->y, i,a_x,&x);
gsl_interp_eval_e (interp_y, data_table_y->x, data_table_y->y, i,a_y,&y);
// output vector f(i)=<x(i),y(i)> in your favorite way
status++;
}

gsl_interp_accel_free (a_x);
gsl_interp_accel_free (a_y);
gsl_interp_free (interp_x);
gsl_interp_free (interp_y);

return status;
}

```