diff options
-rw-r--r-- | src/GipfelWidget.cxx | 38 | ||||
-rw-r--r-- | src/Hill.H | 2 | ||||
-rw-r--r-- | src/Panorama.cxx | 9 | ||||
-rw-r--r-- | src/ProjectionTangentialLSQ.cxx | 95 | ||||
-rw-r--r-- | src/lsq_funcs.mac | 24 |
5 files changed, 67 insertions, 101 deletions
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx index abbc393..f63cc67 100644 --- a/src/GipfelWidget.cxx +++ b/src/GipfelWidget.cxx @@ -214,8 +214,6 @@ void GipfelWidget::draw() { Hills *mnts; Hill *m; - int center_x = w() / 2; - int center_y = h() / 2; int i; if (img == NULL) { @@ -230,6 +228,8 @@ GipfelWidget::draw() { mnts = pan->get_visible_mountains(); for (i=0; i<mnts->get_num(); i++) { m = mnts->get(i); + int m_x = w() / 2 + x() + (int) rint(m->x); + int m_y = h() / 2 + y() + (int) rint(m->y); if (m->flags & (Hill::DUPLICATE|Hill::TRACK_POINT)) { continue; @@ -241,29 +241,17 @@ GipfelWidget::draw() { if (known_hills->contains(m)) { fl_color(FL_RED); - draw_flag(center_x + m->x + x(), center_y + m->y + y(), "1"); + draw_flag(m_x, m_y, "1"); } else if (m->flags & Hill::HIDDEN) { fl_color(FL_BLUE); } else { fl_color(FL_BLACK); } - fl_xyline(center_x + m->x + x() - CROSS_SIZE, center_y + m->y + y(), center_x + m->x + x() + CROSS_SIZE); - fl_yxline(center_x + m->x + x(), center_y + m->y + m->label_y + y() - CROSS_SIZE, center_y + m->y + y() + CROSS_SIZE); - - fl_draw(m->name, - center_x + m->x + x(), - center_y + m->y + m->label_y + y()); - } - - /* markers */ - for (i=0; i<marker->get_num(); i++) { - m = marker->get(i); + fl_xyline(m_x- CROSS_SIZE, m_y, m_x + CROSS_SIZE); + fl_yxline(m_x, m_y + m->label_y - CROSS_SIZE, m_y + CROSS_SIZE); - fl_color(FL_GREEN); - fl_xyline(center_x + m->x + x() - CROSS_SIZE * 2, center_y + m->y + y(), center_x + m->x + x() + CROSS_SIZE * 2); - fl_yxline(center_x + m->x + x(), center_y + m->y + y() - CROSS_SIZE * 2, center_y + m->y + y() + CROSS_SIZE * 2); - draw_flag(center_x + m->x + x(), center_y + m->y + y(), NULL); + fl_draw(m->name, m_x , m_y + m->label_y); } /* track */ @@ -271,6 +259,9 @@ GipfelWidget::draw() { int last_x, last_y, last_initialized = 0; for (i=1; i<track_points->get_num(); i++) { + m = mnts->get(i); + int m_x = w() / 2 + x() + (int) rint(m->x); + int m_y = h() / 2 + y() + (int) rint(m->y); if (!(track_points->get(i)->flags & Hill::VISIBLE)) { continue; } @@ -285,14 +276,13 @@ GipfelWidget::draw() { get_rel_track_width(track_points->get(i))); if (last_initialized) { fl_begin_line(); - fl_vertex(center_x + x() + last_x, center_y + y() + last_y); - fl_vertex(center_x + x() + track_points->get(i)->x, - center_y + y() + track_points->get(i)->y); + fl_vertex(last_x, last_y); + fl_vertex(m_x, m_y); fl_end_line(); } - last_x = track_points->get(i)->x; - last_y = track_points->get(i)->y; + last_x = m_x; + last_y = m_y; last_initialized++; } fl_line_style(0); @@ -303,7 +293,7 @@ GipfelWidget::draw() { static int -overlap(int m1, int n1, int m2, int n2) { +overlap(double m1, double n1, double m2, double n2) { return m1 <= n2 && n1 >= m2; } @@ -27,7 +27,7 @@ class Hill { double a_view; double a_nick; double dist; - int x, y; + double x, y; int label_x, label_y; char *name; int flags; diff --git a/src/Panorama.cxx b/src/Panorama.cxx index b227ebf..a62e2f0 100644 --- a/src/Panorama.cxx +++ b/src/Panorama.cxx @@ -164,7 +164,7 @@ Panorama::guess(Hills *p, Hill *m1) { Hills h; double best = 100000000.0, v; double a_center_best, a_nick_best, a_tilt_best, scale_best; - int x1_sav, y1_sav; + double x1_sav, y1_sav; int i, j; if (m1 == NULL) { @@ -475,7 +475,6 @@ Panorama::update_visible_mountains() { m = close_mountains->get(i); m->a_view = m->alph - parms.a_center; - if (m->a_view > pi_d) { m->a_view -= 2.0*pi_d; } else if (m->a_view < -pi_d) { @@ -499,11 +498,7 @@ Panorama::update_coordinates() { for (int i=0; i<visible_mountains->get_num(); i++) { m = visible_mountains->get(i); - double tmp_x, tmp_y; - - proj->get_coordinates(m->a_view, m->a_nick, &parms, &tmp_x, &tmp_y); - m->x = (int) rint(tmp_x); - m->y = (int) rint(tmp_y); + proj->get_coordinates(m->a_view, m->a_nick, &parms, &m->x, &m->y); } } diff --git a/src/ProjectionTangentialLSQ.cxx b/src/ProjectionTangentialLSQ.cxx index b201d31..059ad5d 100644 --- a/src/ProjectionTangentialLSQ.cxx +++ b/src/ProjectionTangentialLSQ.cxx @@ -22,6 +22,8 @@ static double sec(double a) { return 1.0 / cos(a); } +static double pi_d = asin(1.0) * 2.0, deg2rad = pi_d / 180.0; + #include "lsq_funcs.c" static double @@ -57,24 +59,16 @@ ProjectionTangentialLSQ::comp_params(const Hills *h, ViewParams *parms) { scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x); } - a_center_tmp = comp_center_angle(m1->alph, m2->alph, m1->x, m2->x); - a_nick_tmp = atan ((m1->y + tan(m1->a_nick) * parms->scale) / - (parms->scale - m1->y * tan(m1->a_nick))); - if (isnan(a_center_tmp) || isnan(scale_tmp) || - scale_tmp < 100.0 || isnan(a_nick_tmp)) { + if (isnan(scale_tmp) || scale_tmp < 100.0) { return 1; } else { - parms->a_center = a_center_tmp; + parms->a_center = m1->alph; parms->scale = scale_tmp; - parms->a_nick = a_nick_tmp; + parms->a_nick = 0.0; parms->a_tilt = 0.0; - if (angle_dist(parms->a_center, m1->alph) > pi_d/2.0) { - parms->a_center = parms->a_center + pi_d; - } - lsq(h, parms); return 0; @@ -107,7 +101,7 @@ struct data { const ViewParams *old_params; }; -#define CALL(A) A(c_view, c_nick, c_tilt, scale, k0, k1, u0, v0, m->a_view, m->a_nick) +#define CALL(A) A(c_view, c_nick, c_tilt, scale, k0, k1, m->a_view, m->a_nick) static int lsq_f (const gsl_vector * x, void *data, gsl_vector * f) { @@ -121,30 +115,19 @@ lsq_f (const gsl_vector * x, void *data, gsl_vector * f) { if (x->size >= 6) { k0 = gsl_vector_get (x, 4); k1 = gsl_vector_get (x, 5); -#if 0 - u0 = gsl_vector_get (x, 6); - v0 = gsl_vector_get (x, 7); -#else - u0 = dat->old_params->u0; - v0 = dat->old_params->v0; -#endif - - } else { k0 = dat->old_params->k0; k1 = dat->old_params->k1; - u0 = dat->old_params->u0; - v0 = dat->old_params->v0; } for (int i=0; i<dat->h->get_num(); i++) { Hill *m = dat->h->get(i); - double x = CALL(mac_x); - double y = CALL(mac_y); + double mx = CALL(mac_x); + double my = CALL(mac_y); - gsl_vector_set (f, i*2, x - m->x); - gsl_vector_set (f, i*2+1, y - m->y); + gsl_vector_set (f, i*2, mx - m->x); + gsl_vector_set (f, i*2+1, my - m->y); } return GSL_SUCCESS; @@ -163,19 +146,9 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) { if (x->size >= 6) { k0 = gsl_vector_get (x, 4); k1 = gsl_vector_get (x, 5); -#if 0 - u0 = gsl_vector_get (x, 6); - v0 = gsl_vector_get (x, 7); -#else - u0 = dat->old_params->u0; - v0 = dat->old_params->v0; - -#endif } else { k0 = dat->old_params->k0; k1 = dat->old_params->k1; - u0 = dat->old_params->u0; - v0 = dat->old_params->v0; } for (int i=0; i<dat->h->get_num(); i++) { @@ -188,10 +161,6 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) { if (x->size >= 6) { gsl_matrix_set (J, 2*i, 4, CALL(mac_x_dk0)); gsl_matrix_set (J, 2*i, 5, CALL(mac_x_dk1)); -#if 0 - gsl_matrix_set (J, 2*i, 6, CALL(mac_x_du0)); - gsl_matrix_set (J, 2*i, 7, CALL(mac_x_dv0)); -#endif } gsl_matrix_set (J, 2*i+1, 0, CALL(mac_y_dc_view)); @@ -201,10 +170,6 @@ lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) { if (x->size >= 6) { gsl_matrix_set (J, 2*i+1, 4, CALL(mac_y_dk0)); gsl_matrix_set (J, 2*i+1, 5, CALL(mac_y_dk1)); -#if 0 - gsl_matrix_set (J, 2*i+1, 6, CALL(mac_y_du0)); - gsl_matrix_set (J, 2*i+1, 7, CALL(mac_y_dv0)); -#endif } } @@ -230,6 +195,11 @@ ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms) { int status; int num_params = h->get_num()>3?6:4; + fprintf(stderr, "x %f, y %f\n", + h->get(0)->x, + h->get(0)->y); + + dat.h = h; dat.old_params = parms; @@ -239,8 +209,6 @@ ProjectionTangentialLSQ::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->u0; - x_init[7] = parms->v0; x = gsl_vector_view_array (x_init, num_params); @@ -263,6 +231,7 @@ ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms) { break; } + fprintf(stderr, "%d, |f(x)| = %g\n", i, gsl_blas_dnrm2 (s->f)); } parms->a_center = gsl_vector_get(s->x, 0); @@ -273,16 +242,30 @@ ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms) { if (num_params == 6) { parms->k0 = gsl_vector_get(s->x, 4); parms->k1 = gsl_vector_get(s->x, 5); -#if 0 - parms->u0 = gsl_vector_get(s->x, 6); - parms->v0 = gsl_vector_get(s->x, 7); -#endif } gsl_multifit_fdfsolver_free (s); - fprintf(stderr, "k0 %f k1 %f u0 %f v0 %f\n", - parms->k0, parms->k1, parms->u0, parms->v0); + fprintf(stderr, "center %f, x %f, dx %f, y %f, dy %f\n", + parms->a_center / deg2rad, + h->get(0)->x, + h->get(0)->x - mac_x(parms->a_center, + parms->a_nick, + parms->a_tilt, + parms->scale, + parms->k0, + parms->k1, + h->get(0)->a_view, h->get(0)->a_nick), + h->get(0)->y, + h->get(0)->y - mac_y(parms->a_center, + parms->a_nick, + parms->a_tilt, + parms->scale, + parms->k0, + parms->k1, + h->get(0)->a_view, h->get(0)->a_nick)); + + return 0; } @@ -292,9 +275,11 @@ ProjectionTangentialLSQ::get_coordinates(double a_view, 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, parms->u0, parms->v0, a_view, a_nick); + parms->k0, parms->k1, a_view, a_nick); *y = mac_y(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale, - parms->k0, parms->k1, parms->u0, parms->v0, a_view, a_nick); + parms->k0, parms->k1, a_view, a_nick); + + fprintf(stderr, "==> %f %f\n", *x, *y); } double diff --git a/src/lsq_funcs.mac b/src/lsq_funcs.mac index 3a2cb33..de71198 100644 --- a/src/lsq_funcs.mac +++ b/src/lsq_funcs.mac @@ -2,23 +2,23 @@ * This is the basic pinhole projection model with distortion correction */ -x_undist_unrot : tan(m_view - c_view); -y_undist_unrot : tan(c_nick - m_nick); -d : (u0 - x_undist_unrot) ^ 2 + (v0 - y_undist_unrot) ^ 2; +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_unrot : x_undist_unrot * (1 + dist_fact); -y_unrot : y_undist_unrot * (1 + dist_fact); -x : scale * (y_unrot * sin(c_tilt) + x_unrot * cos(c_tilt)); -y : scale * (y_unrot * cos(c_tilt) - x_unrot * sin(c_tilt)); +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); -y_expand : trigexpand(y); +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 u0, double v0, double m_view, double m_nick)"; +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)), ";}", " "); @@ -32,13 +32,9 @@ 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_x_du0", diff(x_expand, u0)); -printfunc("mac_x_dv0", diff(x_expand, v0)); 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)); -printfunc("mac_y_du0", diff(x_expand, u0)); -printfunc("mac_y_dv0", diff(x_expand, v0)); |