1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
|
//
// Copyright 2005 by Johannes Hofmann
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Library General Public
// License as published by the Free Software Foundation; either
// version 2 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA.
//
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "ProjectionSphaeric.H"
#define BEST_UNDEF 10000000.0
int
ProjectionSphaeric::comp_params(const Hill *m1, const Hill *m2, ViewParams *parms) {
const Hill *m_tmp;
double tmp_x, tmp_y;
double val;
ViewParams best, tmp;
double best_val = BEST_UNDEF;
double d_m1_2, d_m2_2, d_m1_m2_2;
int i, j;
if (m1->x < m2->x) {
m_tmp = m1;
m1 = m2;
m2 = m_tmp;
}
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);
tmp.scale = comp_scale(m1, m2, d_m1_m2_2);
for(i=0; i<2; i++) { // we need to try four possible solutions ...
for(j=0; j<2; j++) {
tmp.a_center = comp_dir_view(m1, m2, d_m1_2, d_m2_2,
tmp.scale, i==0?1.0:-1.0);
tmp.a_nick = comp_nick_view(m1, m2, d_m1_2, tmp.scale, tmp.a_center,
j==0?1.0:-1.0);
// use the point with greater distance from center for tilt computation
if (d_m1_2 > d_m2_2) {
tmp.a_tilt = comp_tilt_view(m1, tmp.scale, tmp.a_center, tmp.a_nick);
} else {
tmp.a_tilt = comp_tilt_view(m2, tmp.scale, tmp.a_center, tmp.a_nick);
}
if (isnan(tmp.a_center) || isnan(tmp.scale) ||
isnan(tmp.a_nick) || isnan(tmp.a_tilt)) {
;
} else {
set_coordinates(m1, &tmp, &tmp_x, &tmp_y);
val = sqrt(pow(tmp_x - m1->x, 2.0) + pow(tmp_y - m1->y, 2.0));
set_coordinates(m2, &tmp, &tmp_x, &tmp_y);
val += sqrt(pow(tmp_x - m2->x, 2.0) + pow(tmp_y - m2->y, 2.0));
if (val < best_val) {
best_val = val;
best = tmp;
}
}
}
}
if (best_val < BEST_UNDEF) {
*parms = best;
return 0;
} else {
return 1;
}
}
void
ProjectionSphaeric::set_coordinates(Hill *m, const ViewParams *parms) {
double x_tmp, y_tmp;
set_coordinates(m, parms, &x_tmp, &y_tmp);
m->x = (int) rint(x_tmp);
m->y = (int) rint(y_tmp);
}
void
ProjectionSphaeric::set_coordinates(const Hill *m, const ViewParams *parms,
double *x, double *y) {
double x_tmp, y_tmp;
x_tmp = m->a_view * parms->scale;
y_tmp = - (m->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
ProjectionSphaeric::comp_scale(const Hill *m1, const Hill *m2, double d_m1_m2_2) {
double sign1 = 1.0;
double nick_m1 = m1->a_nick;
double nick_m2 = m2->a_nick;
double dir_m1 = m1->alph;
double dir_m2 = m2->alph;
return (pow((d_m1_m2_2 / ((dir_m2 * dir_m2) - (2.0 * ((dir_m2 * dir_m1) + (nick_m2 * nick_m1))) + (nick_m2 * nick_m2) + (dir_m1 * dir_m1) + (nick_m1 * nick_m1))), (1.0 / 2.0)) * sign1);
}
// using the sign3 parameter one can choose between the two possible solutions
// sign3 must be 1.0 or -1.0
double
ProjectionSphaeric::comp_dir_view(const Hill *m1, const Hill *m2, double d_m1_2, double d_m2_2, double scale, double sign3) {
double dir_view;
double nick_m1 = m1->a_nick;
double nick_m2 = m2->a_nick;
double dir_m1 = m1->alph;
double dir_m2 = m2->alph;
dir_view = (((pow((16.0 * pow(scale, 4.0) * (((scale * scale) * ((2.0 * (((scale * scale) * ((dir_m1 * ((dir_m1 * ((nick_m1 * ((nick_m1 * ((nick_m2 * ((4.0 * nick_m1) - (6.0 * nick_m2))) - (nick_m1 * nick_m1))) + (nick_m2 * ((4.0 * (nick_m2 * nick_m2)) + (6.0 * (dir_m2 * dir_m2)) + (dir_m1 * dir_m1))))) - pow(nick_m2, 4.0))) + (nick_m1 * nick_m2 * dir_m2 * ((nick_m1 * ((12.0 * nick_m2) - (8.0 * nick_m1))) - (8.0 * (nick_m2 * nick_m2)))))) + ((dir_m2 * dir_m2) * ((nick_m1 * ((nick_m2 * (dir_m2 * dir_m2)) - pow(nick_m1, 3.0))) - pow(nick_m2, 4.0))))) + ((d_m1_2 + d_m2_2) * ((((nick_m1 * nick_m1) + (nick_m2 * nick_m2)) * ((dir_m1 * dir_m1) + (dir_m2 * dir_m2))) + pow(nick_m1, 4.0) + pow(nick_m2, 4.0))))) + (4.0 * ((dir_m1 * ((nick_m1 * ((nick_m1 * dir_m2 * (((scale * scale) * ((nick_m1 * nick_m1) + (dir_m1 * dir_m1) + (dir_m2 * dir_m2))) - d_m1_2 - d_m2_2)) - (dir_m1 * nick_m2 * (d_m2_2 + d_m1_2)))) + ((nick_m2 * nick_m2) * dir_m2 * (((scale * scale) * ((nick_m2 * nick_m2) + (dir_m1 * dir_m1) + (dir_m2 * dir_m2))) - d_m1_2 - d_m2_2)))) - (nick_m1 * nick_m2 * (dir_m2 * dir_m2) * (d_m1_2 + d_m2_2)))) + ((scale * scale) * ((nick_m1 * ((nick_m1 * ((nick_m2 * (((dir_m2 * dir_m2) * ((8.0 * nick_m1) - (12.0 * nick_m2))) + (nick_m1 * ((nick_m1 * ((6.0 * nick_m1) - (15.0 * nick_m2))) + (20.0 * (nick_m2 * nick_m2)))) - (15.0 * pow(nick_m2, 3.0)))) - ((dir_m1 * dir_m1) * ((dir_m1 * dir_m1) + (6.0 * (dir_m2 * dir_m2)))) - pow(dir_m2, 4.0) - pow(nick_m1, 4.0))) + (nick_m2 * ((8.0 * dir_m2 * ((dir_m2 * ((nick_m2 * nick_m2) - (dir_m1 * dir_m2))) - pow(dir_m1, 3.0))) + (6.0 * pow(nick_m2, 4.0)))))) - ((nick_m2 * nick_m2) * (((dir_m1 * dir_m1) * ((6.0 * (dir_m2 * dir_m2)) + (dir_m1 * dir_m1))) + pow(dir_m2, 4.0) + pow(nick_m2, 4.0))))) + (nick_m1 * nick_m2 * (d_m1_2 + d_m2_2) * ((8.0 * ((dir_m1 * dir_m2) - (nick_m1 * nick_m1) - (nick_m2 * nick_m2))) + (12.0 * nick_m1 * nick_m2))))) + (((d_m1_2 * ((2.0 * d_m2_2) - d_m1_2)) - (d_m2_2 * d_m2_2)) * ((nick_m1 * (nick_m1 - (2.0 * nick_m2))) + (nick_m2 * nick_m2))))), (1.0 / 2.0)) * sign3 / 2.0) + (2.0 * (scale * scale) * (((scale * scale) * ((dir_m1 * ((nick_m1 * (nick_m1 - (2.0 * nick_m2))) + (nick_m2 * nick_m2) - (dir_m2 * (dir_m2 + dir_m1)) + (dir_m1 * dir_m1))) + (dir_m2 * ((dir_m2 * dir_m2) + (nick_m1 * (nick_m1 - (2.0 * nick_m2))) + (nick_m2 * nick_m2))))) + ((dir_m2 - dir_m1) * (d_m1_2 - d_m2_2))))) / (4.0 * pow(scale, 4.0) * ((dir_m1 * (dir_m1 - (2.0 * dir_m2))) + (dir_m2 * dir_m2) + (nick_m1 * (nick_m1 - (2.0 * nick_m2))) + (nick_m2 * nick_m2))));
return dir_view;
}
double
ProjectionSphaeric::comp_nick_view(const Hill *m1, const Hill *m2, double d_m1_2, double scale, double dir_view, double sign1) {
double nick_view;
double nick_m1 = m1->a_nick;
double dir_m1 = m1->alph;
nick_view = ((pow(((d_m1_2 / (scale * scale)) - pow((dir_view - dir_m1), 2.0)), (1.0 / 2.0)) * sign1) + nick_m1);
return nick_view;
}
double
ProjectionSphaeric::comp_tilt_view(const Hill *m, double scale, double dir_view, double nick_view) {
double sin_a_tilt1, sin_a_tilt2, sin_a_tilt, res;
double x_tmp = (m->alph - dir_view) * scale;
double y_tmp = (m->a_nick - nick_view) * scale;
double x = m->x;
double y = m->y;
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;
}
|
