merge

25 files changed, 287 insertions, 195 deletions

diff --git a/src/GipfelWidget.H b/src/GipfelWidget.H
index 7d1b5ba..45bfaee 100644
--- a/src/GipfelWidget.H
+++ b/src/GipfelWidget.H
@@ -39,18 +39,18 @@ class GipfelWidget : public Fl_Widget {
 		int get_rel_track_width(Hill *m);
 
 		static int get_pixel_nearest(Fl_Image *img, double x, double y,
-			uchar *r, uchar *g, uchar *b);
+			int *r, int *g, int *b);
 
-		static int get_pixel_bilinear(Fl_Image *img, double x, double y,
-			uchar *r, uchar *g, uchar *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,
-			uchar *r, uchar *g, uchar *b);
+			int *r, int *g, int *b);
 
 	public:
 		typedef enum {
 			NEAREST  = 0,
-			BILINEAR = 1
+			BICUBIC  = 1
 		} sample_mode_t;
 	
 		GipfelWidget(int X,int Y,int W, int H);
@@ -115,16 +115,18 @@ class GipfelWidget : public Fl_Widget {
 
 		void set_projection(ProjectionLSQ::Projection_t p);
 
-		void get_distortion_params(double *k0, double *k1);
+		void get_distortion_params(double *k0, double *k1i, double *x0);
 
-		void set_distortion_params(double k0, double k1);
+		void set_distortion_params(double k0, double k1, double x0);
 
 		Hills *get_mountains();
 
 		int comp_params();
 
 		int get_pixel(GipfelWidget::sample_mode_t,
-			double a_alph, double a_nick, uchar *r, uchar *g, uchar *b);
+			double a_alph, double a_nick, int *r, int *g, int *b);
+
+		double get_angle_off(double view, double nick);
 
 		int get_distortion_profile_name(char *buf, int buflen);
 
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx
index 2127e64..2563d73 100644
--- a/src/GipfelWidget.cxx
+++ b/src/GipfelWidget.cxx
@@ -31,6 +31,7 @@
 
 
 #define MAX(A,B) ((A)>(B)?(A):(B))
+#define MIN(A,B) ((A)<(B)?(A):(B))
 
 static double pi_d, deg2rad;
 
