implement bicubic interpolation

3 files changed, 56 insertions, 2 deletions

diff --git a/src/GipfelWidget.H b/src/GipfelWidget.H
index 8cb6d3f..ed251ef 100644
--- a/src/GipfelWidget.H
+++ b/src/GipfelWidget.H
@@ -44,13 +44,17 @@ class GipfelWidget : public Fl_Widget {
 		static int get_pixel_bilinear(Fl_Image *img, double x, double y,
 			int *r, int *g, int *b);
 
+		static int get_pixel_bicubic(Fl_Image *img, double x, double y,
+			int *r, int *g, int *b);
+
 		static int get_pixel(Fl_Image *img, int x, int y,
 			int *r, int *g, int *b);
 
 	public:
 		typedef enum {
 			NEAREST  = 0,
-			BILINEAR = 1
+			BILINEAR = 1,
+			BICUBIC  = 2
 		} sample_mode_t;
 	
 		GipfelWidget(int X,int Y,int W, int H);
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx
index 1eb8866..2daea2a 100644
--- a/src/GipfelWidget.cxx
+++ b/src/GipfelWidget.cxx
@@ -686,6 +686,9 @@ GipfelWidget::get_pixel(GipfelWidget::sample_mode_t m,
 	if (m == GipfelWidget::BILINEAR) {
 		ret = get_pixel_bilinear(img, px + ((double) img->w()) / 2.0,
 			py + ((double) img->h()) / 2.0, r, g, b);
+	} else if (m == GipfelWidget::BICUBIC) {
+		ret = get_pixel_bicubic(img, px + ((double) img->w()) / 2.0,
+			py + ((double) img->h()) / 2.0, r, g, b);
 	} else {
 		ret = get_pixel_nearest(img, px + ((double) img->w()) / 2.0,
 			py + ((double) img->h()) / 2.0, r, g, b);
@@ -751,6 +754,53 @@ GipfelWidget::get_pixel_bilinear(Fl_Image *img, double x, double y,
 	}
 }
 
+static double
+interp_cubic(double x, double *v) {
+	double a0, a1, a2, a3;
+	double x2 = x * x;
+
+	a0 = v[3] - v[2] - v[0] + v[1];
+	a1 = v[0] - v[1] - a0;
+	a2 = v[2] - v[0];
+	a3 = v[1];
+
+	return a0 * x * x2 +a1 * x2 + a2 * x + a3;
+}
+
+int
+GipfelWidget::get_pixel_bicubic(Fl_Image *img, double x, double y,
+	int *r, int *g, int *b) {
+
+	int fl_x = (int) rint(x);
+	int fl_y = (int) rint(y);
+	int ic[3];
+	double c[3][4];
+	double c1[3][4];
+
+	for (int iy = 0; iy < 4; iy++) {
+		for (int ix = 0; ix < 4; ix++) {
+			if (get_pixel(img, fl_x + ix - 2, fl_y + iy - 2,
+					&ic[0], &ic[1], &ic[2]) != 0) {
+				return 1;
+			}
+
+			for (int l = 0; l < 3; l++) {
+				c[l][ix] = (double) ic[l];
+			}
+		}
+
+		for (int l = 0; l < 3; l++) {
+			c1[l][iy] = interp_cubic(x - fl_x, c[l]);
+		}
+	}
+
+	*r = (int) rint(interp_cubic(x - fl_x, c1[0]));
+	*g = (int) rint(interp_cubic(x - fl_x, c1[1]));
+	*b = (int) rint(interp_cubic(x - fl_x, c1[2]));
+
+	return 0;
+}
+
 int
 GipfelWidget::get_pixel(Fl_Image *img, int x, int y,
                      int *r, int *g, int *b) {
diff --git a/src/gipfel.cxx b/src/gipfel.cxx
index 410471d..8620641 100644
--- a/src/gipfel.cxx
+++ b/src/gipfel.cxx
@@ -497,7 +497,7 @@ int main(int argc, char** argv) {
 			type |= STITCH_VIGNETTE_CALIB;
 		}
 
-		stitch(bilinear_flag?GipfelWidget::BILINEAR:GipfelWidget::NEAREST,
+		stitch(bilinear_flag?GipfelWidget::BICUBIC:GipfelWidget::NEAREST,
 			stitch_w, stitch_h, stitch_from, stitch_to,
 			type, outpath, my_argc, my_argv);