try linear version - does not work yet

3 files changed, 99 insertions, 238 deletions

diff --git a/src/Stitch.H b/src/Stitch.H
index 62175be..29f0ef3 100644
--- a/src/Stitch.H
+++ b/src/Stitch.H
@@ -17,11 +17,13 @@
 class Stitch {
 	private:
 		GipfelWidget *gipf[MAX_PICS];
-		double color_adjust[MAX_PICS][3];
-		double V1, V2;
+		int num_pics;
+		double color_adjust[MAX_PICS][3][3];
 		OutputImage *single_images[MAX_PICS];
 		OutputImage *merged_image;
 
+		uchar color_correct(uchar c, double a, int pic, int color);
+
 	public:
 		Stitch();
 
@@ -36,9 +38,6 @@ class Stitch {
 		int resample(GipfelWidget::sample_mode_t m,
 			int w, int h, double view_start, double view_end);
 
-		int color_calib(GipfelWidget::sample_mode_t m,
-			int w, int h, double view_start, double view_end);
-
 		int vignette_calib(GipfelWidget::sample_mode_t m,
 			int w, int h, double view_start, double view_end);
 };
diff --git a/src/Stitch.cxx b/src/Stitch.cxx
index 1391fff..113dc62 100644
--- a/src/Stitch.cxx
+++ b/src/Stitch.cxx
@@ -28,8 +28,7 @@ Stitch::Stitch() {
 		single_images[i] = NULL;
 	}
 	merged_image = NULL;
-	V1 = 0.213895;
-	V2 = 0.089561;
+	num_pics = 0;
 
 }
 
@@ -52,6 +51,8 @@ Stitch::load_image(char *file) {
 			if (gipf[i]->load_image(file) != 0) {
 				delete gipf[i];
 				gipf[i] = NULL;
+			} else {
+				num_pics++;
 			}
 			break;
 		}
@@ -85,54 +86,11 @@ Stitch::set_output(const char *file, OutputImage *img) {
 	return ret;
 }
 
-typedef struct {
-	int c1[3];
-	int c2[3];
-} sample_t;
-
-static int
-get_color_adjust(const sample_t *samples, int n_samples,
-	const double *color_adjust_1, double *color_adjust_2) {
-
-	gsl_matrix *X, *cov;
-	gsl_vector *yv,  *r;
-	double chisq;
-	
-	X = gsl_matrix_alloc(n_samples, 1);
-	yv = gsl_vector_alloc(n_samples);
-	r = gsl_vector_alloc(1);
-	cov = gsl_matrix_alloc (1, 1);
-
-	for (int c = 0; c<3; c++) {
-		for (int i = 0; i< n_samples; i++) {
-			gsl_matrix_set(X, i, 0, samples[i].c2[c]);
-			gsl_vector_set(yv, i, samples[i].c1[c] * color_adjust_1[c]); 
-		}
-
-		
-        gsl_multifit_linear_workspace * work 
-           = gsl_multifit_linear_alloc (n_samples, 1);
-
-        gsl_multifit_linear (X, yv, r, cov, &chisq, work);
-		gsl_multifit_linear_free (work);
-
-		color_adjust_2[c] = gsl_vector_get(r,0);
-	}
-
-	return 0;
-
-}
-
-
 
 int