@@ -99,15 +100,20 @@ GipfelWidget::load_image(char *file) {
 	// 1. gipfel data in JPEG comment
 	// 2. matching distortion profile
 	// 3. set the to 0.0, 0.0
-	md->get_distortion_params(&pan->parms.k0, &pan->parms.k1);
+	md->get_distortion_params(&pan->parms.k0, &pan->parms.k1, &pan->parms.x0);
 	if (isnan(pan->parms.k0)) {
 		char buf[1024];
 		get_distortion_profile_name(buf, sizeof(buf));
 		load_distortion_params(buf);
 		if (isnan(pan->parms.k0)) {
 			pan->parms.k0 = 0.0;
+		}
+		if (isnan(pan->parms.k1)) {
 			pan->parms.k1 = 0.0;
 		}
+		if (isnan(pan->parms.x0)) {
+			pan->parms.x0 = 0.0;
+		}
 	}
 
 	return 0;
@@ -138,7 +144,7 @@ GipfelWidget::save_image(char *file) {
 	md->set_tilt(get_tilt_angle());
 	md->set_focal_length_35mm(get_focal_length_35mm());
 	md->set_projection_type((int) get_projection());
-	md->set_distortion_params(pan->parms.k0, pan->parms.k1);
+	md->set_distortion_params(pan->parms.k0, pan->parms.k1, pan->parms.x0);
 
 	ret = md->save_image(img_file, file);
 	delete md;
@@ -473,13 +479,13 @@ GipfelWidget::set_projection(ProjectionLSQ::Projection_t p) {
 }
 
 void
-GipfelWidget::get_distortion_params(double *k0, double *k1) {
-	pan->get_distortion_params(k0, k1);
+GipfelWidget::get_distortion_params(double *k0, double *k1, double *x0) {
+	pan->get_distortion_params(k0, k1, x0);
 }
 
 void
-GipfelWidget::set_distortion_params(double k0, double k1) {
-	pan->set_distortion_params(k0, k1);
+GipfelWidget::set_distortion_params(double k0, double k1, double x0) {
+	pan->set_distortion_params(k0, k1, x0);
 	redraw();
 }
 
@@ -669,8 +675,10 @@ GipfelWidget::handle(int event) {
 
 int
 GipfelWidget::get_pixel(GipfelWidget::sample_mode_t m,
-	double a_alph, double a_nick, uchar *r, uchar *g, uchar *b) {
+	double a_alph, double a_nick, int *r, int *g, int *b) {
 	double px, py;
+	int r_tmp, g_tmp, b_tmp;
+	int ret;
 
 	if (img == NULL) {
 		return 1;
@@ -680,18 +688,28 @@ GipfelWidget::get_pixel(GipfelWidget::sample_mode_t m,
 		return 1;
 	}
 
-	if (m == GipfelWidget::BILINEAR) {
-		return get_pixel_bilinear(img, px + ((double) img->w()) / 2.0,
+	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 {
-		return get_pixel_nearest(img, px + ((double) img->w()) / 2.0,
+		ret = get_pixel_nearest(img, px + ((double) img->w()) / 2.0,
 			py + ((double) img->h()) / 2.0, r, g, b);
 	}
+
+	return ret;
+}
+
+double
+GipfelWidget::get_angle_off(double a_alph, double a_nick) {
+	double angle =  atan(pow(pow(tan(a_alph - pan->parms.a_center), 2.0) +
+			pow(tan(a_nick - pan->parms.a_nick), 2.0), 0.5));
+
+	return angle;
 }
 
 int
 GipfelWidget::get_pixel_nearest(Fl_Image *img, double x, double y,
-	uchar *r, uchar *g, uchar *b) {
+	int *r, int *g, int *b) {
 	if (isnan(x) || isnan(y)) {
 		return 1;
 	} else {
@@ -699,48 +717,56 @@ GipfelWidget::get_pixel_nearest(Fl_Image *img, double x, double y,
 	}
 }
 
-int
-GipfelWidget::get_pixel_bilinear(Fl_Image *img, double x, double y,
-	uchar *r, uchar *g, uchar *b) {
-	uchar a_r[4] = {0, 0, 0, 0};
-	uchar a_g[4] = {0, 0, 0, 0};
-	uchar a_b[4] = {0, 0, 0, 0};
-	float v0 , v1;
-	int fl_x = (int) floor(x);
-	int fl_y = (int) floor(y);
-	int i, n;
-
-	i = 0;
-	n = 0;
-	for (int iy = 0; iy <= 1; iy++) {
-		for (int ix = 0; ix <= 1; ix++) {
-			n += get_pixel(img, fl_x + ix, fl_y + iy, &(a_r[i]), &(a_g[i]), &(a_b[i]));
-			i++;
-		}
-	}
+static inline double
+interp_cubic(double x, double x2, double x3, double *v) {
+	double a0, a1, a2, a3;
 
-	v0 = a_r[0] * (1 - (x - fl_x)) + a_r[1] * (x - fl_x);
-	v1 = a_r[2] * (1 - (x - fl_x)) + a_r[3] * (x - fl_x);
-	*r = (uchar) rint(v0 * (1 - (y - fl_y)) + v1 * (y - fl_y));
+	a0 = v[3] - v[2] - v[0] + v[1];
+	a1 = v[0] - v[1] - a0;
+	a2 = v[2] - v[0];
+	a3 = v[1];
 
-	v0 = a_g[0] * (1 - (x - fl_x)) + a_g[1] * (x - fl_x);
-	v1 = a_g[2] * (1 - (x - fl_x)) + a_g[3] * (x - fl_x);
-	*g = (uchar) rint(v0 * (1 - (y - fl_y)) + v1 * (y - fl_y));
+	return a0 * x3 + a1 * x2 + a2 * x + a3;
+}
 
-	v0 = a_b[0] * (1 - (x - fl_x)) + a_b[1] * (x - fl_x);
-	v1 = a_b[2] * (1 - (x - fl_x)) + a_b[3] * (x - fl_x);
-	*b = (uchar) rint(v0 * (1 - (y - fl_y)) + v1 * (y - fl_y));
+int
+GipfelWidget::get_pixel_bicubic(Fl_Image *img, double x, double y,
+	int *r, int *g, int *b) {
+
+	double fl_x = floor(x);
+	double fl_y = floor(y);
+	double dx = x - fl_x, dx2 = dx * dx, dx3 = dx2 * dx;
+	double dy = y - fl_y, dy2 = dy * dy, dy3 = dy2 * dy;
+	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, (int) fl_x + ix - 1, (int) fl_y + iy - 1,
+					&ic[0], &ic[1], &ic[2]) != 0) {
+				return 1;
+			}
 
-	if (n >= 1) {
-		return 1;
-	} else {
-		return 0;
+			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(dx, dx2, dx3, c[l]);
+		}
 	}
+
+	*r = (int) rint(interp_cubic(dy, dy2, dy3, c1[0]));
+	*g = (int) rint(interp_cubic(dy, dy2, dy3, c1[1]));
+	*b = (int) rint(interp_cubic(dy, dy2, dy3, c1[2]));
+	return 0;
 }
 
 int
 GipfelWidget::get_pixel(Fl_Image *img, int x, int y,
-                     uchar *r, uchar *g, uchar *b) {
+                     int *r, int *g, int *b) {
 	if ( img->d() == 0 ) {
 		return 1;
 	}
@@ -752,16 +778,16 @@ GipfelWidget::get_pixel(Fl_Image *img, int x, int y,
 	switch (img->count()) {
 		case 1:
 		{                                            // bitmap
-			const char *buf = img->data()[0];
+			const unsigned char *buf = (const unsigned char*) img->data()[0];
 			switch (img->d())
 			{
 				case 1:
 					*r = *g = *b = *(buf+index);
 					break;
 				case 3:                              // 24bit
-					*r = *(buf+index+0);
-					*g = *(buf+index+1);
-					*b = *(buf+index+2);
+					*r = (int) *(buf+index+0);
+					*g = (int) *(buf+index+1);
+					*b = (int) *(buf+index+2);
 					break;
 				default:                             // ??
 					printf("Not supported: chans=%d\n", img->d());
@@ -774,8 +800,11 @@ GipfelWidget::get_pixel(Fl_Image *img, int x, int y,
 		return 1;
 	}
 
-	return 0;
+	*r = *r * 255;
+	*g = *g * 255;
+	*b = *b * 255;
 
+	return 0;
 }
 
 int
@@ -802,6 +831,7 @@ GipfelWidget::load_distortion_params(const char *prof_name) {
 	Fl_Preferences prof(dist_prefs, prof_name);
 	ret += prof.get("k0", pan->parms.k0, pan->parms.k0);
 	ret += prof.get("k1", pan->parms.k1, pan->parms.k1);
+	ret += prof.get("x0", pan->parms.x0, pan->parms.x0);
 
 	return !ret;
 }
@@ -818,6 +848,7 @@ GipfelWidget::save_distortion_params(const char *prof_name, int force) {
 
 	prof.set("k0", pan->parms.k0);
 	prof.set("k1", pan->parms.k1);
+	prof.set("x0", pan->parms.x0);
 
 	return 0;
 }
diff --git a/src/ImageMetaData.H b/src/ImageMetaData.H
index 27cbbd8..b7a5959 100644
--- a/src/ImageMetaData.H
+++ b/src/ImageMetaData.H
@@ -19,6 +19,7 @@ class ImageMetaData {
 		double tilt;
 		double k0;
 		double k1;
+		double x0;
 		double focal_length;
 		double focal_length_35mm;
 		double scale;
@@ -47,7 +48,7 @@ class ImageMetaData {
 		double get_focal_length();
 		double get_focal_length_35mm();
 		int    get_projection_type();
-		void get_distortion_params(double *_k0, double *_k1);
+		void get_distortion_params(double *_k0, double *_k1, double *_x0);
 
 		void set_longitude(double v);
 		void set_latitude(double v);
@@ -57,6 +58,6 @@ class ImageMetaData {
 		void set_tilt(double v);
 		void set_focal_length_35mm(double v);
 		void set_projection_type(int v);
-		void set_distortion_params(double _k0, double _k1);
+		void set_distortion_params(double _k0, double _k1, double _x0);
 };
 #endif
diff --git a/src/ImageMetaData.cxx b/src/ImageMetaData.cxx
index e9a1241..5f619d2 100644
--- a/src/ImageMetaData.cxx
+++ b/src/ImageMetaData.cxx
@@ -38,6 +38,7 @@ ImageMetaData::clear() {
 	tilt = NAN;
 	k0 = NAN;
 	k1 = NAN;
+	x0 = NAN;
 	focal_length = NAN;
 	focal_length_35mm = NAN;
 	scale = NAN;
@@ -159,7 +160,7 @@ ImageMetaData::load_image_exif(char *name) {
 
 
 #define GIPFEL_FORMAT_1 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, scale %lf, projection type %d"
-#define GIPFEL_FORMAT_2 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, focal_length_35mm %lf, projection type %d, k0 %lf, k1 %lf"
+#define GIPFEL_FORMAT_2 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, focal_length_35mm %lf, projection type %d, k0 %lf, k1 %lf, x0 %lf"
 
 int
 ImageMetaData::load_image_jpgcom(char *name) {
@@ -168,7 +169,7 @@ ImageMetaData::load_image_jpgcom(char *name) {
     pid_t pid;
     int status;
     char buf[1024];
-    double lo, la, he, dir, ni, ti, fr, _k0, _k1;
+    double lo, la, he, dir, ni, ti, fr, _k0, _k1, _x0 = 0.0;
     int pt = 0;
     int n, ret = 1;
 
@@ -181,7 +182,7 @@ ImageMetaData::load_image_jpgcom(char *name) {
     if (p) {
         while (fgets(buf, sizeof(buf), p) != NULL) {
             if ((n = sscanf(buf, GIPFEL_FORMAT_2,
-                    &lo, &la, &he, &dir, &ni, &ti, &fr, &pt, &_k0, &_k1)) >= 8) {
+                    &lo, &la, &he, &dir, &ni, &ti, &fr, &pt, &_k0, &_k1, &_x0)) >= 8) {
 
                 longitude = lo;
                 latitude  = la;
@@ -195,6 +196,7 @@ ImageMetaData::load_image_jpgcom(char *name) {
 				if (n >= 10) {
 					k0 = _k0;
 					k1 = _k1;
+					x0 = _x0;
 				}
 
                 ret = 0;
@@ -267,7 +269,7 @@ ImageMetaData::save_image_jpgcom(char *in_img, char *out_img) {
         tilt,
         focal_length_35mm,
         projection_type,
-		k0, k1);
+		k0, k1, x0);
 
     // try to save gipfel data in JPEG comment section
     args[0] = "wrjpgcom";
@@ -395,13 +397,15 @@ ImageMetaData::set_projection_type(int v) {
 }
 
 void
-ImageMetaData::get_distortion_params(double *_k0, double *_k1) {
+ImageMetaData::get_distortion_params(double *_k0, double *_k1, double *_x0) {
 	*_k0 = k0;	
 	*_k1 = k1;	
+	*_x0 = x0;	
 }
 
 void
-ImageMetaData::set_distortion_params(double _k0, double _k1) {
+ImageMetaData::set_distortion_params(double _k0, double _k1, double _x0) {
 	k0 = _k0;	
 	k1 = _k1;	
+	x0 = _x0;	
 }
diff --git a/src/JPEGOutputImage.H b/src/JPEGOutputImage.H
index b3554c8..06499ca 100644
--- a/src/JPEGOutputImage.H
+++ b/src/JPEGOutputImage.H
@@ -15,7 +15,7 @@ extern "C" {
 
 #include "OutputImage.H"
 
-class JPEGOutputImage : OutputImage {
+class JPEGOutputImage : public OutputImage {
 	private:
 		unsigned char *row;
 		char *file;
@@ -32,7 +32,7 @@ class JPEGOutputImage : OutputImage {
 	protected:
 		int init_internal(int w, int h);
 
-		int set_pixel_internal(int x, char r, char g, char b);	
+		int set_pixel_internal(int x, int r, int g, int b);	
 
 		int next_line_internal();
 
diff --git a/src/JPEGOutputImage.cxx b/src/JPEGOutputImage.cxx
index 6e685d0..0a5eb00 100644
--- a/src/JPEGOutputImage.cxx
+++ b/src/JPEGOutputImage.cxx
@@ -68,10 +68,10 @@ JPEGOutputImage::init_internal(int w1, int h1) {
 }
 
 int
-JPEGOutputImage::set_pixel_internal(int x, char r, char g, char b) {
-	row[x*3+0] = r;
-	row[x*3+1] = g;
-	row[x*3+2] = b;
+JPEGOutputImage::set_pixel_internal(int x, int r, int g, int b) {
+	row[x*3+0] = (unsigned char) (r / 255);
+	row[x*3+1] = (unsigned char) (g / 255);
+	row[x*3+2] = (unsigned char) (b / 255);
 
 	return 0;
 }
diff --git a/src/OutputImage.H b/src/OutputImage.H
index 7c9f0fd..17955b9 100644
--- a/src/OutputImage.H
+++ b/src/OutputImage.H
@@ -19,7 +19,7 @@ class OutputImage {
 		
 		virtual int init(int w1, int h1);
 
-		int set_pixel(int x, char r, char g, char b);
+		int set_pixel(int x, int r, int g, int b);
 
 		int next_line();
 
@@ -30,7 +30,7 @@ class OutputImage {
 
 		virtual int init_internal(int w1, int h1);
 
-		virtual int set_pixel_internal(int x, char r, char g, char b);	
+		virtual int set_pixel_internal(int x, int r, int g, int b);	
 	
 		virtual int next_line_internal();
 
diff --git a/src/OutputImage.cxx b/src/OutputImage.cxx
index ddf8e14..f981f18 100644
--- a/src/OutputImage.cxx
+++ b/src/OutputImage.cxx
@@ -35,7 +35,7 @@ OutputImage::init_internal(int w1, int h1) {
 }
 
 int
-OutputImage::set_pixel(int x, char r, char g, char b) {
+OutputImage::set_pixel(int x, int r, int g, int b) {
 	if (!initialized || x < 0 || x >= W) {
 		return 1;
 	} else {
@@ -44,7 +44,7 @@ OutputImage::set_pixel(int x, char r, char g, char b) {
 }
 
 int
-OutputImage::set_pixel_internal(int x, char r, char g, char b) {
+OutputImage::set_pixel_internal(int x, int r, int g, int b) {
 	return 0;
 }
 
diff --git a/src/Panorama.H b/src/Panorama.H
index 651fe8e..d5ebb37 100644
--- a/src/Panorama.H
+++ b/src/Panorama.H
@@ -41,7 +41,7 @@ class Panorama {
 
 		double alpha(double phi, double lam);
 
-		double nick(double dist, double height);
+		double nick(Hill *m);
 
 		double comp_center_angle(double alph_a, double alph_b, double d1, double d2);
 
@@ -51,8 +51,6 @@ class Panorama {
 
 		int optimize(Hill *m1, Hill *m2);
 
-		double get_value(Hills *p);
-
 		int is_visible(double a_alph);
 
 		double pi_d, deg2rad;
@@ -126,9 +124,9 @@ class Panorama {
 
 		void set_projection(ProjectionLSQ::Projection_t p);
 
-		void get_distortion_params(double *k0, double *k1);
+		void get_distortion_params(double *k0, double *k1, double *x0);
 
-		void set_distortion_params(double k0, double k1);
+		void set_distortion_params(double k0, double k1, double x0);
 
 		int get_coordinates(double a_alph, double a_nick, double *x, double *y);
 };
diff --git a/src/Panorama.cxx b/src/Panorama.cxx
index 1e3015a..27e4579 100644
--- a/src/Panorama.cxx
+++ b/src/Panorama.cxx
@@ -29,12 +29,16 @@ Panorama::Panorama() {
 	parms.scale = 3500.0;
 	parms.k0 = 0.0;
 	parms.k1 = 0.0;
+	parms.x0 = 0.0;
 	view_name = NULL;
 	view_phi = 0.0;
 	view_lam = 0.0;
 	view_height = 0.0;
 	proj = NULL;
 	set_projection(ProjectionLSQ::RECTILINEAR);
+
+
+fprintf(stderr, "=> %f, %f\n", get_earth_radius(0.0), get_earth_radius(pi_d/2.0));
 }
 
 Panorama::~Panorama() {
@@ -129,34 +133,6 @@ Panorama::get_visible_mountains() {
 	return visible_mountains;
 }
 
-double
-Panorama::get_value(Hills *p) {
-	int i, j;
-	double v = 0.0, d_min, d;
-
-	if (isnan(parms.scale) || isnan(parms.a_center) || isnan(parms.a_tilt) || isnan(parms.a_nick) ||
-		parms.scale < 500.0 || parms.scale > 100000.0 || 
-		parms.a_nick > pi_d/4.0 || parms.a_nick < - pi_d/4.0 || 
-		parms.a_tilt > pi_d/16.0 || parms.a_tilt < - pi_d/16.0) {
-		return 10000000.0;
-	}
-
-
-	for (i=0; i<p->get_num(); i++) {
-		d_min = 1000.0;
-		for (j=0; j<visible_mountains->get_num(); j++) {
-			d = pow(p->get(i)->x - visible_mountains->get(j)->x, 2.0) + 
-				pow(p->get(i)->y - visible_mountains->get(j)->y, 2.0);
-			if (d < d_min) {
-				d_min = d;
-			}
-		}
-		v = v + d_min;
-	}
-
-	return v;
-}
-
 int
 Panorama::comp_params(Hills *h) {
 	int ret;
@@ -191,15 +167,17 @@ Panorama::set_scale(double s) {
 }
 
 void
-Panorama::get_distortion_params(double *k0, double *k1) {
+Panorama::get_distortion_params(double *k0, double *k1, double *x0) {
 	*k0 = parms.k0;
 	*k1 = parms.k1;
+	*x0 = parms.x0;
 }
 
 void
-Panorama::set_distortion_params(double k0, double k1) {
+Panorama::set_distortion_params(double k0, double k1, double x0) {
 	parms.k0 = k0;
 	parms.k1 = k1;
+	parms.x0 = x0;
 	update_coordinates();
 }
 
@@ -328,7 +306,7 @@ Panorama::update_angles() {
 		if (m->phi != view_phi || m->lam != view_lam) {
 
 			m->alph = alpha(m->phi, m->lam);
-			m->a_nick = nick(m->dist, m->height);
+			m->a_nick = nick(m);
 		}
 	}
 
@@ -394,7 +372,7 @@ Panorama::update_close_mountains() {
 
 		if (m->flags & Hill::TRACK_POINT ||
 			((m->phi != view_phi || m->lam != view_lam) &&
-			 (m->height / (m->dist * get_earth_radius(m->phi)) 
+			 (m->height / (m->dist * EARTH_RADIUS) 
 			  > height_dist_ratio))) {
 
 			close_mountains->add(m);
@@ -456,29 +434,37 @@ Panorama::alpha(double phi, double lam) {
 
 
 double
-Panorama::nick(double dist, double height) {
+Panorama::nick(Hill *m) {
 	double a, b, c;
 	double beta;
 
-	b = height + get_earth_radius(view_phi);
+	b = m->height + get_earth_radius(m->phi);
 	c = view_height + get_earth_radius(view_phi);
 
-	a = pow(((b * (b - (2.0 * c * cos(dist)))) + (c * c)), (1.0 / 2.0));
+	a = pow(((b * (b - (2.0 * c * cos(m->dist)))) + (c * c)), (1.0 / 2.0));
 	beta = acos((-(b*b) + (a*a) + (c*c))/(2 * a * c));
 
 	return beta - pi_d / 2.0;
 }
 
+// return local distance to center of WGS84 ellipsoid
 double
-Panorama::get_earth_radius(double latitude) {
-	return EARTH_RADIUS;
+Panorama::get_earth_radius(double phi) {
+	double a = 6378137.000;
+	double b = 6356752.315;
+	double r;
+	double ata = tan(phi);
+return EARTH_RADIUS;
+	r = a*pow(pow(ata,2)+1,1.0/2.0)*fabs(b)*pow(pow(b,2)+pow(a,2)*pow(ata,2),-1.0/2.0);
+
+	return r;
 }
 
 double
 Panorama::get_real_distance(Hill *m) {
 	double a, b, c, gam;
 
-	a = view_height + get_earth_radius(m->phi);
+	a = view_height + get_earth_radius(view_phi);
 	b = m->height + get_earth_radius(m->phi); 
 	gam = m->dist;
 
@@ -500,8 +486,6 @@ Panorama::is_visible(double a_alph) {
 
 int
 Panorama::get_coordinates(double a_alph, double a_nick, double *x, double *y) {
-
-
 	if (is_visible(a_alph)) {
 		proj->get_coordinates(a_alph, a_nick, &parms, x, y);
 		return 0;
diff --git a/src/PreviewOutputImage.H b/src/PreviewOutputImage.H
index 8b41684..53b58db 100644
--- a/src/PreviewOutputImage.H
+++ b/src/PreviewOutputImage.H
@@ -14,7 +14,7 @@
 
 #include "OutputImage.H"
 
-class PreviewOutputImage : OutputImage , public Fl_Widget {
+class PreviewOutputImage : public OutputImage , public Fl_Widget {
 	private:
 		uchar *data;
 		int d;
@@ -29,7 +29,7 @@ class PreviewOutputImage : OutputImage , public Fl_Widget {
 	protected:
 		int init_internal(int w, int h);
 
-		int set_pixel_internal(int x, char r, char g, char b);	
+		int set_pixel_internal(int x, int r, int g, int b);	
 
 		int next_line_internal();
 
diff --git a/src/PreviewOutputImage.cxx b/src/PreviewOutputImage.cxx
index 86dadee..0ae182e 100644
--- a/src/PreviewOutputImage.cxx
+++ b/src/PreviewOutputImage.cxx
@@ -35,15 +35,15 @@ PreviewOutputImage::init_internal(int w, int h) {
 
 
 int
-PreviewOutputImage::set_pixel_internal(int x, char r, char g, char b) {
+PreviewOutputImage::set_pixel_internal(int x, int r, int g, int b) {
 	if (!data) {
 		return 1;
 	}
 
 	long index = (line * w() * d + (x * d));
-	*(data+index+0) = r;
-	*(data+index+1) = g;
-	*(data+index+2) = b;
+	*(data+index+0) = (unsigned char) (r / 255);
+	*(data+index+1) = (unsigned char) (g / 255);
+	*(data+index+2) = (unsigned char) (b / 255);
 
 	return 0;
 }
diff --git a/src/ProjectionCylindrical.H b/src/ProjectionCylindrical.H
index 432d868..2cb9687 100644
--- a/src/ProjectionCylindrical.H
+++ b/src/ProjectionCylindrical.H
@@ -14,7 +14,7 @@ class ProjectionCylindrical : public ProjectionLSQ {
 
 		virtual double get_view_angle() {return 6.2831853;}; /* 360 deg */
 
-#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick
 		virtual double mac_x(ARGS);
 		virtual double mac_y(ARGS);
         virtual double mac_x_dc_view(ARGS);
diff --git a/src/ProjectionLSQ.H b/src/ProjectionLSQ.H
index 54ab641..273585f 100644
--- a/src/ProjectionLSQ.H
+++ b/src/ProjectionLSQ.H
@@ -34,7 +34,7 @@ class ProjectionLSQ {
 
 		virtual double get_view_angle();
 
-#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick
 		virtual double mac_x(ARGS);
 		virtual double mac_y(ARGS);
 		virtual double mac_x_dc_view(ARGS);
@@ -43,12 +43,14 @@ class ProjectionLSQ {
 		virtual double mac_x_dscale(ARGS);
 		virtual double mac_x_dk0(ARGS);
 		virtual double mac_x_dk1(ARGS);
+		virtual double mac_x_dx0(ARGS);
 		virtual double mac_y_dc_view(ARGS);
 		virtual double mac_y_dc_nick(ARGS);
 		virtual double mac_y_dc_tilt(ARGS);
 		virtual double mac_y_dscale(ARGS);
 		virtual double mac_y_dk0(ARGS);
 		virtual double mac_y_dk1(ARGS);
+		virtual double mac_y_dx0(ARGS);
 #undef ARGS
 
 };
diff --git a/src/ProjectionLSQ.cxx b/src/ProjectionLSQ.cxx
index 077641b..da5b3da 100644
--- a/src/ProjectionLSQ.cxx
+++ b/src/ProjectionLSQ.cxx
@@ -45,6 +45,7 @@ ProjectionLSQ::comp_params(const Hills *h, ViewParams *parms) {
 		distortion_correct = 1;
 		parms->k0 = 0.0;
 		parms->k1 = 0.0;
+		parms->x0 = 0.0;
 	}
 
 	m1 = h->get(0);
@@ -68,7 +69,7 @@ ProjectionLSQ::comp_params(const Hills *h, ViewParams *parms) {
 	if (distortion_correct) {
 		lsq(h, parms, 1);
 	}
-
+fprintf(stderr, "===> x0 %f\n", parms->x0);
 	return 0;
 }
 
@@ -79,12 +80,12 @@ struct data {
 	const ViewParams *old_params;
 };
 
-#define CALL(A) dat->p->A(c_view, c_nick, c_tilt, scale, k0, k1, m->alph, m->a_nick) 
+#define CALL(A) dat->p->A(c_view, c_nick, c_tilt, scale, k0, k1, x0, m->alph, m->a_nick) 
 
 static int
 lsq_f (const gsl_vector * x, void *data, gsl_vector * f) {
 	struct data *dat = (struct data *) data;
-	double c_view, c_nick, c_tilt, scale, k0, k1;
+	double c_view, c_nick, c_tilt, scale, k0, k1, x0;
 
 	c_view = gsl_vector_get (x, 0);
 	c_nick = gsl_vector_get (x, 1);
@@ -93,9 +94,11 @@ lsq_f (const gsl_vector * x, void *data, gsl_vector * f) {
 	if (dat->distortion_correct) {
 		k0  = gsl_vector_get (x, 4);
 		k1 = gsl_vector_get (x, 5);
+		x0 = gsl_vector_get (x, 6);
 	} else {
 		k0 = dat->old_params->k0;
 		k1 = dat->old_params->k1;
+		x0 = dat->old_params->x0;
 	}
 
 	for (int i=0; i<dat->h->get_num(); i++) {
@@ -115,7 +118,7 @@ lsq_f (const gsl_vector * x, void *data, gsl_vector * f) {
 static int
 lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) {
 	struct data *dat = (struct data *) data;
-    double c_view, c_nick, c_tilt, scale, k0, k1;
+    double c_view, c_nick, c_tilt, scale, k0, k1, x0;
 
     c_view = gsl_vector_get (x, 0);
     c_nick = gsl_vector_get (x, 1);
@@ -124,9 +127,11 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) {
     if (dat->distortion_correct) {
         k0  = gsl_vector_get (x, 4);
         k1 = gsl_vector_get (x, 5);
+        x0 = gsl_vector_get (x, 6);
     } else {
         k0 = dat->old_params->k0;
         k1 = dat->old_params->k1;
+        x0 = dat->old_params->x0;
     }            
 
 	for (int i=0; i<dat->h->get_num(); i++) {
@@ -139,6 +144,7 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) {
 		if (dat->distortion_correct) {
 			gsl_matrix_set (J, 2*i, 4, CALL(mac_x_dk0));
 			gsl_matrix_set (J, 2*i, 5, CALL(mac_x_dk1));
+			gsl_matrix_set (J, 2*i, 6, CALL(mac_x_dx0));
 		}
 
 		gsl_matrix_set (J, 2*i+1, 0, CALL(mac_y_dc_view));
@@ -148,6 +154,7 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) {
 		if (dat->distortion_correct) {
 			gsl_matrix_set (J, 2*i+1, 4, CALL(mac_y_dk0));
 			gsl_matrix_set (J, 2*i+1, 5, CALL(mac_y_dk1));
+			gsl_matrix_set (J, 2*i+1, 6, CALL(mac_y_dx0));
 		}
 	}
 
@@ -173,7 +180,7 @@ ProjectionLSQ::lsq(const Hills *h, ViewParams *parms,
 	double x_init[8];
 	gsl_vector_view x;
 	int status;
-	int num_params = distortion_correct?6:4;
+	int num_params = distortion_correct?7:4;
 
 	dat.p = this;
 	dat.distortion_correct = distortion_correct;
@@ -186,6 +193,7 @@ ProjectionLSQ::lsq(const Hills *h, ViewParams *parms,
 	x_init[3] = parms->scale;
 	x_init[4] = parms->k0;
 	x_init[5] = parms->k1;
+	x_init[6] = parms->x0;
 
 	x = gsl_vector_view_array (x_init, num_params);
 
@@ -216,6 +224,7 @@ ProjectionLSQ::lsq(const Hills *h, ViewParams *parms,
 	if (distortion_correct) {
 		parms->k0 = gsl_vector_get(s->x, 4);
 		parms->k1 = gsl_vector_get(s->x, 5);
+		parms->x0 = gsl_vector_get(s->x, 6);
 	}
 
 	gsl_multifit_fdfsolver_free (s);
@@ -228,9 +237,9 @@ ProjectionLSQ::get_coordinates(double alph, double a_nick,
 	const ViewParams *parms, double *x, double *y) {
 
 	*x = mac_x(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale,
-		parms->k0, parms->k1, alph, a_nick); 
+		parms->k0, parms->k1, parms->x0, alph, a_nick); 
 	*y = mac_y(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale,
-		parms->k0, parms->k1, alph, a_nick); 
+		parms->k0, parms->k1, parms->x0, alph, a_nick); 
 }
 
 double
@@ -238,7 +247,7 @@ ProjectionLSQ::comp_scale(double a1, double a2, double d1, double d2) {
 	return (fabs(d1 - d2) / fabs(a1 - a2));
 }
 
-#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick
 
 double ProjectionLSQ::mac_x(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y(ARGS) {return NAN;}
@@ -248,9 +257,11 @@ double ProjectionLSQ::mac_x_dc_tilt(ARGS) {return NAN;}
 double ProjectionLSQ::mac_x_dscale(ARGS) {return NAN;}
 double ProjectionLSQ::mac_x_dk0(ARGS) {return NAN;}
 double ProjectionLSQ::mac_x_dk1(ARGS) {return NAN;}
+double ProjectionLSQ::mac_x_dx0(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dc_view(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dc_nick(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dc_tilt(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dscale(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dk0(ARGS) {return NAN;}
 double ProjectionLSQ::mac_y_dk1(ARGS) {return NAN;}
+double ProjectionLSQ::mac_y_dx0(ARGS) {return NAN;}
diff --git a/src/ProjectionRectilinear.H b/src/ProjectionRectilinear.H
index e68851e..f3f16bb 100644
--- a/src/ProjectionRectilinear.H
+++ b/src/ProjectionRectilinear.H
@@ -14,7 +14,7 @@ class ProjectionRectilinear : public ProjectionLSQ {
 
 		virtual double get_view_angle() {return 1.0471976;}; /* 60 deg */
 
-#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick
 		virtual double mac_x(ARGS);
 		virtual double mac_y(ARGS);
 		virtual double mac_x_dc_view(ARGS);
@@ -23,12 +23,14 @@ class ProjectionRectilinear : public ProjectionLSQ {
 		virtual double mac_x_dscale(ARGS);
 		virtual double mac_x_dk0(ARGS);
 		virtual double mac_x_dk1(ARGS);
+		virtual double mac_x_dx0(ARGS);
 		virtual double mac_y_dc_view(ARGS);
 		virtual double mac_y_dc_nick(ARGS);
 		virtual double mac_y_dc_tilt(ARGS);
 		virtual double mac_y_dscale(ARGS);
 		virtual double mac_y_dk0(ARGS);
 		virtual double mac_y_dk1(ARGS);
+		virtual double mac_y_dx0(ARGS);
 #undef ARGS
 
 };
diff --git a/src/ProjectionRectilinear.cxx b/src/ProjectionRectilinear.cxx
index fd64b7c..a6ec6d1 100644
--- a/src/ProjectionRectilinear.cxx
+++ b/src/ProjectionRectilinear.cxx
@@ -10,4 +10,3 @@
 #include "ProjectionRectilinear.H"
 
 #include "ProjectionRectilinear_funcs.cxx"
-
diff --git a/src/Stitch.H b/src/Stitch.H
index dcfc691..0632192 100644
--- a/src/Stitch.H
+++ b/src/Stitch.H
@@ -13,13 +13,14 @@
 #define MAX_PICS 256
 
 
+
 class Stitch {
 	private:
 		GipfelWidget *gipf[MAX_PICS];
+		int num_pics;
 		OutputImage *single_images[MAX_PICS];
 		OutputImage *merged_image;
 
-
 	public:
 		Stitch();
 
diff --git a/src/Stitch.cxx b/src/Stitch.cxx
index f21c140..0038b28 100644
--- a/src/Stitch.cxx
+++ b/src/Stitch.cxx
@@ -1,5 +1,5 @@
 //
-// Copyright 2006 Johannes Hofmann <Johannes.Hofmann@gmx.de>
+// Copyright 2007 Johannes Hofmann <Johannes.Hofmann@gmx.de>
 //
 // This software may be used and distributed according to the terms
 // of the GNU General Public License, incorporated herein by reference.
@@ -9,11 +9,17 @@
 #include <string.h>
 #include <math.h>
 
+#include <gsl/gsl_multifit.h>
+
 #include <Fl/Fl.H>
 
 #include "OutputImage.H"
 #include "Stitch.H"
 
+#define MIN(A,B) ((A)<(B)?(A):(B))
+#define MAX(A,B) ((A)>(B)?(A):(B))
+#define MAX_VALUE 65025
+
 static double pi_d = asin(1.0) * 2.0;
 static double deg2rad = pi_d / 180.0;
 
@@ -22,7 +28,9 @@ Stitch::Stitch() {
 		gipf[i] = NULL;
 		single_images[i] = NULL;
 	}
+
 	merged_image = NULL;
+	num_pics = 0;
 }
 
 Stitch::~Stitch() {
@@ -44,6 +52,8 @@ Stitch::load_image(char *file) {
 			if (gipf[i]->load_image(file) != 0) {
 				delete gipf[i];
 				gipf[i] = NULL;
+			} else {
+				num_pics++;
 			}
 			break;
 		}
@@ -85,7 +95,7 @@ Stitch::resample(GipfelWidget::sample_mode_t m,
 	view_end = view_end * deg2rad;
 
 	double step_view = (view_end - view_start) / w;
-	uchar r, g, b;
+	int r, g, b;
 	int y_off = h / 2;
 	int merged_pixel_set;
 	double radius = (double) w / (view_end -view_start);
@@ -93,24 +103,29 @@ Stitch::resample(GipfelWidget::sample_mode_t m,
 	if (merged_image) {
 		merged_image->init(w, h);
 	}
+
 	for (int i=0; i<MAX_PICS; i++) {
 		if (single_images[i]) {
 			single_images[i]->init(w, h);
 		}
 	}
 
-	for (int y=0; y<h; y++) {
+	for (int y = 0; y < h; y++) {
 		double a_nick = atan((double)(y_off - y)/radius);
+		double a_max = 0.0;
 
-		for (int x=0; x<w; x++) {
+		for (int x = 0; x < w; x++) {
 			double a_view;
 			a_view = view_start + x * step_view;
 			merged_pixel_set = 0;
-			for (int i=0; i<MAX_PICS; i++) {
-				if (gipf[i] == NULL) {
-					break;
-				} else if (gipf[i]->get_pixel(m, a_view, a_nick,
+			for (int i = 0; i < num_pics; i++) {
+				if (gipf[i]->get_pixel(m, a_view, a_nick,
 						&r, &g, &b) == 0) {
+
+					r = MAX(MIN(r, MAX_VALUE), 0);
+					g = MAX(MIN(g, MAX_VALUE), 0);
+					b = MAX(MIN(b, MAX_VALUE), 0);
+
 					if (single_images[i]) {
 						single_images[i]->set_pixel(x, r, g, b);
 					}
@@ -142,3 +157,4 @@ Stitch::resample(GipfelWidget::sample_mode_t m,
 
 	return 0;
 }
+
diff --git a/src/TIFFOutputImage.H b/src/TIFFOutputImage.H
index 820a281..9b0f83b 100644
--- a/src/TIFFOutputImage.H
+++ b/src/TIFFOutputImage.H
@@ -12,21 +12,22 @@
 
 #include "OutputImage.H"
 
-class TIFFOutputImage : OutputImage {
+class TIFFOutputImage : public OutputImage {
 	private:
+		int bitspersample;
 		unsigned char *row;
 		char *file;
 		TIFF *tiff;
 
 	public:
-		TIFFOutputImage(const char *file);
+		TIFFOutputImage(const char *file, int b = 8);
 
 		~TIFFOutputImage();
 
 	protected:
 		int init_internal(int w, int h);
 
-		int set_pixel_internal(int x, char r, char g, char b);	
+		int set_pixel_internal(int x, int r, int g, int b);	
 
 		int next_line_internal();
 
diff --git a/src/TIFFOutputImage.cxx b/src/TIFFOutputImage.cxx
index 397caac..b8effed 100644
--- a/src/TIFFOutputImage.cxx
+++ b/src/TIFFOutputImage.cxx
@@ -10,7 +10,8 @@
 
 #include "TIFFOutputImage.H"
 
-TIFFOutputImage::TIFFOutputImage(const char *f) {
+TIFFOutputImage::TIFFOutputImage(const char *f, int b) {
+	bitspersample = (b==16)?16:8;
 	file = strdup(f);
 	tiff = NULL;
 	row = NULL;
@@ -32,12 +33,12 @@ TIFFOutputImage::init_internal(int w1, int h1) {
 		row = NULL;
 	}
 
-	row = (unsigned char*) malloc(sizeof(char) * 4 * w1);
+	row = (unsigned char*) malloc(sizeof(char) * (bitspersample / 8) * 4 * w1);
 	if (!row) {
 		perror("malloc");
 		return 1;
 	}
-	memset(row, 0, sizeof(char) * 4 * w1);
+	memset(row, 0, sizeof(char) * (bitspersample / 8) * 4 * w1);
 
 	if (tiff) {
 		TIFFClose(tiff);
@@ -54,25 +55,33 @@ TIFFOutputImage::init_internal(int w1, int h1) {
 	TIFFSetField(tiff, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
 	TIFFSetField(tiff, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
 	TIFFSetField(tiff, TIFFTAG_ROWSPERSTRIP, 1);
-	TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, 8);
+	TIFFSetField(tiff, TIFFTAG_BITSPERSAMPLE, bitspersample);
 	TIFFSetField(tiff, TIFFTAG_SAMPLESPERPIXEL, 4);
 
 	return 0;
 }
 
 int
-TIFFOutputImage::set_pixel_internal(int x, char r, char g, char b) {
-	row[x*4+0] = r;
-	row[x*4+1] = g;
-	row[x*4+2] = b;
-	row[x*4+3] = 255;
+TIFFOutputImage::set_pixel_internal(int x, int r, int g, int b) {
+	if (bitspersample == 8) {	
+		row[x*4+0] = (unsigned char) (r / 255);
+		row[x*4+1] = (unsigned char) (g / 255);
+		row[x*4+2] = (unsigned char) (b / 255);
+		row[x*4+3] = 255;
+	} else if (bitspersample == 16) {
+		unsigned short *row16 = (unsigned short*) row;
+		row16[x*4+0] = (unsigned short) r;
+		row16[x*4+1] = (unsigned short) g;
+		row16[x*4+2] = (unsigned short) b;
+		row16[x*4+3] = 65025;
+	}
 
 	return 0;
 }
 
 int
 TIFFOutputImage::next_line_internal() {
-	TIFFWriteEncodedStrip(tiff, line -1 , row, W * 4);
+	TIFFWriteEncodedStrip(tiff, line - 1 , row, W * (bitspersample / 8) * 4);
 
 	memset(row, 0, sizeof(char) * 4 * W);
 	return 0;
diff --git a/src/ViewParams.H b/src/ViewParams.H
index 2048fe2..44462c0 100644
--- a/src/ViewParams.H
+++ b/src/ViewParams.H
@@ -15,6 +15,7 @@ class ViewParams {
 		double a_tilt;
 		double k0;
 		double k1;
+		double x0;
 };
 
 #endif
diff --git a/src/gipfel.cxx b/src/gipfel.cxx
index 108c3ad..6df400f 100644
--- a/src/gipfel.cxx
+++ b/src/gipfel.cxx
@@ -46,19 +46,20 @@ Fl_Window *control_win, *view_win;
 Fl_Dial *s_center = NULL;
 Fl_Slider *s_nick, *s_focal_length, *s_tilt, *s_height_dist, *s_track_width;
 Fl_Value_Input *i_view_lat, *i_view_long, *i_view_height;
-Fl_Value_Input *i_distortion_k0, *i_distortion_k1;
+Fl_Value_Input *i_distortion_k0, *i_distortion_k1, *i_distortion_x0;
 Fl_Box *b_viewpoint;
 Fl_Menu_Bar *mb;
 
-#define STITCH_PREVIEW 1
-#define STITCH_JPEG    2
-#define STITCH_TIFF    4
+#define STITCH_PREVIEW           1
+#define STITCH_JPEG              2
+#define STITCH_TIFF              4
 
-static int stitch(GipfelWidget::sample_mode_t m ,int stitch_w, int stitch_h,
+static int stitch(GipfelWidget::sample_mode_t m , int b_16,
+	int stitch_w, int stitch_h,
 	double from, double to, int type, const char *path, int argc, char **argv);
 
 void set_values() {
-	double k0 = 0.0, k1 = 0.0;
+	double k0 = 0.0, k1 = 0.0, x0 = 0.0;
 
 	s_center->value(gipf->get_center_angle());
 	s_nick->value(gipf->get_nick_angle());
@@ -77,9 +78,10 @@ void set_values() {
 		mb->mode(8, FL_MENU_RADIO);
 	}
 
-	gipf->get_distortion_params(&k0, &k1);
+	gipf->get_distortion_params(&k0, &k1, &x0);
 	i_distortion_k0->value(k0);
 	i_distortion_k1->value(k1);
+	i_distortion_x0->value(x0);
 }
 
 void quit_cb() {
@@ -174,14 +176,14 @@ void comp_cb(Fl_Widget *, void *) {
 void save_distortion_cb(Fl_Widget *, void *) {
 	char buf[1024];
 	const char * prof_name;
-	double k0, k1;
+	double k0, k1, x0;
 
 	buf[0] = '\0';
+	gipf->get_distortion_params(&k0, &k1, &x0);
 
-	gipf->get_distortion_params(&k0, &k1);
 	gipf->get_distortion_profile_name(buf, sizeof(buf));
-	prof_name = fl_input("Save Distortion Profile (k0=%f, k1=%f)",
-		buf, k0, k1);
+	prof_name = fl_input("Save Distortion Profile (k0=%f, k1=%f, x0=%f)",
+		buf, k0, k1, x0);
 
 	if (prof_name == NULL) {
 		return;
@@ -215,8 +217,10 @@ void load_distortion_cb(Fl_Widget *, void *) {
 }
 
 void distortion_cb(Fl_Value_Input*, void*) {
-	gipf->set_distortion_params(i_distortion_k0->value(),
-		i_distortion_k1->value());
+	gipf->set_distortion_params(
+		i_distortion_k0->value(),
+		i_distortion_k1->value(),
+		i_distortion_x0->value());
 }
 
 void about_cb() {
@@ -260,6 +264,7 @@ void usage() {
 		"   -d <datafile>   Use <datafile> for GPS data.\n"
 		"   -u <k0>,<k1>    Use distortion correction values k0,k1.\n"
 		"   -s              Stitch mode.\n"
+		"   -4              Create 16bit output (only with TIFF stitching).\n"
 		"   -r <from>,<to>  Stitch range in degrees (e.g. 100.0,200.0).\n"
 		"   -b              Use bilinear interpolation for stitching.\n"
 		"   -w <width>      Width of result image.\n"
@@ -362,23 +367,31 @@ create_control_window() {
 	i_view_height->when(FL_WHEN_ENTER_KEY);
 	i_view_height->callback((Fl_Callback*)view_height_cb);
 
-	i_distortion_k0 = new Fl_Value_Input(235, 220, 80, 20, "Distortion (k0)");
+	i_distortion_k0 = new Fl_Value_Input(250, 200, 80, 20, "k0");
 	i_distortion_k0->labelsize(10);
 	i_distortion_k0->textsize(10);
-	i_distortion_k0->align(FL_ALIGN_TOP);
+	i_distortion_k0->align(FL_ALIGN_LEFT);
 	i_distortion_k0->when(FL_WHEN_ENTER_KEY);
 	i_distortion_k0->callback((Fl_Callback*)distortion_cb);
 
-	i_distortion_k1 = new Fl_Value_Input(315, 220, 80,  20, "Distortion (k1)");
+	i_distortion_k1 = new Fl_Value_Input(250, 225, 80,  20, "k1");
 	i_distortion_k1->labelsize(10);
 	i_distortion_k1->textsize(10);
-	i_distortion_k1->align(FL_ALIGN_TOP);
+	i_distortion_k1->align(FL_ALIGN_LEFT);
 	i_distortion_k1->when(FL_WHEN_ENTER_KEY);
 	i_distortion_k1->callback((Fl_Callback*)distortion_cb);
 
+	i_distortion_x0 = new Fl_Value_Input(250, 250, 80,  20, "x0");
+	i_distortion_x0->labelsize(10);
+	i_distortion_x0->textsize(10);
+	i_distortion_x0->align(FL_ALIGN_LEFT);
+	i_distortion_x0->when(FL_WHEN_ENTER_KEY);
+	i_distortion_x0->callback((Fl_Callback*)distortion_cb);
+
+
 
 	// Buttons
-	Fl_Button *b = new Fl_Button(260, 260, 100, 20, "comp");
+	Fl_Button *b = new Fl_Button(280, 280, 100, 20, "comp");
 	b->color(FL_RED);
 	b->tooltip("compute view parameter from given mountains");
 	b->callback(comp_cb);
@@ -393,15 +406,15 @@ int main(int argc, char** argv) {
 	int err, my_argc;
 	int stitch_flag = 0, stitch_w = 2000, stitch_h = 500;
 	int jpeg_flag = 0, tiff_flag = 0, distortion_flag = 0;
-	int bilinear_flag = 0;
+	int bicubic_flag = 0, b_16_flag = 0;
 	double stitch_from = 0.0, stitch_to = 380.0;
-	double dist_k0 = 0.0, dist_k1 = 0.0;
+	double dist_k0 = 0.0, dist_k1 = 0.0, dist_x0 = 0.0;
 	char *outpath = "/tmp";
 	Fl_Scroll *scroll;
 
 
 	err = 0;
-	while ((c = getopt(argc, argv, ":?d:v:sw:h:j:t:u:br:")) != EOF) {
+	while ((c = getopt(argc, argv, ":?d:v:sw:h:j:t:u:br:4")) != EOF) {
 		switch (c) {  
 			case '?':
 				usage();
@@ -416,6 +429,9 @@ int main(int argc, char** argv) {
 			case 's':
 				stitch_flag++;
 				break;
+			case '4':
+				b_16_flag++;
+				break;
 			case 'r':
 				stitch_flag++;
 				if (optarg && strcmp(optarg, ":")) {
@@ -426,11 +442,17 @@ int main(int argc, char** argv) {
 				}
 				break;
 			case 'u':
+				char *c;
 				distortion_flag++;
 				if (optarg && strcmp(optarg, ":")) {
 					dist_k0 = atof(optarg);
-					if (strchr(optarg, ',')) {
+					c = strchr(optarg, ',');
+					if (c) {
 						dist_k1 = atof(strchr(optarg, ',') + 1);
+						c = strchr(c + 1, ',');
+						if (c) {
+							dist_x0 = atof(strchr(optarg, ',') + 1);
+						}
 					}
 				}
 				break;
@@ -449,7 +471,7 @@ int main(int argc, char** argv) {
 				stitch_h = atoi(optarg);
 				break;
 			case 'b':
-				bilinear_flag++;
+				bicubic_flag++;
 				break;
 			default:
 				err++;
@@ -470,14 +492,17 @@ int main(int argc, char** argv) {
 	}
 
 	if (stitch_flag) {
-		int type = STITCH_PREVIEW;
+		int type = 0;
 		if (jpeg_flag) {
 			type = STITCH_JPEG;
 		} else if (tiff_flag) {
 			type = STITCH_TIFF;
+		} else {
+			type = STITCH_PREVIEW;
 		}
 
-		stitch(bilinear_flag?GipfelWidget::BILINEAR:GipfelWidget::NEAREST,
+		stitch(bicubic_flag?GipfelWidget::BICUBIC:GipfelWidget::NEAREST,
+			b_16_flag,
 			stitch_w, stitch_h, stitch_from, stitch_to,
 			type, outpath, my_argc, my_argv);
 
@@ -508,7 +533,7 @@ int main(int argc, char** argv) {
 	}
 
 	if (distortion_flag) {
-		gipf->set_distortion_params(dist_k0, dist_k1);
+		gipf->set_distortion_params(dist_k0, dist_k1, dist_x0);
 	}
 
 	view_win->size(gipf->w(), gipf->h());
@@ -536,7 +561,7 @@ int main(int argc, char** argv) {
 }
 
 static int
-stitch(GipfelWidget::sample_mode_t m,
+stitch(GipfelWidget::sample_mode_t m, int b_16,
 	int stitch_w, int stitch_h, double from, double to,
 	int type, const char *path, int argc, char **argv) {
 
@@ -548,12 +573,13 @@ stitch(GipfelWidget::sample_mode_t m,
 		st->load_image(argv[i]);
 	}
 
-	if (type == STITCH_JPEG) {
 
-		st->set_output((OutputImage*) new JPEGOutputImage(path, 90));
+	if (type & STITCH_JPEG) {
+
+		st->set_output(new JPEGOutputImage(path, 90));
 		st->resample(m, stitch_w, stitch_h, from, to);
 
-	} else if (type == STITCH_TIFF) {
+	} else if (type & STITCH_TIFF) {
 
 		for (int i=0; i<argc; i++) {
 			char buf[1024];
@@ -566,7 +592,7 @@ stitch(GipfelWidget::sample_mode_t m,
 			}
 			strncat(buf, ".tiff", sizeof(buf));
 
-			st->set_output(argv[i], (OutputImage*) new TIFFOutputImage(buf));
+			st->set_output(argv[i], new TIFFOutputImage(buf, b_16?16:8));
 		}
 
 		st->resample(m, stitch_w, stitch_h, from, to);
@@ -579,8 +605,10 @@ stitch(GipfelWidget::sample_mode_t m,
 
 		win->resizable(scroll);
 		win->show(0, argv); 
-		st->set_output((OutputImage*) img);
+		st->set_output(img);
+
 		st->resample(m, stitch_w, stitch_h, from, to);
+
 		img->redraw();
 		Fl::run();
 	}
diff --git a/src/lsq_cylindrical.mac b/src/lsq_cylindrical.mac
index 0bf9499..a88a2e2 100644
--- a/src/lsq_cylindrical.mac
+++ b/src/lsq_cylindrical.mac
@@ -18,7 +18,7 @@ load("expr2c.mac")$
 x_expand : trigexpand(x_scale)$
 y_expand : trigexpand(y_scale)$
 
-args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick"$
+args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick"$
 
 expr2c("ProjectionCylindrical::mac_x", args, x_expand)$
 expr2c("ProjectionCylindrical::mac_y", args, y_expand)$
diff --git a/src/lsq_rectilinear.mac b/src/lsq_rectilinear.mac
index 69980c0..5bd7ac9 100644
--- a/src/lsq_rectilinear.mac
+++ b/src/lsq_rectilinear.mac
@@ -4,7 +4,7 @@
 
 x : tan(m_view - c_view)$
 y : tan(c_nick - m_nick)$ 
-x_rot : y * sin(c_tilt) + x * cos(c_tilt)$
+x_rot : y * sin(c_tilt) + x * cos(c_tilt)+x0$
 y_rot : y * cos(c_tilt) - x * sin(c_tilt)$
 d : x_rot ^ 2 + y_rot ^ 2$
 dist_fact : d ^ 2 * k1 + d * k0$
@@ -20,7 +20,7 @@ load("expr2c.mac")$
 x_expand : trigexpand(x_dist)$
 y_expand : trigexpand(y_dist)$
 
-args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick"$
+args: "double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double x0, double m_view, double m_nick"$
 
 expr2c("ProjectionRectilinear::mac_x", args, x_expand)$
 expr2c("ProjectionRectilinear::mac_y", args, y_expand)$
@@ -31,9 +31,11 @@ expr2c("ProjectionRectilinear::mac_x_dc_tilt", args, diff(x_expand, c_tilt))$
 expr2c("ProjectionRectilinear::mac_x_dscale", args, diff(x_expand, scale))$
 expr2c("ProjectionRectilinear::mac_x_dk0", args, diff(x_expand, k0))$
 expr2c("ProjectionRectilinear::mac_x_dk1", args, diff(x_expand, k1))$
+expr2c("ProjectionRectilinear::mac_x_dx0", args, diff(x_expand, x0))$
 expr2c("ProjectionRectilinear::mac_y_dc_view", args, diff(y_expand, c_view))$
 expr2c("ProjectionRectilinear::mac_y_dc_nick", args, diff(y_expand, c_nick))$
 expr2c("ProjectionRectilinear::mac_y_dc_tilt", args, diff(y_expand, c_tilt))$
 expr2c("ProjectionRectilinear::mac_y_dscale", args, diff(y_expand, scale))$
 expr2c("ProjectionRectilinear::mac_y_dk0", args, diff(y_expand, k0))$
 expr2c("ProjectionRectilinear::mac_y_dk1", args, diff(y_expand, k1))$
+expr2c("ProjectionRectilinear::mac_y_dx0", args, diff(y_expand, x0))$