text: 
  #ifndef BICUBIC_C
  #define BICUBIC_C
  
  
  #include "getpixel.c"
  
  
  static float cubic_interpolation(float v[4], float x)
  {
  	return v[1] + 0.5 * x*(v[2] - v[0]
  			+ x*(2.0*v[0] - 5.0*v[1] + 4.0*v[2] - v[3]
  			+ x*(3.0*(v[1] - v[2]) + v[3] - v[0])));
  }
  
  static float bicubic_interpolation_cell(float p[4][4], float x, float y)
  {
  	float v[4];
  	v[0] = cubic_interpolation(p[0], y);
  	v[1] = cubic_interpolation(p[1], y);
  	v[2] = cubic_interpolation(p[2], y);
  	v[3] = cubic_interpolation(p[3], y);
  	return cubic_interpolation(v, x);
  }
  
  static
  void bicubic_interpolation(float *result,
  		float *img, int w, int h, int pd, float x, float y)
  {
  	x -= 1;
  	y -= 1;
  
  	getsample_operator p = getsample_1;
  
  	int ix = floor(x);
  	int iy = floor(y);
  	for (int l = 0; l < pd; l++) {
  		float c[4][4];
  		for (int j = 0; j < 4; j++)
  			for (int i = 0; i < 4; i++)
  				c[i][j] = p(img, w, h, pd, ix + i, iy + j, l);
  		float r = bicubic_interpolation_cell(c, x - ix, y - iy);
  		result[l] = r;
  	}
  }
  
  static
  void bicubic_interpolation_nans(float *result,
  		float *img, int w, int h, int pd, float x, float y)
  {
  	x -= 1;
  	y -= 1;
  
  	getsample_operator p = getsample_nan;
  
  	int ix = floor(x);
  	int iy = floor(y);
  	for (int l = 0; l < pd; l++) {
  		float c[4][4];
  		for (int j = 0; j < 4; j++)
  			for (int i = 0; i < 4; i++)
  				c[i][j] = p(img, w, h, pd, ix + i, iy + j, l);
  		float r = bicubic_interpolation_cell(c, x - ix, y - iy);
  		result[l] = r;
  	}
  }
  
  
  static
  void bicubic_interpolation_boundary(float *result,
  		float *img, int w, int h, int pd, float x, float y,
  		int boundary)
  {
  	x -= 1;
  	y -= 1;
  
  	getsample_operator p;
  	switch(boundary)
  	{
  	default:
  	case 0: p = getsample_0; break;
  	case 1: p = getsample_1; break;
  	case 2: p = getsample_2; break;
  	case -1: p = getsample_error; break;
  	}
  
  	int ix = floor(x);
  	int iy = floor(y);
  	for (int l = 0; l < pd; l++) {
  		float c[4][4];
  		for (int j = 0; j < 4; j++)
  			for (int i = 0; i < 4; i++)
  				c[i][j] = p(img, w, h, pd, ix + i, iy + j, l);
  		float r = bicubic_interpolation_cell(c, x - ix, y - iy);
  		result[l] = r;
  	}
  }
  
  static
  void bicubic_interpolation_boundary2(float *result,
  		float *img, int w, int h, int pd, float x, float y,
  		getsample_operator p)
  {
  	x -= 1;
  	y -= 1;
  
  	int ix = floor(x);
  	int iy = floor(y);
  	for (int l = 0; l < pd; l++) {
  		float c[4][4];
  		for (int j = 0; j < 4; j++)
  			for (int i = 0; i < 4; i++)
  				c[i][j] = p(img, w, h, pd, ix + i, iy + j, l);
  		float r = bicubic_interpolation_cell(c, x - ix, y - iy);
  		result[l] = r;
  	}
  }
  
  #endif//BICUBIC_C
  