-Stitch::color_calib(GipfelWidget::sample_mode_t m,
+Stitch::vignette_calib(GipfelWidget::sample_mode_t m,
 	int w, int h, double view_start, double view_end) {
 
-	sample_t * sample_points[MAX_PICS][MAX_PICS];
-
-	memset(sample_points, 0, sizeof(sample_points));
-
 	view_start = view_start * deg2rad;
 	view_end = view_end * deg2rad;
 
@@ -142,18 +100,20 @@ Stitch::color_calib(GipfelWidget::sample_mode_t m,
 	int merged_pixel_set;
 	double radius = (double) w / (view_end -view_start);
 
-	if (merged_image) {
-		merged_image->init(w, h);
-	}
+	int max_samples = 10 * 3, n_samples = 0;
+	int ret;
+	int n_vars = (num_pics-1) * 9;
+	gsl_matrix *X, *cov;
+	gsl_vector *yv,  *c;
+	double chisq;
 
-	for(int j=0; j<3; j++) {
-		color_adjust[0][j] = 1.0;
-	}
+	X = gsl_matrix_alloc(max_samples, n_vars);
+	yv = gsl_vector_alloc(max_samples);
+	c = gsl_vector_alloc(n_vars);
+	cov = gsl_matrix_alloc (n_vars, n_vars);
 
-	for (int i=1; i<MAX_PICS; i++) {
-		for(int j=0; j<3; j++) {
-			color_adjust[i][j] = 0.0;
-		}
+	if (merged_image) {
+		merged_image->init(w, h);
 	}
 
 	for (int y=0; y<h; y++) {
@@ -163,72 +123,78 @@ Stitch::color_calib(GipfelWidget::sample_mode_t m,
 			double a_view;
 			a_view = view_start + x * step_view;
 			merged_pixel_set = 0;
+			double a1, a2;
 			uchar c1[3], c2[3];
+			double c1d[3], c2d[3];
 
-			for (int p1=0; p1<MAX_PICS; p1++) {
-
+			for (int p1=0; p1<num_pics; p1++) {
 				if (gipf[p1] == NULL) {
 					break;
 				} else if (gipf[p1]->get_pixel(m, a_view, a_nick,
 						&c1[0], &c1[1], &c1[2]) == 0) {
 
-					double a1 = gipf[p1]->get_angle_off(a_view, a_nick);
-					double devign = ( 1 + a1 * a1 * V2 + a1*a1*a1*a1 * V1);
-					for (int l=0; l<3;l++) {
-						c1[l]=((uchar)MIN(rint((double) c1[l] * devign), 255));
+					for (int l = 0; l<3; l++) {	
+						c1d[l] = (double) c1[l];
 					}
 
-					for (int p2=0; p2<MAX_PICS; p2++) {
+					a1 = gipf[p1]->get_angle_off(a_view, a_nick);
 
+					for (int p2=1; p2<num_pics; p2++) {
 						if (p1 == p2) {
 							continue;
 						} else if (gipf[p2] == NULL) {
 							break;
 						} else if (gipf[p2]->get_pixel(m, a_view, a_nick,
 								&c2[0], &c2[1], &c2[2]) == 0) {
-							double a2 = gipf[p2]->get_angle_off(a_view, a_nick);
-							devign = ( 1 + a2 * a2 * V2 + a2*a2*a2*a2 * V1);
-							for (int l=0; l<3;l++) {
-								c2[l]=((uchar)MIN(rint((double) c2[l] * devign), 255));
-							}
 
-							if (fabs(a1-a2) < 0.01 &&
-								fabs(((double) c1[1]) / ((double) c1[0]) -
-									((double) c2[1]) / ((double) c2[0])) < 1.0 && 
-								fabs(((double) c1[2]) / ((double) c1[0]) -
-									((double) c2[2]) / ((double) c2[0])) < 1.0) {
-
-								if (sample_points[p1][p2] == NULL) {
-									fprintf(stderr, "allocated sample_points [%d][%d]\n", p1, p2);
-									sample_points[p1][p2] = (sample_t*)
-									calloc(500, sizeof(sample_t));
-								}
+							for (int l = 0; l<3; l++) {	
+								c2d[l] = (double) c2[l];
+							}
 
-								for (int k=0; k < 500; k++) {
-									if (sample_points[p1][p2][k].c1[0] == 0) {
-										for (int l=0; l < 3; l++) {
-											sample_points[p1][p2][k].c1[l] = c1[l];
-											sample_points[p1][p2][k].c2[l] = c2[l];
-										}
-										break;
+							a2 = gipf[p2]->get_angle_off(a_view, a_nick);
+
+							if (n_samples < max_samples &&
+								c1[0] < 200 && c1[1] < 200 && c1[2] < 200 &&
+								c2[0] < 200 && c2[1] < 200 && c2[2] < 200 &&
+								c1[0] > 20 && c1[1] > 20 && c1[2] > 20 &&
+								c2[0] > 20 && c2[1] > 20 && c2[2] > 20 &&
+								fabs(c1d[1] / c1d[0] - c2d[1] / c2d[0]) < 0.01 && 
+								fabs(c1d[2] / c1d[0] - c2d[2] / c2d[0]) < 0.01 ) {
+
+								for (int l = 0; l<3; l++) {	
+									gsl_matrix_set_zero(X);
+									if (p1 == 0) {
+fprintf(stderr, "==> p1 = 0 p2 = %d\n", p2);
+									//	gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 2, c2d[l] * a2 *a2 * a2 * a2);
+									//	gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 1, c2d[l] * a2 * a2);
+										gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 0, c2d[l]);
+										gsl_vector_set(yv, n_samples, c1d[l]);
+									} else {
+									//	gsl_matrix_set(X, n_samples, (p1-1) * 9 + 3 * l + 2, - c1d[l] * a1 *a1 * a1 * a1);
+									//	gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 2, c2d[l] * a2 * a2 * a2 * a2);
+									//	gsl_matrix_set(X, n_samples, (p1-1) * 9 + 3 * l + 1, - c1d[l] * a1 * a1);
+									//	gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 1, c2d[l] * a2 * a2);
+										gsl_matrix_set(X, n_samples, (p1-1) * 9 + 3 * l + 0, - c1d[l]);
+										gsl_matrix_set(X, n_samples, (p2-1) * 9 + 3 * l + 0, c2d[l]);
+										gsl_vector_set(yv, n_samples, 0.0);
 									}
-
+									n_samples++;
 								}
-			
+
 								if (merged_image) {
 									merged_image->set_pixel(x, 255, 0, 0);	
 									merged_pixel_set++;
 								}
-
-
 							}
+
 						}
 					}
-					
 
 					if (!merged_pixel_set && merged_image) {
 						merged_image->set_pixel(x, c1[0], c1[1], c1[2]);
+						merged_pixel_set++;
 					}
+
 				}
 			}
 		}
@@ -240,141 +206,47 @@ Stitch::color_calib(GipfelWidget::sample_mode_t m,
 	if (merged_image) {
 		merged_image->done();
 	}
-
-	for (int n=0; n<20; n++) {	
-		for (int i=0; i<MAX_PICS; i++) {
-			if (color_adjust[i][0] > 0.0) {
-				for (int j=1; j<MAX_PICS; j++) {
-					if (sample_points[i][j] != NULL && color_adjust[j][0] <= 0.0) {
-						get_color_adjust(sample_points[i][j], 500, color_adjust[i], color_adjust[j]);
-
-
-fprintf(stderr, "===> Adjust %d: %lf %lf %lf\n", j, color_adjust[j][0], color_adjust[j][1], color_adjust[j][2]);
-					}
+gsl_vector_fprintf(stderr, yv, "%f");
+	gsl_multifit_linear_workspace * work 
+           = gsl_multifit_linear_alloc (n_samples, n_vars);
+
+	ret = gsl_multifit_linear (X, yv, c, cov, &chisq, work);
+	gsl_multifit_linear_free (work);
+
+	fprintf(stderr, "gsl_multifit_linear returned %d\n", ret);
+
+	for (int p1=0; p1 < num_pics; p1++) {
+		for (int l = 0; l<3; l++) {	
+			if (p1 == 0) {
+				color_adjust[p1][l][0] = 1.0;
+				color_adjust[p1][l][1] = 0.0;
+				color_adjust[p1][l][2] = 0.0;
+			} else {
+				for (int k = 0; k<3; k++) {	
+					color_adjust[p1][l][k] = gsl_vector_get(c, (p1-1) * 9 + 3 * l + k);
 				}
-
+				color_adjust[p1][l][1] = 0.0;
+				color_adjust[p1][l][2] = 0.0;
 			}
 		}
-	}
-
+fprintf(stderr, "==> color_adjust(%d) %f %f %f\n", p1, color_adjust[p1][0][0], color_adjust[p1][1][0],color_adjust[p1][2][0]);
+	}		
 
 	return 0;
 }
 
+uchar
+Stitch::color_correct(uchar c, double a, int pic, int color) {
+	double cd = (double) c;
+	
+	cd = cd * (
+		color_adjust[pic][color][0] + 
+		color_adjust[pic][color][1] * a * a + 
+		color_adjust[pic][color][2] * a * a * a * a);
 
-int
-Stitch::vignette_calib(GipfelWidget::sample_mode_t m,
-	int w, int h, double view_start, double view_end) {
-
-	view_start = view_start * deg2rad;
-	view_end = view_end * deg2rad;
-
-	double step_view = (view_end - view_start) / w;
-	uchar r, g, b;
-	int y_off = h / 2;
-	int merged_pixel_set;
-	double radius = (double) w / (view_end -view_start);
-
-	int max_samples = 10000, n_samples = 0;
-	gsl_matrix *X, *cov;
-	gsl_vector *yv,  *c;
-	double chisq;
-
-	X = gsl_matrix_alloc(max_samples, 2);
-	yv = gsl_vector_alloc(max_samples);
-	c = gsl_vector_alloc(2);
-	cov = gsl_matrix_alloc (2, 2);
-
-	if (merged_image) {
-		merged_image->init(w, h);
-	}
-
-	for (int y=0; y<h; y++) {
-		double a_nick = atan((double)(y_off - y)/radius);
-
-		for (int x=0; x<w; x++) {
-			double a_view;
-			a_view = view_start + x * step_view;
-			merged_pixel_set = 0;
-			double l1 = 0.0, l2 = 0.0;
-			double a1, a2;
-			uchar c1[3], c2[3];
-
-			for (int i=0; i<MAX_PICS; i++) {
-
-				if (gipf[i] == NULL) {
-					break;
-				} else if (gipf[i]->get_pixel(m, a_view, a_nick,
-						&c2[0], &c2[1], &c2[2]) == 0) {
-					l2 = (double) c2[0] * color_adjust[i][0] +
-						(double) c2[1] * color_adjust[i][1] +
-						(double) c2[2] * color_adjust[i][2];
-					a2 = gipf[i]->get_angle_off(a_view, a_nick);
-
-					if (l1 > 0.0) {
-
-						if (n_samples < max_samples &&
-							fabs(a1 - a2) > 0.02 &&
-							c1[0] < 200 && c1[1] < 200 && c1[2] < 200 &&
-							c2[0] < 200 && c2[1] < 200 && c2[2] < 200 &&
-							fabs(((double) c1[1]) / ((double) c1[0]) -
-								((double) c2[1]) / ((double) c2[0])) < 0.05 && 
-							fabs(((double) c1[2]) / ((double) c1[0]) -
-								((double) c2[2]) / ((double) c2[0])) < 0.05) {
-
-							gsl_matrix_set(X, n_samples, 0, l2 *  a2 * a2 *a2 * a2 - l1 * a1 *a1 *a1 * a1);
-							gsl_matrix_set(X, n_samples, 1, l2 * a2 * a2 - l1 * a1 * a1);
-							gsl_vector_set(yv, n_samples, l1 - l2);
-							n_samples++;
-
-							if (merged_image) {
-								double d = fabs(V2 * (l2 * a2 * a2 - l1 * a1 * a1) + V1 * (l2 *  a2 * a2 *a2 - l1 * a1 * a1 *a1) - (l1 -l2)); 
-
-								merged_image->set_pixel(x, (uchar) rint(d * 2), 0, 0);	
-								merged_pixel_set++;
-							}
-
-
-						}
-					}
-						
-					l1 = l2;
-					a1 = a2;
-					c1[0] = c2[0];
-					c1[1] = c2[1];
-					c1[2] = c2[2];
-
-					if (!merged_pixel_set && merged_image) {
-						merged_image->set_pixel(x, c1[0], c1[1], c1[2]);
-					}
-				}
-			}
-		}
-		if (merged_image) {
-			merged_image->next_line();
-		}
-	}
-
-	if (merged_image) {
-		merged_image->done();
-	}
-
-	gsl_multifit_linear_workspace * work 
-           = gsl_multifit_linear_alloc (n_samples, 2);
-
-	gsl_multifit_linear (X, yv, c, cov, &chisq, work);
-         gsl_multifit_linear_free (work);
-
-	V1 = gsl_vector_get(c,0);
-	V2 = gsl_vector_get(c,1);
-
-	fprintf(stderr, "===>  v1 %lf v2 %lf, (i %d)\n", 
-		V1, V2, n_samples);
-
-	return 0;
+	return (uchar) MIN(rint(cd), 255.0);
 }
 
-
 int
 Stitch::resample(GipfelWidget::sample_mode_t m,
 	int w, int h, double view_start, double view_end) {
@@ -410,20 +282,11 @@ Stitch::resample(GipfelWidget::sample_mode_t m,
 					break;
 				} else if (gipf[i]->get_pixel(m, a_view, a_nick,
 						&r, &g, &b) == 0) {
-					double a2 = gipf[i]->get_angle_off(a_view, a_nick);
-					double devign = ( 1 + a2 * a2 * V2 + a2*a2*a2*a2 * V1);
-
-#if 1 
-					r = (uchar) MIN(rint(((double) r) * devign * color_adjust[i][0]), 255);
-					g = (uchar) MIN(rint(((double) g) * devign * color_adjust[i][1]), 255);
-					b = (uchar) MIN(rint(((double) b) * devign * color_adjust[i][2]), 255);
-#else 
-
-					r = (uchar) MIN(rint((double) 100 * devign), 255);
-					g = (uchar) MIN(rint((double) 100 * devign), 255);
-					b = (uchar) MIN(rint((double) 100 * devign), 255);
-#endif
+					double a = gipf[i]->get_angle_off(a_view, a_nick);
 
+					r = color_correct(r, a, i, 0);
+					g = color_correct(g, a, i, 1);
+					b = color_correct(b, a, i, 2);
 
 					if (single_images[i]) {
 						single_images[i]->set_pixel(x, r, g, b);
diff --git a/src/gipfel.cxx b/src/gipfel.cxx
index 23a5c5f..e9743d8 100644
--- a/src/gipfel.cxx
+++ b/src/gipfel.cxx
@@ -568,10 +568,6 @@ stitch(GipfelWidget::sample_mode_t m,
 		st->load_image(argv[i]);
 	}
 
-	if (type & STITCH_VIGNETTE_CALIB) {
-		st->color_calib(m, stitch_w, stitch_h, from, to);
-		//st->vignette_calib(m, stitch_w, stitch_h, from, to);
-	}
 
 	if (type & STITCH_JPEG) {
 
@@ -604,6 +600,9 @@ stitch(GipfelWidget::sample_mode_t m,
 		win->show(0, argv); 
 		st->set_output((OutputImage*) img);
 
+		if (type & STITCH_VIGNETTE_CALIB) {
+			st->vignette_calib(m, stitch_w, stitch_h, from, to);
+		}
 		st->resample(m, stitch_w, stitch_h, from, to);
 
 		img->redraw();