From 23900d0b19e5c978934108f6a54aa82f6c7f27cb Mon Sep 17 00:00:00 2001 From: Johannes Hofmann Date: Sun, 17 Dec 2006 16:13:32 +0100 Subject: use LSQ for all projection types --- src/GipfelWidget.H | 4 +- src/GipfelWidget.cxx | 6 +- src/Makefile.am | 12 +- src/Makefile.lsq_funcs | 9 +- src/Panorama.H | 10 +- src/Panorama.cxx | 18 +-- src/Projection.H | 30 ----- src/Projection.cxx | 26 ----- src/ProjectionCylindrical.H | 35 ++++++ src/ProjectionCylindrical.cxx | 12 ++ src/ProjectionLSQ.H | 54 +++++++++ src/ProjectionLSQ.cxx | 249 ++++++++++++++++++++++++++++++++++++++++ src/ProjectionRectilinear.H | 35 ++++++ src/ProjectionRectilinear.cxx | 13 +++ src/ProjectionTangentialLSQ.H | 25 ---- src/ProjectionTangentialLSQ.cxx | 248 --------------------------------------- src/gipfel.cxx | 6 +- src/lsq_cylindrical.mac | 40 +++++++ src/lsq_funcs.mac | 40 ------- src/lsq_rectilinear.mac | 40 +++++++ 20 files changed, 512 insertions(+), 400 deletions(-) delete mode 100644 src/Projection.H delete mode 100644 src/Projection.cxx create mode 100644 src/ProjectionCylindrical.H create mode 100644 src/ProjectionCylindrical.cxx create mode 100644 src/ProjectionLSQ.H create mode 100644 src/ProjectionLSQ.cxx create mode 100644 src/ProjectionRectilinear.H create mode 100644 src/ProjectionRectilinear.cxx delete mode 100644 src/ProjectionTangentialLSQ.H delete mode 100644 src/ProjectionTangentialLSQ.cxx create mode 100644 src/lsq_cylindrical.mac delete mode 100644 src/lsq_funcs.mac create mode 100644 src/lsq_rectilinear.mac 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 #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 -// -// 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 -// -// This software may be used and distributed according to the terms -// of the GNU General Public License, incorporated herein by reference. - -#include -#include - -#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 +// +// 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 +// +// This software may be used and distributed according to the terms +// of the GNU General Public License, incorporated herein by reference. + +#include +#include + +#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 +// +// 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/ProjectionLSQ.cxx b/src/ProjectionLSQ.cxx new file mode 100644 index 0000000..3b953c8 --- /dev/null +++ b/src/ProjectionLSQ.cxx @@ -0,0 +1,249 @@ +// +// Copyright 2006 Johannes Hofmann +// +// This software may be used and distributed according to the terms +// of the GNU General Public License, incorporated herein by reference. + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "ProjectionLSQ.H" + + +double +ProjectionLSQ::sec(double a) { + return 1.0 / cos(a); +} + +int +ProjectionLSQ::comp_params(const Hills *h, ViewParams *parms) { + const Hill *tmp, *m1, *m2; + double scale_tmp; + int distortion_correct = 0; + + if (h->get_num() < 2) { + fprintf(stderr, "Please position at least 2 hills\n"); + return 1; + } else if (h->get_num() > 3) { + fprintf(stderr, "Performing calibration\n"); + distortion_correct = 1; + parms->k0 = 0.0; + parms->k1 = 0.0; + } + + m1 = h->get(0); + m2 = h->get(1); + + scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x); + if (isnan(scale_tmp) || scale_tmp < 100.0) { + // try again with mountains swapped + tmp = m1; + m1 = m2; + m2 = tmp; + scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x); + } + + + if (isnan(scale_tmp) || scale_tmp < 100.0) { + return 1; + } else { + + parms->a_center = m1->alph; + parms->scale = scale_tmp; + parms->a_nick = 0.0; + parms->a_tilt = 0.0; + + lsq(h, parms, 0); + + if (distortion_correct) { + lsq(h, parms, 1); + } + + return 0; + } +} + +struct data { + ProjectionLSQ *p; + int distortion_correct; + const Hills *h; + const ViewParams *old_params; +}; + +#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) { + struct data *dat = (struct data *) data; + double c_view, c_nick, c_tilt, scale, k0, k1; + + c_view = gsl_vector_get (x, 0); + c_nick = gsl_vector_get (x, 1); + c_tilt = gsl_vector_get (x, 2); + scale = gsl_vector_get (x, 3); + if (dat->distortion_correct) { + k0 = gsl_vector_get (x, 4); + k1 = gsl_vector_get (x, 5); + } else { + k0 = dat->old_params->k0; + k1 = dat->old_params->k1; + } + + for (int i=0; ih->get_num(); i++) { + Hill *m = dat->h->get(i); + + double mx = CALL(mac_x); + double my = CALL(mac_y); + + gsl_vector_set (f, i*2, mx - m->x); + gsl_vector_set (f, i*2+1, my - m->y); + } + + return GSL_SUCCESS; +} + + +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; + + c_view = gsl_vector_get (x, 0); + c_nick = gsl_vector_get (x, 1); + c_tilt = gsl_vector_get (x, 2); + scale = gsl_vector_get (x, 3); + if (dat->distortion_correct) { + k0 = gsl_vector_get (x, 4); + k1 = gsl_vector_get (x, 5); + } else { + k0 = dat->old_params->k0; + k1 = dat->old_params->k1; + } + + for (int i=0; ih->get_num(); i++) { + Hill *m = dat->h->get(i); + + gsl_matrix_set (J, 2*i, 0, CALL(mac_x_dc_view)); + gsl_matrix_set (J, 2*i, 1, CALL(mac_x_dc_nick)); + gsl_matrix_set (J, 2*i, 2, CALL(mac_x_dc_tilt)); + gsl_matrix_set (J, 2*i, 3, CALL(mac_x_dscale)); + 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+1, 0, CALL(mac_y_dc_view)); + gsl_matrix_set (J, 2*i+1, 1, CALL(mac_y_dc_nick)); + gsl_matrix_set (J, 2*i+1, 2, CALL(mac_y_dc_tilt)); + gsl_matrix_set (J, 2*i+1, 3, CALL(mac_y_dscale)); + 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)); + } + } + + return GSL_SUCCESS; +} + +static int +lsq_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) { + lsq_f (x, data, f); + lsq_df (x, data, J); + + return GSL_SUCCESS; +} + +int +ProjectionLSQ::lsq(const Hills *h, ViewParams *parms, + int distortion_correct) { + + const gsl_multifit_fdfsolver_type *T; + gsl_multifit_fdfsolver *s; + gsl_multifit_function_fdf f; + struct data dat; + double x_init[8]; + gsl_vector_view x; + int status; + int num_params = distortion_correct?6:4; + + dat.p = this; + dat.distortion_correct = distortion_correct; + dat.h = h; + dat.old_params = parms; + + x_init[0] = parms->a_center; + x_init[1] = parms->a_nick; + x_init[2] = parms->a_tilt; + x_init[3] = parms->scale; + x_init[4] = parms->k0; + x_init[5] = parms->k1; + + x = gsl_vector_view_array (x_init, num_params); + + f.f = &lsq_f; + f.df = &lsq_df; + f.fdf = &lsq_fdf; + f.n = h->get_num() * 2; + f.p = num_params; + f.params = &dat; + + T = gsl_multifit_fdfsolver_lmsder; + s = gsl_multifit_fdfsolver_alloc (T, h->get_num() * 2, num_params); + gsl_multifit_fdfsolver_set (s, &f, &x.vector); + + for (int i=0; i<100; i++) { + + status = gsl_multifit_fdfsolver_iterate (s); + if (status) { + fprintf(stderr, "gsl_multifit_fdfsolver_iterate: %d\n", status); + break; + } + } + + parms->a_center = gsl_vector_get(s->x, 0); + parms->a_nick = gsl_vector_get(s->x, 1); + parms->a_tilt = gsl_vector_get(s->x, 2); + parms->scale = gsl_vector_get(s->x, 3); + + if (distortion_correct) { + parms->k0 = gsl_vector_get(s->x, 4); + parms->k1 = gsl_vector_get(s->x, 5); + } + + gsl_multifit_fdfsolver_free (s); + + fprintf(stderr, "k0 %f, k1 %f\n", parms->k0, parms->k1); + + return 0; +} + +void +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); + *y = mac_y(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale, + parms->k0, parms->k1, alph, a_nick); +} + +double +ProjectionLSQ::comp_scale(double a1, double a2, double d1, double d2) { + double sign1 = 1.0; + double sc, tan_a1, tan_a2; + + tan_a1 = tan(a1); + tan_a2 = tan(a2); + + sc = ((((1.0 + (tan_a1 * tan_a2)) * (d1 - d2)) - (sign1 * pow((((1.0 + pow((tan_a1 * tan_a2), 2.0)) * ((d1 * d1) + (d2 * d2))) + (2.0 * ((tan_a1 * tan_a2 * pow((d1 + d2), 2.0)) - (d1 * d2 * (((tan_a1 * tan_a1) * (2.0 + (tan_a2 * tan_a2))) + 1.0 + (2.0 * (tan_a2 * tan_a2))))))), (1.0 / 2.0)))) / (2.0 * (tan_a1 - tan_a2))); + + return sc; +} 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 +// +// 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 +// +// This software may be used and distributed according to the terms +// of the GNU General Public License, incorporated herein by reference. + +#include +#include + +#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 -// -// 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/ProjectionTangentialLSQ.cxx b/src/ProjectionTangentialLSQ.cxx deleted file mode 100644 index 76ec2fa..0000000 --- a/src/ProjectionTangentialLSQ.cxx +++ /dev/null @@ -1,248 +0,0 @@ -// -// Copyright 2006 Johannes Hofmann -// -// This software may be used and distributed according to the terms -// of the GNU General Public License, incorporated herein by reference. - -#include -#include -#include -#include - -#include -#include -#include -#include -#include - -#include "ProjectionTangentialLSQ.H" - - -static double sec(double a) { - return 1.0 / cos(a); -} - -#include "lsq_funcs.c" - -int -ProjectionTangentialLSQ::comp_params(const Hills *h, ViewParams *parms) { - const Hill *tmp, *m1, *m2; - double scale_tmp; - int distortion_correct = 0; - - if (h->get_num() < 2) { - fprintf(stderr, "Please position at least 2 hills\n"); - return 1; - } else if (h->get_num() > 3) { - fprintf(stderr, "Performing calibration\n"); - distortion_correct = 1; - parms->k0 = 0.0; - parms->k1 = 0.0; - } - - m1 = h->get(0); - m2 = h->get(1); - - scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x); - if (isnan(scale_tmp) || scale_tmp < 100.0) { - // try again with mountains swapped - tmp = m1; - m1 = m2; - m2 = tmp; - scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x); - } - - - if (isnan(scale_tmp) || scale_tmp < 100.0) { - return 1; - } else { - - parms->a_center = m1->alph; - parms->scale = scale_tmp; - parms->a_nick = 0.0; - parms->a_tilt = 0.0; - - lsq(h, parms, 0); - - if (distortion_correct) { - lsq(h, parms, 1); - } - - return 0; - } -} - -struct data { - 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) - -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; - - c_view = gsl_vector_get (x, 0); - c_nick = gsl_vector_get (x, 1); - c_tilt = gsl_vector_get (x, 2); - scale = gsl_vector_get (x, 3); - if (dat->distortion_correct) { - k0 = gsl_vector_get (x, 4); - k1 = gsl_vector_get (x, 5); - } else { - k0 = dat->old_params->k0; - k1 = dat->old_params->k1; - } - - for (int i=0; ih->get_num(); i++) { - Hill *m = dat->h->get(i); - - double mx = CALL(mac_x); - double my = CALL(mac_y); - - gsl_vector_set (f, i*2, mx - m->x); - gsl_vector_set (f, i*2+1, my - m->y); - } - - return GSL_SUCCESS; -} - - -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; - - c_view = gsl_vector_get (x, 0); - c_nick = gsl_vector_get (x, 1); - c_tilt = gsl_vector_get (x, 2); - scale = gsl_vector_get (x, 3); - if (dat->distortion_correct) { - k0 = gsl_vector_get (x, 4); - k1 = gsl_vector_get (x, 5); - } else { - k0 = dat->old_params->k0; - k1 = dat->old_params->k1; - } - - for (int i=0; ih->get_num(); i++) { - Hill *m = dat->h->get(i); - - gsl_matrix_set (J, 2*i, 0, CALL(mac_x_dc_view)); - gsl_matrix_set (J, 2*i, 1, CALL(mac_x_dc_nick)); - gsl_matrix_set (J, 2*i, 2, CALL(mac_x_dc_tilt)); - gsl_matrix_set (J, 2*i, 3, CALL(mac_x_dscale)); - 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+1, 0, CALL(mac_y_dc_view)); - gsl_matrix_set (J, 2*i+1, 1, CALL(mac_y_dc_nick)); - gsl_matrix_set (J, 2*i+1, 2, CALL(mac_y_dc_tilt)); - gsl_matrix_set (J, 2*i+1, 3, CALL(mac_y_dscale)); - 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)); - } - } - - return GSL_SUCCESS; -} - -static int -lsq_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) { - lsq_f (x, data, f); - lsq_df (x, data, J); - - return GSL_SUCCESS; -} - -int -ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms, - int distortion_correct) { - - const gsl_multifit_fdfsolver_type *T; - gsl_multifit_fdfsolver *s; - gsl_multifit_function_fdf f; - struct data dat; - double x_init[8]; - gsl_vector_view x; - int status; - int num_params = distortion_correct?6:4; - - dat.distortion_correct = distortion_correct; - dat.h = h; - dat.old_params = parms; - - x_init[0] = parms->a_center; - x_init[1] = parms->a_nick; - x_init[2] = parms->a_tilt; - x_init[3] = parms->scale; - x_init[4] = parms->k0; - x_init[5] = parms->k1; - - x = gsl_vector_view_array (x_init, num_params); - - f.f = &lsq_f; - f.df = &lsq_df; - f.fdf = &lsq_fdf; - f.n = h->get_num() * 2; - f.p = num_params; - f.params = &dat; - - T = gsl_multifit_fdfsolver_lmsder; - s = gsl_multifit_fdfsolver_alloc (T, h->get_num() * 2, num_params); - gsl_multifit_fdfsolver_set (s, &f, &x.vector); - - for (int i=0; i<100; i++) { - - status = gsl_multifit_fdfsolver_iterate (s); - if (status) { - fprintf(stderr, "gsl_multifit_fdfsolver_iterate: %d\n", status); - break; - } - } - - parms->a_center = gsl_vector_get(s->x, 0); - parms->a_nick = gsl_vector_get(s->x, 1); - parms->a_tilt = gsl_vector_get(s->x, 2); - parms->scale = gsl_vector_get(s->x, 3); - - if (distortion_correct) { - parms->k0 = gsl_vector_get(s->x, 4); - parms->k1 = gsl_vector_get(s->x, 5); - } - - gsl_multifit_fdfsolver_free (s); - - fprintf(stderr, "k0 %f, k1 %f\n", parms->k0, parms->k1); - - return 0; -} - -void -ProjectionTangentialLSQ::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); - *y = mac_y(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale, - parms->k0, parms->k1, alph, a_nick); -} - -double -ProjectionTangentialLSQ::comp_scale(double a1, double a2, double d1, double d2) { - double sign1 = 1.0; - double sc, tan_a1, tan_a2; - - tan_a1 = tan(a1); - tan_a2 = tan(a2); - - sc = ((((1.0 + (tan_a1 * tan_a2)) * (d1 - d2)) - (sign1 * pow((((1.0 + pow((tan_a1 * tan_a2), 2.0)) * ((d1 * d1) + (d2 * d2))) + (2.0 * ((tan_a1 * tan_a2 * pow((d1 + d2), 2.0)) - (d1 * d2 * (((tan_a1 * tan_a1) * (2.0 + (tan_a2 * tan_a2))) + 1.0 + (2.0 * (tan_a2 * tan_a2))))))), (1.0 / 2.0)))) / (2.0 * (tan_a1 - tan_a2))); - - return sc; -} 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))$ -- cgit v1.2.3