diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/ProjectionTangentialLSQ.H | 29 | ||||
| -rw-r--r-- | src/ProjectionTangentialLSQ.cxx | 225 | ||||
| -rw-r--r-- | src/lsq_funcs.mac | 26 |
4 files changed, 269 insertions, 13 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index 8ba69c1..1cae9b7 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -7,6 +7,7 @@ gipfel_SOURCES = \ Panorama.cxx \ Projection.cxx \ ProjectionTangential.cxx \ + ProjectionTangentialLSQ.cxx \ ProjectionSphaeric.cxx \ Hill.cxx \ Fl_Value_Dial.cxx \ @@ -24,6 +25,7 @@ noinst_HEADERS = \ Panorama.H \ Projection.H \ ProjectionTangential.H \ + ProjectionTangentialLSQ.H \ ProjectionSphaeric.H \ Hill.H \ ViewParams.H \ diff --git a/src/ProjectionTangentialLSQ.H b/src/ProjectionTangentialLSQ.H new file mode 100644 index 0000000..492860d --- /dev/null +++ b/src/ProjectionTangentialLSQ.H @@ -0,0 +1,29 @@ +// +// 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_center_angle(double alph_a, double alph_b, double d1, double d2); + + double comp_scale(double alph_a, double alph_b, double d1, double d2); + + double angle_dist(double a1, double a2); + + int optimize(const Hill *m1, const Hill *m2, ViewParams *parms); + + public: + void get_coordinates(double a_view, double a_nick, const ViewParams *parms, + double *x, double *y); + + int comp_params(const Hill *m1, const Hill *m2, ViewParams *parms); +}; +#endif diff --git a/src/ProjectionTangentialLSQ.cxx b/src/ProjectionTangentialLSQ.cxx new file mode 100644 index 0000000..6ae7757 --- /dev/null +++ b/src/ProjectionTangentialLSQ.cxx @@ -0,0 +1,225 @@ +// +// 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 <stdlib.h> +#include <string.h> +#include <math.h> + +#include "ProjectionTangentialLSQ.H" + + +static double sec(double a) { + return 1.0 / cos(a); +} + +#include "lsq_funcs.c" + +static double +comp_tilt(double tan_nick_view, double tan_dir_view, double n_scale, + double tan_nick_m, double tan_dir_m, + double x, double y, double pi_d); + +int +ProjectionTangentialLSQ::comp_params(const Hill *m1, const Hill *m2, ViewParams *parms) { + const Hill *tmp; + double a_center_tmp, scale_tmp, a_nick_tmp; + + 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); + } + + 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)) { + return 1; + } else { + + parms->a_center = a_center_tmp; + parms->scale = scale_tmp; + parms->a_nick = a_nick_tmp; + + optimize(m1, m2, parms); + + return 0; + } +} + +double +ProjectionTangentialLSQ::angle_dist(double a1, double a2) { + double ret; + + a1 = fmod(a1, 2.0 * pi_d); + if (a1 < 0.0) { + a1 = a1 + 2.0 * pi_d; + } + a2 = fmod(a2, 2.0 * pi_d); + if (a2 < 0.0) { + a2 = a2 + 2.0 * pi_d; + } + + ret = fabs(a1 - a2); + if (ret > pi_d) { + ret = 2.0 * pi_d - ret; + } + + return ret; +} + +int +ProjectionTangentialLSQ::optimize(const Hill *m1, const Hill *m2, ViewParams *parms) { + int i; + double tan_nick_view, tan_dir_view, n_scale; + double tan_nick_m1, tan_dir_m1; + double tan_nick_m2, tan_dir_m2; + double d_m1_2, d_m2_2, d_m1_m2_2; + + d_m1_2 = pow(m1->x, 2.0) + pow(m1->y, 2.0); + d_m2_2 = pow(m2->x, 2.0) + pow(m2->y, 2.0); + d_m1_m2_2 = pow(m1->x - m2->x, 2.0) + pow(m1->y - m2->y, 2.0); + + tan_nick_view = tan(parms->a_nick); + tan_dir_view = tan(parms->a_center); + n_scale = parms->scale; + tan_dir_m1 = tan(m1->alph); + tan_nick_m1 = tan(m1->a_nick); + tan_dir_m2 = tan(m2->alph); + tan_nick_m2 = tan(m2->a_nick); +#if 0 + for (i=0; i<5; i++) { + opt_step(&tan_nick_view, &tan_dir_view, &n_scale, + tan_dir_m1, tan_nick_m1, tan_dir_m2, tan_nick_m2, + d_m1_2, d_m2_2, d_m1_m2_2); + } +#endif + + if (isnan(tan_dir_view) || isnan(tan_nick_view) || isnan(n_scale)) { + fprintf(stderr, "No solution found.\n"); + return 1; + } + + parms->a_center = atan(tan_dir_view); + parms->a_nick = atan(tan_nick_view); + + if (parms->a_center > 2.0 * pi_d) { + parms->a_center = parms->a_center - 2.0 * pi_d; + } else if (parms->a_center < 0.0) { + parms->a_center = parms->a_center + 2.0 * pi_d; + } + + // atan(tan_dir_view) is not the only possible solution. + // Choose the one which is close to m1->alph. + if (angle_dist(parms->a_center, m1->alph) > pi_d/2.0) { + parms->a_center = parms->a_center + pi_d; + } + + parms->scale = n_scale; + + // use the point with greater distance from center for tilt computation + if (d_m1_2 > d_m2_2) { + parms->a_tilt = comp_tilt(tan_nick_view, tan_dir_view, n_scale, + tan_nick_m1, tan_dir_m1, + (double) m1->x, (double) m1->y, pi_d); + } else { + parms->a_tilt = comp_tilt(tan_nick_view, tan_dir_view, n_scale, + tan_nick_m2, tan_dir_m2, + (double) m2->x, (double) m2->y, pi_d); + } + + return 0; +} + +void +ProjectionTangentialLSQ::get_coordinates(double a_view, double a_nick, + const ViewParams *parms, double *x, double *y) { + double x_tmp, y_tmp; + + x_tmp = tan(a_view) * parms->scale; + y_tmp = - (tan(a_nick - parms->a_nick) * parms->scale); + + // rotate by a_tilt; + *x = x_tmp * cos(parms->a_tilt) - y_tmp * sin(parms->a_tilt); + *y = x_tmp * sin(parms->a_tilt) + y_tmp * cos(parms->a_tilt); +} + +double +ProjectionTangentialLSQ::comp_center_angle(double a1, double a2, double d1, double d2) { + double sign1 = 1.0; + double tan_acenter, tan_a1, tan_a2, a_center; + + tan_a1 = tan(a1); + tan_a2 = tan(a2); + + tan_acenter = (((pow(((pow((1.0 + (tan_a1 * tan_a2)), 2.0) * ((d1 * d1) + (d2 * d2))) + (2.0 * d1 * d2 * ((2.0 * ((tan_a2 * tan_a1) - (tan_a2 * tan_a2))) - ((tan_a1 * tan_a1) * (2.0 + (tan_a2 * tan_a2))) - 1.0))), (1.0 / 2.0)) * sign1) + ((1.0 - (tan_a1 * tan_a2)) * (d1 - d2))) / (2.0 * ((d2 * tan_a2) - (d1 * tan_a1)))); + + a_center = atan(tan_acenter); + + if (a_center > 2.0 * pi_d) { + a_center = a_center - 2.0 * pi_d; + } else if (a_center < 0.0) { + a_center = a_center + 2.0 * pi_d; + } + + // atan(tan_dir_view) is not the only possible solution. + // Choose the one which is close to m1->alph. + if (fabs(a_center - a1) > pi_d/2.0) { + a_center = a_center + pi_d; + } + + return a_center; +} + +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; +} + +static double +comp_tilt(double tan_nick_view, double tan_dir_view, double n_scale, + double tan_nick_m, double tan_dir_m, + double x, double y, double pi_d) { + double y_tmp, x_tmp, sin_a_tilt1, sin_a_tilt2, sin_a_tilt, res; + + y_tmp = - (((tan_nick_view - tan_nick_m) * n_scale) / + (tan_nick_m * tan_nick_view + 1)); + x_tmp = - (((tan_dir_view - tan_dir_m) * n_scale) / + (tan_dir_m * tan_dir_view + 1)); + + + sin_a_tilt1 = - (y * - pow(x*x + y*y - y_tmp*y_tmp, 0.5) - x * y_tmp) / + (x*x + y*y); + + sin_a_tilt2 = - (y * pow(x*x + y*y - y_tmp*y_tmp, 0.5) - x * y_tmp) / + (x*x + y*y); + + sin_a_tilt = fabs(sin_a_tilt1) < fabs(sin_a_tilt2)?sin_a_tilt1:sin_a_tilt2; + + res = asin(sin_a_tilt); + + if (res > pi_d / 4.0) { + res = res - pi_d / 2.0; + } else if (res < -pi_d / 4.0) { + res = res + pi_d / 2.0; + } + + return res; +} diff --git a/src/lsq_funcs.mac b/src/lsq_funcs.mac index f2a70b9..3e1bc18 100644 --- a/src/lsq_funcs.mac +++ b/src/lsq_funcs.mac @@ -15,29 +15,29 @@ y : y_unrot * cos(c_tilt) - x_unrot * sin(c_tilt)$ x_expand : trigexpand(x)$ y_expand : trigexpand(y)$ -args: "(double c_view, double c_nick, 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 m_view, double m_nick)"$ -sprint("double mac_x", args, "{ return",string(subst(pow, "^", x_expand)),"}"," +sprint("double mac_x", args, "{ return",string(subst(pow, "^", x_expand)),";}"," ")$ -sprint("double mac_y", args, "{ return", string(subst(pow, "^",y_expand)),"}"," +sprint("double mac_y", args, "{ return", string(subst(pow, "^",y_expand)),";}"," ")$ -sprint("double mac_x_dc_view", args, "{ return", string(subst(pow, "^",diff(x_expand, c_view))),"}"," +sprint("double mac_x_dc_view", args, "{ return", string(subst(pow, "^",diff(x_expand, c_view))),";}"," ")$ -sprint("double mac_x_dc_nick", args, "{ return", string(subst(pow, "^",diff(x_expand, c_nick))),"}"," +sprint("double mac_x_dc_nick", args, "{ return", string(subst(pow, "^",diff(x_expand, c_nick))),";}"," ")$ -sprint("double mac_x_dscale", args, "{ return",string(subst(pow, "^",diff(x_expand, scale))),"}"," +sprint("double mac_x_dscale", args, "{ return",string(subst(pow, "^",diff(x_expand, scale))),";}"," ")$ -sprint("double mac_x_dk0", args, "{ return",string(subst(pow, "^",diff(x_expand, k0))),"}"," +sprint("double mac_x_dk0", args, "{ return",string(subst(pow, "^",diff(x_expand, k0))),";}"," ")$ -sprint("double mac_x_dk1",args, "{ return", string(subst(pow, "^",diff(x_expand, k1))),"}"," +sprint("double mac_x_dk1",args, "{ return", string(subst(pow, "^",diff(x_expand, k1))),";}"," ")$ -sprint("double mac_y_dc_view",args, "{ return", string(subst(pow, "^",diff(y_expand, c_view))),"}"," +sprint("double mac_y_dc_view",args, "{ return", string(subst(pow, "^",diff(y_expand, c_view))),";}"," ")$ -sprint("double mac_y_dc_nick",args, "{ return", string(subst(pow, "^",diff(y_expand, c_nick))),"}"," +sprint("double mac_y_dc_nick",args, "{ return", string(subst(pow, "^",diff(y_expand, c_nick))),";}"," ")$ -sprint("double mac_y_dscale", args, "{ return",string(subst(pow, "^",diff(y_expand, scale))),"}"," +sprint("double mac_y_dscale", args, "{ return",string(subst(pow, "^",diff(y_expand, scale))),";}"," ")$ -sprint("double mac_y_dk0", args, "{ return",string(subst(pow, "^",diff(y_expand, k0))),"}"," +sprint("double mac_y_dk0", args, "{ return",string(subst(pow, "^",diff(y_expand, k0))),";}"," ")$ -sprint("double mac_y_dk1",args, "{ return", string(subst(pow, "^",diff(y_expand, k1))),"}"," +sprint("double mac_y_dk1",args, "{ return", string(subst(pow, "^",diff(y_expand, k1))),";}"," ")$ |