use LSQ for all projection types

19 files changed, 273 insertions, 161 deletions

diff --git a/src/GipfelWidget.H b/src/GipfelWidget.H
index d504de0..b6bb901 100644
--- a/src/GipfelWidget.H
+++ b/src/GipfelWidget.H
@@ -104,9 +104,9 @@ class GipfelWidget : public Fl_Widget {
 
 		void set_track_width(double w);
 
-		Projection::Projection_t get_projection();
+		ProjectionLSQ::Projection_t get_projection();
 
-		void set_projection(Projection::Projection_t p);
+		void set_projection(ProjectionLSQ::Projection_t p);
 
 		void set_distortion_params(double k0, double k1);
 
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx
index e8ca448..5f75693 100644
--- a/src/GipfelWidget.cxx
+++ b/src/GipfelWidget.cxx
@@ -98,7 +98,7 @@ GipfelWidget::load_image(char *file) {
 	set_center_angle(md->get_direction());
 	set_nick_angle(md->get_nick());
 	set_tilt_angle(md->get_tilt());
-	set_projection((Projection::Projection_t) md->get_projection_type());
+	set_projection((ProjectionLSQ::Projection_t) md->get_projection_type());
 	set_focal_length_35mm(md->get_focal_length_35mm());
 	md->get_distortion_params(&pan->parms.k0, &pan->parms.k1);
 	fprintf(stderr, "%f %f\n", pan->parms.k0, pan->parms.k1);
@@ -470,7 +470,7 @@ GipfelWidget::set_focal_length_35mm(double s) {
 }
 
 void
-GipfelWidget::set_projection(Projection::Projection_t p) {
+GipfelWidget::set_projection(ProjectionLSQ::Projection_t p) {
 	pan->set_projection(p);
 	set_labels(pan->get_visible_mountains());
 	redraw();
@@ -482,7 +482,7 @@ GipfelWidget::set_distortion_params(double k0, double k1) {
 	redraw();
 }
 
-Projection::Projection_t
+ProjectionLSQ::Projection_t
 GipfelWidget::get_projection() {
 	return pan->get_projection();
 }
diff --git a/src/Makefile.am b/src/Makefile.am
index 36c5b80..e55ebf4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -5,9 +5,9 @@ gipfel_SOURCES = \
 	util.c \
 	GipfelWidget.cxx \
 	Panorama.cxx \
-	Projection.cxx \
-	ProjectionTangentialLSQ.cxx \
-	ProjectionSphaeric.cxx \
+	ProjectionLSQ.cxx \
+	ProjectionRectilinear.cxx \
+	ProjectionCylindrical.cxx \
 	Hill.cxx \
 	Fl_Value_Dial.cxx \
 	Fl_Search_Chooser.cxx \
@@ -22,9 +22,9 @@ gipfel_SOURCES = \
 noinst_HEADERS = \
 	GipfelWidget.H \
 	Panorama.H \
-	Projection.H \
-	ProjectionTangentialLSQ.H \
-	ProjectionSphaeric.H \
+	ProjectionLSQ.H \
+	ProjectionRectilinear.H \
+	ProjectionCylindrical.H \
 	Hill.H \
 	ViewParams.H \
 	Fl_Value_Dial.H \
diff --git a/src/Makefile.lsq_funcs b/src/Makefile.lsq_funcs
index f117e1e..8f92113 100644
--- a/src/Makefile.lsq_funcs
+++ b/src/Makefile.lsq_funcs
@@ -1,4 +1,7 @@
-all: lsq_funcs.c
+all: ProjectionRectilinear_funcs.cxx ProjectionCylindrical_funcs.cxx
 
-lsq_funcs.c: lsq_funcs.mac
-	maxima -b lsq_funcs.mac | grep "^static" > lsq_funcs.c
+ProjectionRectilinear_funcs.cxx: lsq_rectilinear.mac
+	maxima -b lsq_rectilinear.mac | grep "^double" > ProjectionRectilinear_funcs.cxx
+
+ProjectionCylindrical_funcs.cxx: lsq_cylindrical.mac
+	maxima -b lsq_cylindrical.mac | grep "^double" > ProjectionCylindrical_funcs.cxx
diff --git a/src/Panorama.H b/src/Panorama.H
index ca000f3..5dd670f 100644
--- a/src/Panorama.H
+++ b/src/Panorama.H
@@ -8,7 +8,7 @@
 #define PANORAMA_H
 
 #include "Hill.H"
-#include "Projection.H"
+#include "ProjectionLSQ.H"
 #include "ViewParams.H"
 
 
@@ -23,8 +23,8 @@ class Panorama {
 		Hills *mountains;
 		Hills *close_mountains;
 		Hills *visible_mountains;
-		Projection::Projection *proj;
-		Projection::Projection_t projection_type;
+		ProjectionLSQ *proj;
+		ProjectionLSQ::Projection_t projection_type;
 
 		Hill * get_pos(const char *name);
 
@@ -125,9 +125,9 @@ class Panorama {
 
 		int comp_params(Hills *h);
 
-		Projection::Projection_t get_projection();
+		ProjectionLSQ::Projection_t get_projection();
 
-		void set_projection(Projection::Projection_t p);
+		void set_projection(ProjectionLSQ::Projection_t p);
 
 		void set_distortion_params(double k0, double k1);
 
diff --git a/src/Panorama.cxx b/src/Panorama.cxx
index 7ec3ccd..784c6ed 100644
--- a/src/Panorama.cxx
+++ b/src/Panorama.cxx
@@ -10,8 +10,8 @@
 #include <math.h>
 
 #include "Panorama.H"
-#include "ProjectionTangentialLSQ.H"
-#include "ProjectionSphaeric.H"
+#include "ProjectionRectilinear.H"
+#include "ProjectionCylindrical.H"
 
 #define EARTH_RADIUS 6371010.0
 
@@ -32,7 +32,7 @@ Panorama::Panorama() {
 	view_lam = 0.0;
 	view_height = 0.0;
 	proj = NULL;
-	set_projection(Projection::TANGENTIAL);
+	set_projection(ProjectionLSQ::RECTILINEAR);
 }
 
 Panorama::~Panorama() {
@@ -220,7 +220,7 @@ Panorama::set_view_height(double v) {
 }
 
 void
-Panorama::set_projection(Projection::Projection_t p) {
+Panorama::set_projection(ProjectionLSQ::Projection_t p) {
 	projection_type = p;
 
 	if (proj) {
@@ -228,12 +228,12 @@ Panorama::set_projection(Projection::Projection_t p) {
 	}
 
 	switch(projection_type) {
-		case Projection::TANGENTIAL:
-			proj = new ProjectionTangentialLSQ();
+		case ProjectionLSQ::RECTILINEAR:
+			proj = new ProjectionRectilinear();
 			view_angle = pi_d / 3.0;
 			break;
-		case Projection::SPHAERIC:
-			proj = new ProjectionSphaeric();
+		case ProjectionLSQ::CYLINDRICAL:
+			proj = new ProjectionCylindrical();
 			view_angle = pi_d * 2.0;
 			break;
 	}
@@ -285,7 +285,7 @@ Panorama::get_view_height() {
 	return view_height;
 }
 
-Projection::Projection_t
+ProjectionLSQ::Projection_t
 Panorama::get_projection() {
 	return projection_type;
 }
diff --git a/src/Projection.H b/src/Projection.H
deleted file mode 100644
index 2a953b2..0000000
--- a/src/Projection.H
+++ /dev/null
@@ -1,30 +0,0 @@
-//
-// Copyright 2006 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.
-
-#ifndef PROJECTION_H
-#define PROJECTION_H
-
-#include "Hill.H"
-#include "ViewParams.H"
-
-class Projection {
-	protected:
-		double pi_d;
-
-	public:
-		typedef enum {
-			TANGENTIAL = 0,
-			SPHAERIC   = 1
-		} Projection_t;
-
-		Projection();
-
-		virtual void get_coordinates(double a_view, double a_nick,
-			const ViewParams *parms, double *x, double *y);
-
-		virtual int comp_params(const Hills *h, ViewParams *parms);
-};
-#endif
diff --git a/src/Projection.cxx b/src/Projection.cxx
deleted file mode 100644
index bb577a0..0000000
--- a/src/Projection.cxx
+++ /dev/null
@@ -1,26 +0,0 @@
-//
-// Copyright 2006 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.
-
-#include <stdio.h>
-#include <math.h>
-
-#include "Projection.H"
-
-Projection::Projection() {
-	pi_d = asin(1.0) * 2.0;
-};
-
-void
-Projection::get_coordinates(double a_view, double a_nick,
-	const ViewParams *parms, double *x, double *y) {
-	fprintf(stderr, "Error: Projection::set_coordinates()\n");
-}
-
-int 
-Projection::comp_params(const Hills *h, ViewParams *parms) {
-	fprintf(stderr, "Error: Projection::comp_params()\n");
-	return 1;
-}
diff --git a/src/ProjectionCylindrical.H b/src/ProjectionCylindrical.H
new file mode 100644
index 0000000..e64a911
--- /dev/null
+++ b/src/ProjectionCylindrical.H
@@ -0,0 +1,35 @@
+//
+// Copyright 2006 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.
+
+#ifndef PROJECTIONCYLINDRICAL_H
+#define PROJECTIONCYLINDRICAL_H
+
+#include "ProjectionLSQ.H"
+
+class ProjectionCylindrical : public ProjectionLSQ {
+	public:
+
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+
+		virtual double mac_x(ARGS);
+		virtual double mac_y(ARGS);
+        virtual double mac_x_dc_view(ARGS);
+		virtual double mac_x_dc_nick(ARGS);
+		virtual double mac_x_dc_tilt(ARGS);
+		virtual double mac_x_dscale(ARGS);
+		virtual double mac_x_dk0(ARGS);
+		virtual double mac_x_dk1(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);
+
+#undef ARGS
+
+};
+#endif
diff --git a/src/ProjectionCylindrical.cxx b/src/ProjectionCylindrical.cxx
new file mode 100644
index 0000000..d41e5f8
--- /dev/null
+++ b/src/ProjectionCylindrical.cxx
@@ -0,0 +1,12 @@
+//
+// Copyright 2006 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.
+
+#include <stdio.h>
+#include <math.h>
+
+#include "ProjectionCylindrical.H"
+
+#include "ProjectionCylindrical_funcs.cxx"
diff --git a/src/ProjectionLSQ.H b/src/ProjectionLSQ.H
new file mode 100644
index 0000000..0834108
--- /dev/null
+++ b/src/ProjectionLSQ.H
@@ -0,0 +1,54 @@
+//
+// Copyright 2006 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.
+
+#ifndef PROJECTIONLSQ_H
+#define PROJECTIONLSQ_H
+
+#include "Hill.H"
+#include "ViewParams.H"
+
+class ProjectionLSQ {
+	private:
+		double comp_scale(double alph_a, double alph_b, double d1, double d2);
+
+		int lsq(const Hills *m, ViewParams *parms, int distortion_correct);
+
+	protected:
+		double sec(double a);
+
+	public:
+
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+
+		virtual double mac_x(ARGS);
+		virtual double mac_y(ARGS);
+		virtual double mac_x_dc_view(ARGS);
+		virtual double mac_x_dc_nick(ARGS);
+		virtual double mac_x_dc_tilt(ARGS);
+		virtual double mac_x_dscale(ARGS);
+		virtual double mac_x_dk0(ARGS);
+		virtual double mac_x_dk1(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);
+
+#undef ARGS
+
+
+		typedef enum {
+			RECTILINEAR = 0,
+			CYLINDRICAL = 1
+		} Projection_t;
+
+		void get_coordinates(double a_view, double a_nick,
+			const ViewParams *parms, double *x, double *y);
+
+		int comp_params(const Hills *h, ViewParams *parms);
+};
+#endif
diff --git a/src/ProjectionTangentialLSQ.cxx b/src/ProjectionLSQ.cxx
index 76ec2fa..3b953c8 100644
--- a/src/ProjectionTangentialLSQ.cxx
+++ b/src/ProjectionLSQ.cxx
@@ -15,17 +15,16 @@
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_multifit_nlin.h>
 
-#include "ProjectionTangentialLSQ.H"
+#include "ProjectionLSQ.H"
 
 
-static double sec(double a) {
+double
+ProjectionLSQ::sec(double a) {
 	return 1.0 / cos(a);
 }
 
-#include "lsq_funcs.c"
-
 int
-ProjectionTangentialLSQ::comp_params(const Hills *h, ViewParams *parms) {
+ProjectionLSQ::comp_params(const Hills *h, ViewParams *parms) {
 	const Hill *tmp, *m1, *m2;
 	double scale_tmp;
 	int distortion_correct = 0;
@@ -73,12 +72,13 @@ ProjectionTangentialLSQ::comp_params(const Hills *h, ViewParams *parms) {
 }
 
 struct data {
+	ProjectionLSQ *p;
 	int distortion_correct;
 	const Hills *h;
 	const ViewParams *old_params;
 };
 
-#define CALL(A) 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, m->alph, m->a_nick) 
 
 static int
 lsq_f (const gsl_vector * x, void *data, gsl_vector * f) {
@@ -162,7 +162,7 @@ lsq_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) {
 }
 
 int
-ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms,
+ProjectionLSQ::lsq(const Hills *h, ViewParams *parms,
 	int distortion_correct) {
 
 	const gsl_multifit_fdfsolver_type *T;
@@ -174,6 +174,7 @@ ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms,
 	int status;
 	int num_params = distortion_correct?6:4;
 
+	dat.p = this;
 	dat.distortion_correct = distortion_correct;
 	dat.h = h;
 	dat.old_params = parms;
@@ -225,7 +226,7 @@ ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms,
 }
 
 void 
-ProjectionTangentialLSQ::get_coordinates(double alph, double a_nick,
+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,
@@ -235,7 +236,7 @@ ProjectionTangentialLSQ::get_coordinates(double alph, double a_nick,
 }
 
 double
-ProjectionTangentialLSQ::comp_scale(double a1, double a2, double d1, double d2) {
+ProjectionLSQ::comp_scale(double a1, double a2, double d1, double d2) {
 	double sign1 = 1.0;
 	double sc, tan_a1, tan_a2;
 
diff --git a/src/ProjectionRectilinear.H b/src/ProjectionRectilinear.H
new file mode 100644
index 0000000..1fc7417
--- /dev/null
+++ b/src/ProjectionRectilinear.H
@@ -0,0 +1,35 @@
+//
+// Copyright 2006 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.
+
+#ifndef PROJECTIONRECTILINEAR_H
+#define PROJECTIONRECTILINEAR_H
+
+#include "ProjectionLSQ.H"
+
+class ProjectionRectilinear : public ProjectionLSQ {
+	public:
+
+#define ARGS double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick
+
+		virtual double mac_x(ARGS);
+		virtual double mac_y(ARGS);
+		virtual double mac_x_dc_view(ARGS);
+		virtual double mac_x_dc_nick(ARGS);
+		virtual double mac_x_dc_tilt(ARGS);
+		virtual double mac_x_dscale(ARGS);
+		virtual double mac_x_dk0(ARGS);
+		virtual double mac_x_dk1(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);
+
+#undef ARGS
+
+};
+#endif
diff --git a/src/ProjectionRectilinear.cxx b/src/ProjectionRectilinear.cxx
new file mode 100644
index 0000000..fd64b7c
--- /dev/null
+++ b/src/ProjectionRectilinear.cxx
@@ -0,0 +1,13 @@
+//
+// Copyright 2006 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.
+
+#include <stdio.h>
+#include <math.h>
+
+#include "ProjectionRectilinear.H"
+
+#include "ProjectionRectilinear_funcs.cxx"
+
diff --git a/src/ProjectionTangentialLSQ.H b/src/ProjectionTangentialLSQ.H
deleted file mode 100644
index 4e3f19b..0000000
--- a/src/ProjectionTangentialLSQ.H
+++ /dev/null
@@ -1,25 +0,0 @@
-//
-// Copyright 2006 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.
-
-#ifndef PROJECTIONTANGENTIALLSQ_H
-#define PROJECTIONTANGENTIALLSQ_H
-
-#include "Hill.H"
-#include "Projection.H"
-
-class ProjectionTangentialLSQ : public Projection {
-	private:
-		double comp_scale(double alph_a, double alph_b, double d1, double d2);
-
-		int lsq(const Hills *m, ViewParams *parms, int distortion_correct);
-
-	public:
-		void get_coordinates(double a_view, double a_nick,
-			const ViewParams *parms, double *x, double *y);
-
-		int comp_params(const Hills *h, ViewParams *parms);
-};
-#endif
diff --git a/src/gipfel.cxx b/src/gipfel.cxx
index 89a5f26..491020d 100644
--- a/src/gipfel.cxx
+++ b/src/gipfel.cxx
@@ -65,7 +65,7 @@ void set_values() {
 	i_view_long->value(gipf->get_view_long());
 	i_view_height->value(gipf->get_view_height());
 	b_viewpoint->label(gipf->get_viewpoint());
-	if (gipf->get_projection() == Projection::TANGENTIAL) {
+	if (gipf->get_projection() == ProjectionLSQ::RECTILINEAR) {
 		mb->mode(8, FL_MENU_RADIO|FL_MENU_VALUE);
 		mb->mode(9, FL_MENU_RADIO);
 	} else {
@@ -148,9 +148,9 @@ void viewpoint_cb(Fl_Value_Input* o, void*) {
 
 void proj_cb(Fl_Value_Input* o, void*d) {
 	if(d == NULL) {
-		gipf->set_projection(Projection::TANGENTIAL);
+		gipf->set_projection(ProjectionLSQ::RECTILINEAR);
 	} else {
-		gipf->set_projection(Projection::SPHAERIC);
+		gipf->set_projection(ProjectionLSQ::CYLINDRICAL);
 	}
 }
 
diff --git a/src/lsq_cylindrical.mac b/src/lsq_cylindrical.mac
new file mode 100644
index 0000000..9a1d61a
--- /dev/null
+++ b/src/lsq_cylindrical.mac
@@ -0,0 +1,40 @@
+/*
+ * This is the basic pinhole projection model with distortion correction
+ */
+
+x : tan(m_view - c_view)$
+y : tan(c_nick - m_nick)$ 
+x_rot : y * sin(c_tilt) + x * cos(c_tilt)$
+y_rot : y * cos(c_tilt) - x * sin(c_tilt)$
+d : x_rot ^ 2 + y_rot ^ 2$
+dist_fact : d ^ 2 * k1 + d * k0$
+x_dist : x_rot * (1 + dist_fact) * scale$
+y_dist : y_rot * (1 + dist_fact) * scale$
+
+/*
+ * Some mangling for C code generation
+ */
+
+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)"$
+
+printfunc(name, expression) := sprint("double ProjectionCylindrical::", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
+")$
+
+printfunc("mac_x", x_expand)$
+printfunc("mac_y", y_expand)$
+
+printfunc("mac_x_dc_view", diff(x_expand, c_view))$
+printfunc("mac_x_dc_nick", diff(x_expand, c_nick))$
+printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt))$
+printfunc("mac_x_dscale", diff(x_expand, scale))$
+printfunc("mac_x_dk0", diff(x_expand, k0))$
+printfunc("mac_x_dk1", diff(x_expand, k1))$
+printfunc("mac_y_dc_view", diff(y_expand, c_view))$
+printfunc("mac_y_dc_nick", diff(y_expand, c_nick))$
+printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt))$
+printfunc("mac_y_dscale", diff(y_expand, scale))$
+printfunc("mac_y_dk0", diff(y_expand, k0))$
+printfunc("mac_y_dk1", diff(y_expand, k1))$
diff --git a/src/lsq_funcs.mac b/src/lsq_funcs.mac
deleted file mode 100644
index de71198..0000000
--- a/src/lsq_funcs.mac
+++ /dev/null
@@ -1,40 +0,0 @@
-/*
- * This is the basic pinhole projection model with distortion correction
- */
-
-x : tan(m_view - c_view);
-y : tan(c_nick - m_nick); 
-x_rot : y * sin(c_tilt) + x * cos(c_tilt);
-y_rot : y * cos(c_tilt) - x * sin(c_tilt);
-d : x_rot ^ 2 + y_rot ^ 2;
-dist_fact : d ^ 2 * k1 + d * k0;
-x_dist : x_rot * (1 + dist_fact) * scale;
-y_dist : y_rot * (1 + dist_fact) * scale;
-
-/*
- * Some mangling for C code generation
- */
-
-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)";
-
-printfunc(name, expression) := sprint("static double", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
-");
-
-printfunc("mac_x", x_expand);
-printfunc("mac_y", y_expand);
-
-printfunc("mac_x_dc_view", diff(x_expand, c_view));
-printfunc("mac_x_dc_nick", diff(x_expand, c_nick));
-printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt));
-printfunc("mac_x_dscale", diff(x_expand, scale));
-printfunc("mac_x_dk0", diff(x_expand, k0));
-printfunc("mac_x_dk1", diff(x_expand, k1));
-printfunc("mac_y_dc_view", diff(y_expand, c_view));
-printfunc("mac_y_dc_nick", diff(y_expand, c_nick));
-printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt));
-printfunc("mac_y_dscale", diff(y_expand, scale));
-printfunc("mac_y_dk0", diff(y_expand, k0));
-printfunc("mac_y_dk1", diff(y_expand, k1));
diff --git a/src/lsq_rectilinear.mac b/src/lsq_rectilinear.mac
new file mode 100644
index 0000000..662680a
--- /dev/null
+++ b/src/lsq_rectilinear.mac
@@ -0,0 +1,40 @@
+/*
+ * This is the basic pinhole projection model with distortion correction
+ */
+
+x : tan(m_view - c_view)$
+y : tan(c_nick - m_nick)$ 
+x_rot : y * sin(c_tilt) + x * cos(c_tilt)$
+y_rot : y * cos(c_tilt) - x * sin(c_tilt)$
+d : x_rot ^ 2 + y_rot ^ 2$
+dist_fact : d ^ 2 * k1 + d * k0$
+x_dist : x_rot * (1 + dist_fact) * scale$
+y_dist : y_rot * (1 + dist_fact) * scale$
+
+/*
+ * Some mangling for C code generation
+ */
+
+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)"$
+
+printfunc(name, expression) := sprint("double ProjectionRectilinear::", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
+")$
+
+printfunc("mac_x", x_expand)$
+printfunc("mac_y", y_expand)$
+
+printfunc("mac_x_dc_view", diff(x_expand, c_view))$
+printfunc("mac_x_dc_nick", diff(x_expand, c_nick))$
+printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt))$
+printfunc("mac_x_dscale", diff(x_expand, scale))$
+printfunc("mac_x_dk0", diff(x_expand, k0))$
+printfunc("mac_x_dk1", diff(x_expand, k1))$
+printfunc("mac_y_dc_view", diff(y_expand, c_view))$
+printfunc("mac_y_dc_nick", diff(y_expand, c_nick))$
+printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt))$
+printfunc("mac_y_dscale", diff(y_expand, scale))$
+printfunc("mac_y_dk0", diff(y_expand, k0))$
+printfunc("mac_y_dk1", diff(y_expand, k1))$