/[zanavi_public1]/navit/navit/transform.c
ZANavi

Contents of /navit/navit/transform.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 50 - (show annotations) (download)
Wed Jun 22 07:33:35 2016 UTC (7 years, 9 months ago) by zoff99
File MIME type: text/plain
File size: 44702 byte(s)
v2.0.51
1 /**
2 * ZANavi, Zoff Android Navigation system.
3 * Copyright (C) 2011-2014 Zoff <zoff@zoff.cc>
4 *
5 * This program is free software; you can redistribute it and/or
6 * modify it under the terms of the GNU General Public License
7 * version 2 as published by the Free Software Foundation.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the
16 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
17 * Boston, MA 02110-1301, USA.
18 */
19
20 /**
21 * Navit, a modular navigation system.
22 * Copyright (C) 2005-2008 Navit Team
23 *
24 * This program is free software; you can redistribute it and/or
25 * modify it under the terms of the GNU General Public License
26 * version 2 as published by the Free Software Foundation.
27 *
28 * This program is distributed in the hope that it will be useful,
29 * but WITHOUT ANY WARRANTY; without even the implied warranty of
30 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 * GNU General Public License for more details.
32 *
33 * You should have received a copy of the GNU General Public License
34 * along with this program; if not, write to the
35 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
36 * Boston, MA 02110-1301, USA.
37 */
38
39 #define _USE_MATH_DEFINES 1
40 #include <assert.h>
41 #include <stdio.h>
42 #include <math.h>
43 #include <limits.h>
44 #include <glib.h>
45 #include <string.h>
46 #include "config.h"
47 #include "coord.h"
48 #include "debug.h"
49 #include "item.h"
50 #include "map.h"
51 #include "transform.h"
52 #include "projection.h"
53 #include "point.h"
54 #include "navit.h"
55
56 #define POST_SHIFT 8
57
58 #ifdef ENABLE_ROLL
59 #define HOG(t) ((t).hog)
60 #else
61 #define HOG(t) 0
62 #endif
63
64
65
66
67
68
69
70 /*
71 * ---------------------------------------------
72 *
73 * single point of GEO transformation functions
74 *
75 * ---------------------------------------------
76 */
77
78 /* ZZ GEO TRANS ZZ */
79
80 #ifndef NAVIT_TRANS_LAT_LON_GEO_NOFUNCS
81
82 __inline__ double transform_to_geo_lat(int y)
83 {
84 double lat = TO_GEO_LAT_(y);
85 return lat;
86 }
87
88 __inline__ double transform_to_geo_lon(int x)
89 {
90 double lon = TO_GEO_LON_(x);
91 return lon;
92 }
93
94 __inline__ int transform_from_geo_lat(double lat)
95 {
96 int y = FROM_GEO_LAT_(lat);
97
98 return y;
99 }
100
101 __inline__ int transform_from_geo_lon(double lon)
102 {
103 int x = FROM_GEO_LON_(lon);
104
105 return x;
106 }
107
108
109 // --- FASTER ---
110 // --- FASTER ---
111 __inline__ double transform_to_geo_lat_fast(int y)
112 {
113 double lat = TO_GEO_LAT_FAST_(y);
114 return lat;
115 }
116
117 __inline__ double transform_to_geo_lon_fast(int x)
118 {
119 double lon = TO_GEO_LON_FAST_(x);
120 return lon;
121 }
122
123 __inline__ int transform_from_geo_lat_fast(double lat)
124 {
125 int y = FROM_GEO_LAT_FAST_(lat);
126
127 return y;
128 }
129
130 __inline__ int transform_from_geo_lon_fast(double lon)
131 {
132 int x = FROM_GEO_LON_FAST_(lon);
133
134 return x;
135 }
136 // --- FASTER ---
137 // --- FASTER ---
138
139 #endif
140
141 /*
142 * ---------------------------------------------
143 *
144 * single point of GEO transformation functions
145 *
146 * ---------------------------------------------
147 */
148
149
150
151
152 static void transform_set_screen_dist(struct transformation *t, int dist)
153 {
154 t->screen_dist = dist;
155 t->xscale3d = dist;
156 t->yscale3d = dist;
157 t->wscale3d = dist << POST_SHIFT;
158 }
159
160 static void transform_setup_matrix(struct transformation *t)
161 {
162 navit_float det;
163 navit_float fac;
164
165 /* this is turning (direction!) ----- */
166 navit_float yawc = navit_cos(-M_PI * t->yaw / 180);
167 navit_float yaws = navit_sin(-M_PI * t->yaw / 180);
168
169 dbg(0, "yaw=%d\n", t->yaw);
170 dbg(0, "yawc=%f\n", yawc);
171 dbg(0, "yaws=%f\n", yaws);
172 /* this is turning (direction!) ----- */
173
174 navit_float pitchc = navit_cos(-M_PI * t->pitch / 180);
175 navit_float pitchs = navit_sin(-M_PI * t->pitch / 180);
176
177 dbg(0, "pitch=%d\n", t->pitch);
178 dbg(0, "pitchc=%f\n", pitchc);
179 dbg(0, "pitchs=%f\n", pitchs);
180
181 #ifdef ENABLE_ROLL
182 navit_float rollc=navit_cos(M_PI*t->roll/180);
183 navit_float rolls=navit_sin(M_PI*t->roll/180);
184 #else
185 navit_float rollc = 1;
186 navit_float rolls = 0;
187 #endif
188
189 int scale = t->scale;
190 int order_dir = -1;
191
192 //dbg(1,"yaw=%d pitch=%d center=0x%x,0x%x\n", t->yaw, t->pitch, t->map_center.x, t->map_center.y);
193 t->znear = 1 << POST_SHIFT;
194 t->zfar = 300 * t->znear;
195 t->scale_shift = 0;
196 t->order = t->order_base;
197
198 if (t->scale >= 1)
199 {
200 scale = t->scale;
201 }
202 else
203 {
204 scale = 1.0 / t->scale;
205 order_dir = 1;
206 }
207
208 while (scale > 1)
209 {
210 if (order_dir < 0)
211 {
212 t->scale_shift++;
213 }
214
215 t->order += order_dir;
216 scale >>= 1;
217 }
218
219 // fac = 256 * 1 / t->scale
220 fac = (1 << POST_SHIFT) * (1 << t->scale_shift) / t->scale;
221 dbg(0,"fac=%f = (1_post_shift=%d) * (1_scale_shift=%d) / (scale=%f) ## order=%d\n", fac, (1 << POST_SHIFT), (1 << t->scale_shift), t->scale, t->order);
222
223 // -----------------------------
224 // rollc = 1, rolls = 0
225 t->m00 = rollc * yawc * fac;
226 t->m01 = rollc * yaws * fac;
227 t->m02 = -rolls * fac;
228 t->m10 = (pitchs * rolls * yawc - pitchc * yaws) * (-fac);
229 t->m11 = (pitchs * rolls * yaws + pitchc * yawc) * (-fac);
230 t->m12 = pitchs * rollc * (-fac);
231 t->m20 = (pitchc * rolls * yawc + pitchs * yaws) * fac;
232 t->m21 = (pitchc * rolls * yaws - pitchs * yawc) * fac;
233 t->m22 = pitchc * rollc * fac;
234
235 // HOG(x) = 0 !!
236 // m00 = fac
237 // m01 = 0
238 // m10 = 0
239 // m11 = -fac
240
241 // ** // xcn = xc * t->m00 + yc * t->m01 + HOG(*t) * t->m02;
242 // ** // ycn = xc * t->m10 + yc * t->m11 + HOG(*t) * t->m12;
243
244 // -----------------------------
245
246 t->offx = t->screen_center.x;
247 t->offy = t->screen_center.y;
248
249 if (t->pitch)
250 {
251 t->ddd = 1;
252 t->offz = t->screen_dist;
253 //dbg(1,"near %d far %d\n",t->znear,t->zfar);
254 t->xscale = t->xscale3d;
255 t->yscale = t->yscale3d;
256 t->wscale = t->wscale3d;
257 }
258 else
259 {
260 t->ddd = 0;
261 t->offz = 0;
262 t->xscale = 1;
263 t->yscale = 1;
264 t->wscale = 1;
265 }
266
267 det = (navit_float) t->m00 * (navit_float) t->m11 * (navit_float) t->m22 + (navit_float) t->m01 * (navit_float) t->m12 * (navit_float) t->m20 + (navit_float) t->m02 * (navit_float) t->m10 * (navit_float) t->m21 - (navit_float) t->m02 * (navit_float) t->m11 * (navit_float) t->m20 - (navit_float) t->m01 * (navit_float) t->m10 * (navit_float) t->m22 - (navit_float) t->m00 * (navit_float) t->m12 * (navit_float) t->m21;
268
269 t->im00 = (t->m11 * t->m22 - t->m12 * t->m21) / det;
270 t->im01 = (t->m02 * t->m21 - t->m01 * t->m22) / det;
271 t->im02 = (t->m01 * t->m12 - t->m02 * t->m11) / det;
272 t->im10 = (t->m12 * t->m20 - t->m10 * t->m22) / det;
273 t->im11 = (t->m00 * t->m22 - t->m02 * t->m20) / det;
274 t->im12 = (t->m02 * t->m10 - t->m00 * t->m12) / det;
275 t->im20 = (t->m10 * t->m21 - t->m11 * t->m20) / det;
276 t->im21 = (t->m01 * t->m20 - t->m00 * t->m21) / det;
277 t->im22 = (t->m00 * t->m11 - t->m01 * t->m10) / det;
278 }
279
280 struct transformation *
281 transform_new(void)
282 {
283 struct transformation *this_;
284
285 this_=g_new0(struct transformation, 1);
286 transform_set_screen_dist(this_, 100);
287 this_->order_base = 14;
288 #if 0
289 this_->pitch=20;
290 #endif
291 #if 0
292 this_->roll=30;
293 this_->hog=1000;
294 #endif
295 transform_setup_matrix(this_);
296 return this_;
297 }
298
299 int transform_get_hog(struct transformation *this_)
300 {
301 return HOG(*this_);
302 }
303
304 void transform_set_hog(struct transformation *this_, int hog)
305 {
306 #ifdef ENABLE_ROLL
307 this_->hog=hog;
308 #else
309 dbg(0, "not supported\n");
310 #endif
311
312 }
313
314 int transform_get_attr(struct transformation *this_, enum attr_type type, struct attr *attr, struct attr_iter *iter)
315 {
316 switch (type)
317 {
318 #ifdef ENABLE_ROLL
319 case attr_hog:
320 attr->u.num=this_->hog;
321 break;
322 #endif
323 default:
324 return 0;
325 }
326
327 attr->type = type;
328
329 return 1;
330 }
331
332 int transform_set_attr(struct transformation *this_, struct attr *attr)
333 {
334 switch (attr->type)
335 {
336 #ifdef ENABLE_ROLL
337 case attr_hog:
338 this_->hog=attr->u.num;
339 return 1;
340 #endif
341 default:
342 return 0;
343 }
344 }
345
346 int transformation_get_order_base(struct transformation *this_)
347 {
348 return this_->order_base;
349 }
350
351 void transform_set_order_base(struct transformation *this_, int order_base)
352 {
353 this_->order_base = order_base;
354 }
355
356 struct transformation *
357 transform_dup(struct transformation *t)
358 {
359 struct transformation *ret=g_new0(struct transformation, 1);
360 *ret = *t;
361 return ret;
362 }
363
364 static const navit_float gar2geo_units = 360.0 / (1 << 24);
365 static const navit_float geo2gar_units = 1 / (360.0 / (1 << 24));
366
367 void transform_to_geo(enum projection pro, struct coord *c, struct coord_geo *g)
368 {
369 // dbg(0,"enter\n");
370
371 int x, y, northern, zone;
372
373 #if 0
374 int hash_id;
375 int s;
376 struct hash_entry_transform *v = NULL;
377 #endif
378
379 switch (pro)
380 {
381 case projection_mg:
382 /* ZZ GEO TRANS ZZ */
383 g->lng = TO_GEO_LON_FAST_(c->x);
384 g->lat = TO_GEO_LAT_FAST_(c->y);
385 /* ZZ GEO TRANS ZZ */
386 break;
387 case projection_garmin:
388 g->lng = c->x * gar2geo_units;
389 g->lat = c->y * gar2geo_units;
390 break;
391 case projection_utm:
392 x = c->x;
393 y = c->y;
394 northern = y >= 0;
395 if (!northern)
396 {
397 y += 10000000;
398 }
399 zone = (x / 1000000);
400 x = x % 1000000;
401 transform_utm_to_geo(x, y, zone, northern, g);
402 break;
403 default:
404 break;
405 }
406 }
407
408 void transform_from_geo(enum projection pro, struct coord_geo *g, struct coord *c)
409 {
410 // dbg(0,"enter\n");
411
412 #if 0
413 int hash_id;
414 int s;
415 struct hash_entry_transform *v = NULL;
416 #endif
417
418 switch (pro)
419 {
420 case projection_mg:
421 /* ZZ GEO TRANS ZZ */
422 c->x = FROM_GEO_LON_FAST_(g->lng);
423 c->y = FROM_GEO_LAT_FAST_(g->lat);
424 /* ZZ GEO TRANS ZZ */
425 break;
426 case projection_garmin:
427 c->x = g->lng * geo2gar_units;
428 c->y = g->lat * geo2gar_units;
429 break;
430 default:
431 break;
432 }
433 }
434
435 void transform_from_to_count(struct coord *cfrom, enum projection from, struct coord *cto, enum projection to, int count)
436 {
437 struct coord_geo g;
438 int i;
439
440 for (i = 0; i < count; i++)
441 {
442 transform_to_geo(from, cfrom, &g);
443 transform_from_geo(to, &g, cto);
444 cfrom++;
445 cto++;
446 }
447 }
448
449 void transform_from_to(struct coord *cfrom, enum projection from, struct coord *cto, enum projection to)
450 {
451 struct coord_geo g;
452
453 transform_to_geo(from, cfrom, &g);
454 transform_from_geo(to, &g, cto);
455 }
456
457 void transform_geo_to_cart(struct coord_geo *geo, navit_float a, navit_float b, struct coord_geo_cart *cart)
458 {
459 navit_float n, ee = 1 - b * b / (a * a);
460 n = a / sqrtf(1 - ee * navit_sin(geo->lat) * navit_sin(geo->lat));
461 cart->x = n * navit_cos(geo->lat) * navit_cos(geo->lng);
462 cart->y = n * navit_cos(geo->lat) * navit_sin(geo->lng);
463 cart->z = n * (1 - ee) * navit_sin(geo->lat);
464 }
465
466 void transform_cart_to_geo(struct coord_geo_cart *cart, navit_float a, navit_float b, struct coord_geo *geo)
467 {
468 navit_float lat, lati, n, ee = 1 - b * b / (a * a), lng = navit_tan(cart->y / cart->x);
469
470 lat = navit_tan(cart->z / navit_sqrt((cart->x * cart->x) + (cart->y * cart->y)));
471 do
472 {
473 lati = lat;
474
475 n = a / navit_sqrt(1 - ee * navit_sin(lat) * navit_sin(lat));
476 lat = navit_atan((cart->z + ee * n * navit_sin(lat)) / navit_sqrt(cart->x * cart->x + cart->y * cart->y));
477 }
478 while (fabs(lat - lati) >= 0.000000000000001);
479
480 geo->lng = lng / M_PI * 180;
481 geo->lat = lat / M_PI * 180;
482 }
483
484 void transform_utm_to_geo(const double UTMEasting, const double UTMNorthing, int ZoneNumber, int NorthernHemisphere, struct coord_geo *geo)
485 {
486 //converts UTM coords to lat/long. Equations from USGS Bulletin 1532
487 //East Longitudes are positive, West longitudes are negative.
488 //North latitudes are positive, South latitudes are negative
489 //Lat and Long are in decimal degrees.
490 //Written by Chuck Gantz- chuck.gantz@globalstar.com
491
492 double Lat, Long;
493 double k0 = 0.99960000000000004;
494 double a = 6378137;
495 double eccSquared = 0.0066943799999999998;
496 double eccPrimeSquared;
497 double e1 = (1 - sqrt(1 - eccSquared)) / (1 + sqrt(1 - eccSquared));
498 double N1, T1, C1, R1, D, M;
499 double LongOrigin;
500 double mu, phi1, phi1Rad;
501 double x, y;
502 double rad2deg = 180 / M_PI;
503
504 x = UTMEasting - 500000.0; //remove 500,000 meter offset for longitude
505 y = UTMNorthing;
506
507 if (!NorthernHemisphere)
508 {
509 y -= 10000000.0;//remove 10,000,000 meter offset used for southern hemisphere
510 }
511
512 LongOrigin = (ZoneNumber - 1) * 6 - 180 + 3; //+3 puts origin in middle of zone
513
514 eccPrimeSquared = (eccSquared) / (1 - eccSquared);
515
516 M = y / k0;
517 mu = M / (a * (1 - eccSquared / 4 - 3 * eccSquared * eccSquared / 64 - 5 * eccSquared * eccSquared * eccSquared / 256));
518 phi1Rad = mu + (3 * e1 / 2 - 27 * e1 * e1 * e1 / 32) * sin(2 * mu) + (21 * e1 * e1 / 16 - 55 * e1 * e1 * e1 * e1 / 32) * sin(4 * mu) + (151 * e1 * e1 * e1 / 96) * sin(6 * mu);
519 phi1 = phi1Rad * rad2deg;
520
521 N1 = a / sqrt(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad));
522 T1 = tan(phi1Rad) * tan(phi1Rad);
523 C1 = eccPrimeSquared * cos(phi1Rad) * cos(phi1Rad);
524 R1 = a * (1 - eccSquared) / pow(1 - eccSquared * sin(phi1Rad) * sin(phi1Rad), 1.5);
525 D = x / (N1 * k0);
526
527 Lat = phi1Rad - (N1 * tan(phi1Rad) / R1) * (D * D / 2 - (5 + 3 * T1 + 10 * C1 - 4 * C1 * C1 - 9 * eccPrimeSquared) * D * D * D * D / 24 + (61 + 90 * T1 + 298 * C1 + 45 * T1 * T1 - 252 * eccPrimeSquared - 3 * C1 * C1) * D * D * D * D * D * D / 720);
528 Lat = Lat * rad2deg;
529
530 Long = (D - (1 + 2 * T1 + C1) * D * D * D / 6 + (5 - 2 * C1 + 28 * T1 - 3 * C1 * C1 + 8 * eccPrimeSquared + 24 * T1 * T1) * D * D * D * D * D / 120) / cos(phi1Rad);
531 Long = LongOrigin + Long * rad2deg;
532
533 geo->lat = Lat;
534 geo->lng = Long;
535 }
536
537 void transform_datum(struct coord_geo *from, enum map_datum from_datum, struct coord_geo *to, enum map_datum to_datum)
538 {
539 }
540
541 int transform(struct transformation *t, enum projection pro, struct coord *c, struct point *p, int count, int mindist, int width, int *width_return)
542 {
543 struct coord c1;
544 int xcn, ycn;
545 struct coord_geo g;
546 int xc, yc, zc = 0, xco = 0, yco = 0, zco = 0;
547 int xm, ym, zct;
548 int zlimit = t->znear;
549 int visible, visibleo = -1;
550 int i, j = 0, k = 0;
551
552 /* ZZ GEO PX ZZ */
553 int mindist2 = TO_SCREEN_(mindist);
554 /* ZZ GEO PX ZZ */
555
556 //dbg(0,"count=%d\n", count);
557 for (i = 0; i < count; i++) // how many coords to calculate?
558 {
559 if (pro == t->pro)
560 {
561 /* ZZ GEO PX ZZ */
562 xc = TO_SCREEN_(c[i].x);
563 yc = TO_SCREEN_(c[i].y);
564 /* ZZ GEO PX ZZ */
565 }
566 else
567 {
568 // if not in "t->pro" (usually "projection_mg") than calc to geo and back to "t->pro" (usually "projection_mg")
569
570 //dbg(0,"to from geo\n");
571 transform_to_geo(pro, &c[i], &g);
572 transform_from_geo(t->pro, &g, &c1);
573 /* ZZ GEO PX ZZ */
574 xc = TO_SCREEN_(c1.x);
575 yc = TO_SCREEN_(c1.y);
576 /* ZZ GEO PX ZZ */
577 }
578
579
580
581 // if next coord is closer than "mindist" to prev coord -> leave it out, continute to next coord ----------------------
582 if (i != 0 && i != count - 1 && mindist2)
583 {
584 /* ZZ GEO PX ZZ */
585 if (xc > TO_SCREEN_(c[k].x) - mindist2 && xc < TO_SCREEN_(c[k].x) + mindist2 && yc > TO_SCREEN_(c[k].y) - mindist2
586 && yc < TO_SCREEN_(c[k].y) + mindist2 && (c[i + 1].x != c[0].x || c[i + 1].y != c[0].y))
587 {
588 continue;
589 }
590 k = i;
591 }
592 // if next coord is closer than "mindist2" to prev coord -> leave it out, continute to next coord ----------------------
593
594
595 #if 0
596 // useless !?
597 xm = xc;
598 ym = yc;
599 // useless !?
600 #endif
601
602 // dbg(2,"0x%x, 0x%x - 0x%x,0x%x contains 0x%x,0x%x\n", t->r.lu.x, t->r.lu.y, t->r.rl.x, t->r.rl.y, c->x, c->y);
603 // ret=coord_rect_contains(&t->r, c);
604
605 /* ZZ GEO PX ZZ */
606 xc -= TO_SCREEN_(t->map_center.x); // relative to map center
607 yc -= TO_SCREEN_(t->map_center.y); // relative to map center
608 /* ZZ GEO PX ZZ */
609
610 xc >>= t->scale_shift; // apply zoom level
611 yc >>= t->scale_shift; // apply zoom level
612
613 // ------------------------
614 // xm, ym never used?
615 xm = xc;
616 ym = yc;
617 // xm, ym never used?
618 // ------------------------
619
620
621 // -- Matrix transform --
622 #if 0
623 xcn = xc * t->m00 + yc * t->m01 + HOG(*t) * t->m02;
624 ycn = xc * t->m10 + yc * t->m11 + HOG(*t) * t->m12;
625 #endif
626
627 #if 1
628 xcn = xc * t->m00 + yc * t->m01;
629 ycn = xc * t->m10 + yc * t->m11;
630 #endif
631
632 //dbg(0, "cxn = xc * t->m00 + yc * t->m01 # %d %d %d %d %d\n", xcn, xc, t->m00, yc, t->m01);
633 //dbg(0, "ycn = xc * t->m10 + yc * t->m11 # %d %d %d %d %d\n", ycn, xc, t->m10, yc, t->m11);
634 // -- Matrix transform --
635
636
637
638 #if 0
639 // -------- NEVER USED NOW !! -----------
640 // -------- NEVER USED NOW !! -----------
641 // -------- NEVER USED NOW !! -----------
642 // -------- NEVER USED NOW !! -----------
643 if (t->ddd)
644 {
645 zc = (xc * t->m20 + yc * t->m21 + HOG(*t) * t->m22);
646 zct = zc;
647 zc += t->offz << POST_SHIFT;
648 //dbg(1,"zc=%d\n", zc);
649 //dbg(1,"zc(%d)=xc(%d)*m20(%d)+yc(%d)*m21(%d)\n", (xc*t->m20+yc*t->m21), xc, t->m20, yc, t->m21);
650 /* visibility */
651 visible = (zc < zlimit ? 0 : 1);
652 //dbg(1,"visible=%d old %d\n", visible, visibleo);
653 if (visible != visibleo && visibleo != -1)
654 {
655 //dbg(1,"clipping (%d,%d,%d)-(%d,%d,%d) (%d,%d,%d)\n", xcn, ycn, zc, xco, yco, zco, xco-xcn, yco-ycn, zco-zc);
656 if (zco != zc)
657 {
658 xcn = xcn + (long long) (xco - xcn) * (zlimit - zc) / (zco - zc);
659 ycn = ycn + (long long) (yco - ycn) * (zlimit - zc) / (zco - zc);
660 }
661 //dbg(1,"result (%d,%d,%d) * %d / %d\n", xcn,ycn,zc,zlimit-zc,zco-zc);
662 zc = zlimit;
663 xco = xcn;
664 yco = ycn;
665 zco = zc;
666 if (visible)
667 {
668 i--;
669 }
670 visibleo = visible;
671 }
672 else
673 {
674 xco = xcn;
675 yco = ycn;
676 zco = zc;
677 visibleo = visible;
678
679 if (!visible)
680 {
681 continue;
682 }
683 }
684 //dbg(1,"zc=%d\n", zc);
685 //dbg(1,"xcn %d ycn %d\n", xcn, ycn);
686 //dbg(1,"%d,%d %d\n",xc,yc,zc);
687 //#if 0
688 // dbg(0,"%d/%d=%d %d/%d=%d\n",xcn,xc,xcn/xc,ycn,yc,ycn/yc);
689 //#endif
690
691 //#if 1
692 xc = (long long) xcn * t->xscale / zc;
693 yc = (long long) ycn * t->yscale / zc;
694 //#else
695 // xc=xcn/(1000+zc);
696 // yc=ycn/(1000+zc);
697 //#endif
698
699
700 //#if 0
701 // dbg(1,"%d,%d %d\n",xc,yc,zc);
702 //#endif
703 }
704 // -------- NEVER USED NOW !! -----------
705 // -------- NEVER USED NOW !! -----------
706 // -------- NEVER USED NOW !! -----------
707 // -------- NEVER USED NOW !! -----------
708 else
709 #endif
710 {
711
712 // -------- normal 2D -> used now !! -----------
713 // -------- normal 2D -> used now !! -----------
714
715 xc = xcn;
716 yc = ycn;
717 xc >>= POST_SHIFT; // divide by 256 (2^8)
718 yc >>= POST_SHIFT; // divide by 256 (2^8)
719
720 // -------- normal 2D -> used now !! -----------
721 // -------- normal 2D -> used now !! -----------
722
723 }
724
725 xc += t->offx; // relative to screen center (before was map center!)
726 yc += t->offy; // relative to screen center (before was map center!)
727 p[j].x = xc;
728 p[j].y = yc;
729
730 if (width_return)
731 {
732 if (t->ddd)
733 {
734 width_return[j] = width * t->wscale / zc;
735 }
736 else
737 {
738 width_return[j] = width;
739 }
740 }
741
742
743 j++;
744
745 } // ------- END for loop ----------
746
747 return j;
748 }
749
750 static void transform_apply_inverse_matrix(struct transformation *t, struct coord_geo_cart *in, struct coord_geo_cart *out)
751 {
752 out->x = in->x * t->im00 + in->y * t->im01 + in->z * t->im02;
753 out->y = in->x * t->im10 + in->y * t->im11 + in->z * t->im12;
754 out->z = in->x * t->im20 + in->y * t->im21 + in->z * t->im22;
755 }
756
757 static int transform_zplane_intersection(struct coord_geo_cart *p1, struct coord_geo_cart *p2, navit_float z, struct coord_geo_cart *result)
758 {
759 navit_float dividend = z - p1->z;
760 navit_float divisor = p2->z - p1->z;
761 navit_float q;
762 if (!divisor)
763 {
764 if (dividend)
765 return 0; /* no intersection */
766 else
767 return 3; /* identical planes */
768 }
769
770 q = dividend / divisor;
771 result->x = p1->x + q * (p2->x - p1->x);
772 result->y = p1->y + q * (p2->y - p1->y);
773 result->z = z;
774
775 if (q >= 0 && q <= 1)
776 return 1; /* intersection within [p1,p2] */
777
778 return 2; /* intersection without [p1,p2] */
779 }
780
781 static void transform_screen_to_3d(struct transformation *t, struct point *p, navit_float z, struct coord_geo_cart *cg)
782 {
783 double xc, yc;
784 double offz = t->offz << POST_SHIFT;
785
786 xc = p->x - t->offx;
787 yc = p->y - t->offy;
788 cg->x = xc * z / t->xscale;
789 cg->y = yc * z / t->yscale;
790 cg->z = z - offz;
791 }
792
793 static int transform_reverse_near_far(struct transformation *t, struct point *p, struct coord *c, int near, int far)
794 {
795 double xc, yc;
796
797 //dbg(1,"%d,%d\n",p->x,p->y);
798
799 #if 0
800 // NEVER USED NOW !!!!!!!! ----------------
801 // NEVER USED NOW !!!!!!!! ----------------
802 // NEVER USED NOW !!!!!!!! ----------------
803 if (t->ddd)
804 {
805 struct coord_geo_cart nearc, farc, nears, fars, intersection;
806 transform_screen_to_3d(t, p, near, &nearc);
807 transform_screen_to_3d(t, p, far, &farc);
808 transform_apply_inverse_matrix(t, &nearc, &nears);
809 transform_apply_inverse_matrix(t, &farc, &fars);
810 if (transform_zplane_intersection(&nears, &fars, HOG(*t), &intersection) != 1)
811 {
812 return 0;
813 }
814 xc = intersection.x;
815 yc = intersection.y;
816 }
817 // NEVER USED NOW !!!!!!!! ----------------
818 // NEVER USED NOW !!!!!!!! ----------------
819 // NEVER USED NOW !!!!!!!! ----------------
820 else
821 #endif
822 {
823 double xcn, ycn;
824 xcn = p->x - t->offx; // relative to screen center
825 ycn = p->y - t->offy; // relative to screen center
826 xc = (xcn * t->im00 + ycn * t->im01) * (1 << POST_SHIFT);
827 yc = (xcn * t->im10 + ycn * t->im11) * (1 << POST_SHIFT);
828 }
829
830 /* ZZ GEO PX ZZ */
831 c->x = xc * (1 << t->scale_shift) + TO_SCREEN_(t->map_center.x); // zoom level, and relative to map center
832 c->y = yc * (1 << t->scale_shift) + TO_SCREEN_(t->map_center.y); // zoom level, and relative to map center
833 c->x = FROM_SCREEN_(c->x);
834 c->y = FROM_SCREEN_(c->y);
835 /* ZZ GEO PX ZZ */
836
837 return 1;
838 }
839
840 int transform_reverse(struct transformation *t, struct point *p, struct coord *c)
841 {
842 return transform_reverse_near_far(t, p, c, t->znear, t->zfar);
843 }
844
845 enum projection transform_get_projection(struct transformation *this_)
846 {
847 return this_->pro;
848 }
849
850 void transform_set_projection(struct transformation *this_, enum projection pro)
851 {
852 this_->pro = pro;
853 }
854
855 static int min4(int v1, int v2, int v3, int v4)
856 {
857 int res = v1;
858 if (v2 < res)
859 res = v2;
860 if (v3 < res)
861 res = v3;
862 if (v4 < res)
863 res = v4;
864 return res;
865 }
866
867 static int max4(int v1, int v2, int v3, int v4)
868 {
869 int res = v1;
870 if (v2 > res)
871 res = v2;
872 if (v3 > res)
873 res = v3;
874 if (v4 > res)
875 res = v4;
876 return res;
877 }
878
879 struct map_selection *
880 transform_get_selection(struct transformation *this_, enum projection pro, int order)
881 {
882
883 struct map_selection *ret, *curri, *curro;
884 struct coord_geo g;
885
886 ret = map_selection_dup(this_->map_sel);
887 curri = this_->map_sel;
888 curro = ret;
889 while (curri)
890 {
891 if (this_->pro != pro)
892 {
893 transform_to_geo(this_->pro, &curri->u.c_rect.lu, &g);
894 transform_from_geo(pro, &g, &curro->u.c_rect.lu);
895 // dbg(1,"%f,%f", g.lat, g.lng);
896 transform_to_geo(this_->pro, &curri->u.c_rect.rl, &g);
897 transform_from_geo(pro, &g, &curro->u.c_rect.rl);
898 // dbg(1,": - %f,%f\n", g.lat, g.lng);
899 }
900 // dbg(1,"transform rect for %d is %d,%d - %d,%d\n", pro, curro->u.c_rect.lu.x, curro->u.c_rect.lu.y, curro->u.c_rect.rl.x, curro->u.c_rect.rl.y);
901 curro->order += order;
902 #if 0
903 curro->u.c_rect.lu.x-=500;
904 curro->u.c_rect.lu.y+=500;
905 curro->u.c_rect.rl.x+=500;
906 curro->u.c_rect.rl.y-=500;
907 #endif
908 curro->range = item_range_all;
909 curri = curri->next;
910 curro = curro->next;
911 }
912 return ret;
913 }
914
915 struct coord *
916 transform_center(struct transformation *this_)
917 {
918 return &this_->map_center;
919 }
920
921 struct coord *
922 transform_get_center(struct transformation *this_)
923 {
924 return &this_->map_center;
925 }
926
927 void transform_set_center(struct transformation *this_, struct coord *c)
928 {
929 this_->map_center = *c;
930 }
931
932 void transform_set_yaw(struct transformation *t, int yaw)
933 {
934 t->yaw = yaw;
935 transform_setup_matrix(t);
936 }
937
938 int transform_get_yaw(struct transformation *this_)
939 {
940 return this_->yaw;
941 }
942
943 void transform_set_pitch(struct transformation *this_, int pitch)
944 {
945 this_->pitch = pitch;
946 transform_setup_matrix(this_);
947 }
948 int transform_get_pitch(struct transformation *this_)
949 {
950 return this_->pitch;
951 }
952
953 void transform_set_roll(struct transformation *this_, int roll)
954 {
955 #ifdef ENABLE_ROLL
956 this_->roll=roll;
957 transform_setup_matrix(this_);
958 #else
959 //dbg(0, "not supported\n");
960 #endif
961 }
962
963 int transform_get_roll(struct transformation *this_)
964 {
965 #ifdef ENABLE_ROLL
966 return this_->roll;
967 #else
968 return 0;
969 #endif
970 }
971
972 void transform_set_distance(struct transformation *this_, int distance)
973 {
974 transform_set_screen_dist(this_, distance);
975 transform_setup_matrix(this_);
976 }
977
978 int transform_get_distance(struct transformation *this_)
979 {
980 return this_->screen_dist;
981 }
982
983 void transform_set_scales(struct transformation *this_, int xscale, int yscale, int wscale)
984 {
985 this_->xscale3d = xscale;
986 this_->yscale3d = yscale;
987 this_->wscale3d = wscale;
988 }
989
990 void transform_set_screen_selection(struct transformation *t, struct map_selection *sel)
991 {
992 map_selection_destroy(t->screen_sel);
993 t->screen_sel = map_selection_dup(sel);
994 if (sel)
995 {
996 t->screen_center.x = (sel->u.p_rect.rl.x - sel->u.p_rect.lu.x) / 2;
997 t->screen_center.y = (sel->u.p_rect.rl.y - sel->u.p_rect.lu.y) / 2;
998 transform_setup_matrix(t);
999 }
1000 }
1001
1002 void transform_set_screen_center(struct transformation *t, struct point *p)
1003 {
1004 t->screen_center = *p;
1005 }
1006
1007 #if 0
1008 void
1009 transform_set_size(struct transformation *t, int width, int height)
1010 {
1011 t->width=width;
1012 t->height=height;
1013 }
1014 #endif
1015
1016 void transform_get_size(struct transformation *t, int *width, int *height)
1017 {
1018 struct point_rect *r;
1019 if (t->screen_sel)
1020 {
1021 r = &t->screen_sel->u.p_rect;
1022 *width = r->rl.x - r->lu.x;
1023 *height = r->rl.y - r->lu.y;
1024 }
1025 }
1026
1027 void transform_setup(struct transformation *t, struct pcoord *c, int scale, int yaw)
1028 {
1029 t->pro = c->pro;
1030 t->map_center.x = c->x;
1031 t->map_center.y = c->y;
1032 t->scale = scale / 16.0;
1033 transform_set_yaw(t, yaw);
1034 }
1035
1036 #if 0
1037
1038 void
1039 transform_setup_source_rect_limit(struct transformation *t, struct coord *center, int limit)
1040 {
1041 t->center=*center;
1042 t->scale=1;
1043 t->angle=0;
1044 t->r.lu.x=center->x-limit;
1045 t->r.rl.x=center->x+limit;
1046 t->r.rl.y=center->y-limit;
1047 t->r.lu.y=center->y+limit;
1048 }
1049 #endif
1050
1051 void transform_setup_source_rect(struct transformation *t)
1052 {
1053 int i;
1054 struct coord screen[4];
1055 struct point screen_pnt[4];
1056 struct point_rect *pr;
1057 struct map_selection *ms, *msm, *next, **msm_last;
1058 ms = t->map_sel;
1059 while (ms)
1060 {
1061 next = ms->next;
1062 g_free(ms);
1063 ms = next;
1064 }
1065 t->map_sel = NULL;
1066 msm_last = &t->map_sel;
1067 ms = t->screen_sel;
1068 while (ms)
1069 {
1070 msm=g_new0(struct map_selection, 1);
1071 *msm = *ms;
1072 pr = &ms->u.p_rect;
1073 screen_pnt[0].x = pr->lu.x; /* left upper */
1074 screen_pnt[0].y = pr->lu.y;
1075 screen_pnt[1].x = pr->rl.x; /* right upper */
1076 screen_pnt[1].y = pr->lu.y;
1077 screen_pnt[2].x = pr->rl.x; /* right lower */
1078 screen_pnt[2].y = pr->rl.y;
1079 screen_pnt[3].x = pr->lu.x; /* left lower */
1080 screen_pnt[3].y = pr->rl.y;
1081 if (t->ddd)
1082 {
1083 struct coord_geo_cart tmp, cg[8];
1084 struct coord c;
1085 int valid = 0;
1086 unsigned char edgenodes[] = { 0, 1, 1, 2, 2, 3, 3, 0, 4, 5, 5, 6, 6, 7, 7, 4, 0, 4, 1, 5, 2, 6, 3, 7 };
1087 for (i = 0; i < 8; i++)
1088 {
1089 transform_screen_to_3d(t, &screen_pnt[i % 4], (i >= 4 ? t->zfar : t->znear), &tmp);
1090 transform_apply_inverse_matrix(t, &tmp, &cg[i]);
1091 }
1092 msm->u.c_rect.lu.x = 0;
1093 msm->u.c_rect.lu.y = 0;
1094 msm->u.c_rect.rl.x = 0;
1095 msm->u.c_rect.rl.y = 0;
1096 for (i = 0; i < 12; i++)
1097 {
1098 if (transform_zplane_intersection(&cg[edgenodes[i * 2]], &cg[edgenodes[i * 2 + 1]], HOG(*t), &tmp) == 1)
1099 {
1100 c.x = tmp.x * (1 << t->scale_shift) + t->map_center.x;
1101 c.y = tmp.y * (1 << t->scale_shift) + t->map_center.y;
1102 //dbg(1,"intersection with edge %d at 0x%x,0x%x\n",i,c.x,c.y);
1103 if (valid)
1104 coord_rect_extend(&msm->u.c_rect, &c);
1105 else
1106 {
1107 msm->u.c_rect.lu = c;
1108 msm->u.c_rect.rl = c;
1109 valid = 1;
1110 }
1111 //dbg(1,"rect 0x%x,0x%x - 0x%x,0x%x\n",msm->u.c_rect.lu.x,msm->u.c_rect.lu.y,msm->u.c_rect.rl.x,msm->u.c_rect.rl.y);
1112 }
1113 }
1114 }
1115 else
1116 {
1117 for (i = 0; i < 4; i++)
1118 {
1119 transform_reverse(t, &screen_pnt[i], &screen[i]);
1120 //dbg(1,"map(%d) %d,%d=0x%x,0x%x\n", i,screen_pnt[i].x, screen_pnt[i].y, screen[i].x, screen[i].y);
1121 }
1122 msm->u.c_rect.lu.x = min4(screen[0].x, screen[1].x, screen[2].x, screen[3].x);
1123 msm->u.c_rect.rl.x = max4(screen[0].x, screen[1].x, screen[2].x, screen[3].x);
1124 msm->u.c_rect.rl.y = min4(screen[0].y, screen[1].y, screen[2].y, screen[3].y);
1125 msm->u.c_rect.lu.y = max4(screen[0].y, screen[1].y, screen[2].y, screen[3].y);
1126 }
1127 //dbg(1,"%dx%d\n", msm->u.c_rect.rl.x-msm->u.c_rect.lu.x,
1128 // msm->u.c_rect.lu.y-msm->u.c_rect.rl.y);
1129 *msm_last = msm;
1130 msm_last = &msm->next;
1131 ms = ms->next;
1132 }
1133 }
1134
1135 long transform_get_scale(struct transformation *t)
1136 {
1137 return (int) (t->scale * 16);
1138 }
1139
1140 void transform_set_scale(struct transformation *t, long scale)
1141 {
1142 t->scale = scale / 16.0;
1143 transform_setup_matrix(t);
1144 }
1145
1146 int transform_get_order(struct transformation *t)
1147 {
1148 //dbg(1,"order %d\n", t->order);
1149 return t->order;
1150 }
1151
1152 #define TWOPI (M_PI*2)
1153 #define GC2RAD(c) ((c) * TWOPI/(1<<24))
1154 #define minf(a,b) ((a) < (b) ? (a) : (b))
1155
1156 static double transform_distance_garmin(struct coord *c1, struct coord *c2)
1157 {
1158 #ifdef USE_HALVESINE
1159 static const int earth_radius = 6371*1000; //m change accordingly
1160 // static const int earth_radius = 3960; //miles
1161
1162 //Point 1 cords
1163 navit_float lat1 = GC2RAD(c1->y);
1164 navit_float long1 = GC2RAD(c1->x);
1165
1166 //Point 2 cords
1167 navit_float lat2 = GC2RAD(c2->y);
1168 navit_float long2 = GC2RAD(c2->x);
1169
1170 //Haversine Formula
1171 navit_float dlong = long2-long1;
1172 navit_float dlat = lat2-lat1;
1173
1174 navit_float sinlat = navit_sin(dlat/2);
1175 navit_float sinlong = navit_sin(dlong/2);
1176
1177 navit_float a=(sinlat*sinlat)+navit_cos(lat1)*navit_cos(lat2)*(sinlong*sinlong);
1178 navit_float c=2*navit_asin(minf(1,navit_sqrt(a)));
1179 #ifdef AVOID_FLOAT
1180 return round(earth_radius*c);
1181 #else
1182 return earth_radius*c;
1183 #endif
1184 #else
1185 #define GMETER 2.3887499999999999
1186 navit_float dx, dy;
1187 dx = c1->x - c2->x;
1188 dy = c1->y - c2->y;
1189 return navit_sqrt(dx * dx + dy * dy) * GMETER;
1190 #undef GMETER
1191 #endif
1192 }
1193
1194 double transform_scale(int y)
1195 {
1196 struct coord c;
1197 struct coord_geo g;
1198 c.x = 0;
1199 c.y = y;
1200 transform_to_geo(projection_mg, &c, &g);
1201 return 1 / navit_cos(g.lat / 180 * M_PI);
1202 }
1203
1204 #ifdef AVOID_FLOAT
1205 static int
1206 tab_sqrt[]=
1207 { 14142,13379,12806,12364,12018,11741,11517,11333,11180,11051,10943,10850,10770,10701,10640,10587,10540,10499,10462,10429,10400,10373,10349,10327,10307,10289,10273,10257,10243,10231,10219,10208};
1208
1209 static int tab_int_step = 0x20000;
1210
1211 static int tab_int_scale[]=
1212 { 10000,10002,10008,10019,10033,10052,10076,10103,10135,10171,10212,10257,10306,10359,10417,10479,10546,10617,10693,10773,10858,10947,11041,11140,11243,11352,11465,11582,11705,11833,11965,12103,12246,12394,12547,12706,12870,13039,13214,13395,13581,13773,13971,14174,14384,14600,14822,15050,15285,15526,15774,16028,16289,16557,16832,17114,17404,17700,18005,18316,18636,18964,19299,19643,19995,20355,20724,21102,21489,21885,22290,22705,23129,23563,24007,24461,24926,25401,25886,26383,26891};
1213
1214 int transform_int_scale(int y)
1215 {
1216 int i,size = sizeof(tab_int_scale)/sizeof(int);
1217
1218 if (y < 0)
1219 {
1220 y=-y;
1221 }
1222
1223 i=y/tab_int_step;
1224
1225 if (i < size-1)
1226 {
1227 return tab_int_scale[i]+((tab_int_scale[i+1]-tab_int_scale[i])*(y-i*tab_int_step))/tab_int_step;
1228 }
1229
1230 return tab_int_scale[size-1];
1231 }
1232 #endif
1233
1234 double transform_distance(enum projection pro, struct coord *c1, struct coord *c2)
1235 {
1236 if (pro == projection_mg)
1237 {
1238
1239 /* ZZ GEO PX ZZ */
1240 struct coord c1_corr;
1241 struct coord *c1c = &c1_corr;
1242 c1_corr.x = TO_SCREEN_(c1->x);
1243 c1_corr.y = TO_SCREEN_(c1->y);
1244 struct coord c2_corr;
1245 struct coord *c2c = &c2_corr;
1246 c2_corr.x = TO_SCREEN_(c2->x);
1247 c2_corr.y = TO_SCREEN_(c2->y);
1248 /* ZZ GEO PX ZZ */
1249
1250 #ifndef AVOID_FLOAT
1251 double dx, dy, scale = transform_scale((c1c->y + c2c->y) / 2);
1252 dx = c1c->x - c2c->x;
1253 dy = c1c->y - c2c->y;
1254 return sqrt(dx * dx + dy * dy) / scale;
1255 #else
1256 int dx,dy,f,scale=transform_int_scale((c1c->y+c2c->y)/2);
1257 dx=c1c->x-c2c->x;
1258 dy=c1c->y-c2c->y;
1259
1260 if (dx < 0)
1261 dx=-dx;
1262
1263 if (dy < 0)
1264 dy=-dy;
1265
1266 while (dx > 20000 || dy > 20000)
1267 {
1268 dx/=10;
1269 dy/=10;
1270 scale/=10;
1271 }
1272
1273 if (! dy)
1274 return dx*10000/scale;
1275
1276 if (! dx)
1277 return dy*10000/scale;
1278
1279 if (dx > dy)
1280 {
1281 f=dx*8/dy-8;
1282
1283 if (f >= 32)
1284 return dx*10000/scale;
1285
1286 return dx*tab_sqrt[f]/scale;
1287 }
1288 else
1289 {
1290 f=dy*8/dx-8;
1291 if (f >= 32)
1292 return dy*10000/scale;
1293
1294 return dy*tab_sqrt[f]/scale;
1295 }
1296 #endif
1297 }
1298 else if (pro == projection_garmin)
1299 {
1300 return transform_distance_garmin(c1, c2);
1301 }
1302 else
1303 {
1304 dbg(0, "Unknown projection: %d\n", pro);
1305 return 0;
1306 }
1307 }
1308
1309 void transform_project(enum projection pro, struct coord *c, int distance, int angle, struct coord *res)
1310 {
1311 double scale;
1312 switch (pro)
1313 {
1314 case projection_mg:
1315 scale = transform_scale(c->y);
1316 res->x = c->x + distance * sin(angle * M_PI / 180) * scale;
1317 res->y = c->y + distance * cos(angle * M_PI / 180) * scale;
1318 break;
1319 default:
1320 dbg(0, "Unsupported projection: %d\n", pro);
1321 return;
1322 }
1323
1324 }
1325
1326 double transform_polyline_length(enum projection pro, struct coord *c, int count)
1327 {
1328 double ret = 0;
1329 int i;
1330
1331 for (i = 0; i < count - 1; i++)
1332 {
1333 ret += transform_distance(pro, &c[i], &c[i + 1]);
1334 }
1335
1336 return ret;
1337 }
1338
1339 // calc the distance (squared) of a point (p) to a line segment (l1 .. l2)
1340 // return (int) distance squared
1341 int transform_distance_point2line_sq(struct coord *p, struct coord *l1, struct coord *l2)
1342 {
1343 /* ZZ GEO PX ZZ */
1344 struct coord p_corr;
1345 struct coord *pc = &p_corr;
1346 p_corr.x = TO_SCREEN_(p->x);
1347 p_corr.y = TO_SCREEN_(p->y);
1348 struct coord l1_corr;
1349 struct coord *l1c = &l1_corr;
1350 l1_corr.x = TO_SCREEN_(l1->x);
1351 l1_corr.y = TO_SCREEN_(l1->y);
1352 struct coord l2_corr;
1353 struct coord *l2c = &l2_corr;
1354 l2_corr.x = TO_SCREEN_(l2->x);
1355 l2_corr.y = TO_SCREEN_(l2->y);
1356 /* ZZ GEO PX ZZ */
1357
1358
1359 int A = pc->x - l1c->x;
1360 int B = pc->y - l1c->y;
1361 float C = l2c->x - l1c->x;
1362 float D = l2c->y - l1c->y;
1363
1364 int dot = A * C + B * D;
1365 int len_sq = C * C + D * D;
1366 float param = (float) dot / (float) len_sq;
1367
1368 int xx, yy;
1369
1370 if (param < 0 || (l1c->x == l2c->x && l1c->y == l2c->y))
1371 {
1372 xx = l1c->x;
1373 yy = l1c->y;
1374 }
1375 else if (param > 1)
1376 {
1377 xx = l2c->x;
1378 yy = l2c->y;
1379 }
1380 else
1381 {
1382 xx = l1c->x + param * C;
1383 yy = l1c->y + param * D;
1384 }
1385
1386 int dx = pc->x - xx;
1387 int dy = pc->y - yy;
1388
1389 if (dx > 32767 || dy > 32767 || dx < -32767 || dy < -32767)
1390 {
1391 return INT_MAX;
1392 }
1393
1394 return (dx * dx + dy * dy);
1395
1396 }
1397
1398 int transform_distance_sq(struct coord *c1, struct coord *c2)
1399 {
1400 /* ZZ GEO PX ZZ */
1401 struct coord c1_corr;
1402 struct coord *c1c = &c1_corr;
1403 c1_corr.x = TO_SCREEN_(c1->x);
1404 c1_corr.y = TO_SCREEN_(c1->y);
1405 struct coord c2_corr;
1406 struct coord *c2c = &c2_corr;
1407 c2_corr.x = TO_SCREEN_(c2->x);
1408 c2_corr.y = TO_SCREEN_(c2->y);
1409 /* ZZ GEO PX ZZ */
1410
1411 int dx = c1c->x - c2c->x;
1412 int dy = c1c->y - c2c->y;
1413
1414 if (dx > 32767 || dy > 32767 || dx < -32767 || dy < -32767)
1415 {
1416 return INT_MAX;
1417 }
1418 else
1419 {
1420 return dx * dx + dy * dy;
1421 }
1422 }
1423
1424 navit_float transform_distance_sq_float(struct coord *c1, struct coord *c2)
1425 {
1426 /* ZZ GEO PX ZZ */
1427 struct coord c1_corr;
1428 struct coord *c1c = &c1_corr;
1429 c1_corr.x = TO_SCREEN_(c1->x);
1430 c1_corr.y = TO_SCREEN_(c1->y);
1431 struct coord c2_corr;
1432 struct coord *c2c = &c2_corr;
1433 c2_corr.x = TO_SCREEN_(c2->x);
1434 c2_corr.y = TO_SCREEN_(c2->y);
1435 /* ZZ GEO PX ZZ */
1436
1437
1438 int dx = c1c->x - c2c->x;
1439 int dy = c1c->y - c2c->y;
1440 return (navit_float) dx * dx + dy * dy;
1441 }
1442
1443 int transform_distance_sq_pc(struct pcoord *c1, struct pcoord *c2)
1444 {
1445 struct coord p1, p2;
1446 p1.x = c1->x;
1447 p1.y = c1->y;
1448 p2.x = c2->x;
1449 p2.y = c2->y;
1450 return transform_distance_sq(&p1, &p2);
1451 }
1452
1453 int transform_distance_line_sq(struct coord *l0, struct coord *l1, struct coord *ref, struct coord *lpnt)
1454 {
1455 int vx, vy, wx, wy;
1456 int c1, c2;
1457 int climit = 1000000;
1458 struct coord l;
1459
1460 vx = l1->x - l0->x;
1461 vy = l1->y - l0->y;
1462 wx = ref->x - l0->x;
1463 wy = ref->y - l0->y;
1464
1465 c1 = vx * wx + vy * wy;
1466
1467 if (c1 <= 0)
1468 {
1469 if (lpnt)
1470 *lpnt = *l0;
1471
1472 return transform_distance_sq(l0, ref);
1473 }
1474
1475 c2 = vx * vx + vy * vy;
1476
1477 if (c2 <= c1)
1478 {
1479 if (lpnt)
1480 *lpnt = *l1;
1481
1482 return transform_distance_sq(l1, ref);
1483 }
1484
1485 while (c1 > climit || c2 > climit)
1486 {
1487 c1 /= 256;
1488 c2 /= 256;
1489 }
1490
1491 l.x = l0->x + vx * c1 / c2;
1492 l.y = l0->y + vy * c1 / c2;
1493
1494 if (lpnt)
1495 *lpnt = l;
1496
1497 return transform_distance_sq(&l, ref);
1498 }
1499
1500 navit_float transform_distance_line_sq_float(struct coord *l0, struct coord *l1, struct coord *ref, struct coord *lpnt)
1501 {
1502 navit_float vx, vy, wx, wy;
1503 navit_float c1, c2;
1504 struct coord l;
1505
1506 vx = l1->x - l0->x;
1507 vy = l1->y - l0->y;
1508 wx = ref->x - l0->x;
1509 wy = ref->y - l0->y;
1510
1511 c1 = vx * wx + vy * wy;
1512 if (c1 <= 0)
1513 {
1514 if (lpnt)
1515 *lpnt = *l0;
1516 return transform_distance_sq_float(l0, ref);
1517 }
1518 c2 = vx * vx + vy * vy;
1519 if (c2 <= c1)
1520 {
1521 if (lpnt)
1522 *lpnt = *l1;
1523 return transform_distance_sq_float(l1, ref);
1524 }
1525 l.x = l0->x + vx * c1 / c2;
1526 l.y = l0->y + vy * c1 / c2;
1527 if (lpnt)
1528 *lpnt = l;
1529 return transform_distance_sq_float(&l, ref);
1530 }
1531
1532 int transform_distance_polyline_sq__v2(struct coord *c, int count, struct coord *ref)
1533 {
1534 int i, dist, distn;
1535
1536 if (count < 2)
1537 {
1538 int d;
1539 d = transform_distance_sq(&c[0], ref);
1540 //dbg(0,"d=%d\n", d);
1541 return d;
1542 }
1543
1544 dist = transform_distance_point2line_sq(ref, &c[0], &c[1]);
1545 //dbg(0,"dist1:%d\n", dist);
1546
1547 for (i = 2; i < count; i++)
1548 {
1549 distn = transform_distance_point2line_sq(ref, &c[i - 1], &c[i]);
1550 //dbg(0,"dist2:%d\n", dist);
1551 if (distn < dist)
1552 {
1553 dist = distn;
1554 }
1555 }
1556 //dbg(0,"dist final:%d\n", dist);
1557 return dist;
1558 }
1559
1560 int transform_distance_polyline_sq(struct coord *c, int count, struct coord *ref, struct coord *lpnt, int *pos)
1561 {
1562 int i, dist, distn;
1563 struct coord lp;
1564 if (count < 2)
1565 {
1566 // dbg(0,"1\n");
1567 return INT_MAX;
1568 }
1569 if (pos)
1570 {
1571 *pos = 0;
1572 }
1573
1574 dist = transform_distance_line_sq(&c[0], &c[1], ref, lpnt);
1575 // dbg(0,"dist:%d\n", dist);
1576
1577 for (i = 2; i < count; i++)
1578 {
1579 distn = transform_distance_line_sq(&c[i - 1], &c[i], ref, &lp);
1580 if (distn < dist)
1581 {
1582 dist = distn;
1583 if (lpnt)
1584 {
1585 *lpnt = lp;
1586 }
1587 if (pos)
1588 {
1589 *pos = i - 1;
1590 }
1591 }
1592 }
1593 return dist;
1594 }
1595
1596 int transform_douglas_peucker(struct coord *in, int count, int dist_sq, struct coord *out)
1597 {
1598 int ret = 0;
1599 int i, d, dmax = 0, idx = 0;
1600 for (i = 1; i < count - 2; i++)
1601 {
1602 d = transform_distance_line_sq(&in[0], &in[count - 1], &in[i], NULL);
1603 if (d > dmax)
1604 {
1605 idx = i;
1606 dmax = d;
1607 }
1608 }
1609 if (dmax > dist_sq)
1610 {
1611 ret = transform_douglas_peucker(in, idx, dist_sq, out) - 1;
1612 ret += transform_douglas_peucker(in + idx, count - idx, dist_sq, out + ret);
1613 }
1614 else
1615 {
1616 if (count > 0)
1617 out[ret++] = in[0];
1618 if (count > 1)
1619 out[ret++] = in[count - 1];
1620 }
1621 return ret;
1622 }
1623
1624 int transform_douglas_peucker_float(struct coord *in, int count, navit_float dist_sq, struct coord *out)
1625 {
1626 int ret = 0;
1627 int i, idx = 0;
1628 navit_float d, dmax = 0;
1629 for (i = 1; i < count - 2; i++)
1630 {
1631 d = transform_distance_line_sq_float(&in[0], &in[count - 1], &in[i], NULL);
1632 if (d > dmax)
1633 {
1634 idx = i;
1635 dmax = d;
1636 }
1637 }
1638 if (dmax > dist_sq)
1639 {
1640 ret = transform_douglas_peucker_float(in, idx, dist_sq, out) - 1;
1641 ret += transform_douglas_peucker_float(in + idx, count - idx, dist_sq, out + ret);
1642 }
1643 else
1644 {
1645 if (count > 0)
1646 out[ret++] = in[0];
1647 if (count > 1)
1648 out[ret++] = in[count - 1];
1649 }
1650 return ret;
1651 }
1652
1653 void transform_print_deg(double deg)
1654 {
1655 printf("%2.0f:%2.0f:%2.4f", floor(deg), fmod(deg * 60, 60), fmod(deg * 3600, 60));
1656 }
1657
1658 #ifdef AVOID_FLOAT
1659 static int tab_atan[]=
1660 { 0,262,524,787,1051,1317,1584,1853,2126,2401,2679,2962,3249,3541,3839,4142,4452,4770,5095,5430,5774,6128,6494,6873,7265,7673,8098,8541,9004,9490,10000,10538};
1661
1662 static int
1663 atan2_int_lookup(int val)
1664 {
1665 int len = sizeof(tab_atan)/sizeof(int);
1666 int i = len/2;
1667 int p = i-1;
1668
1669 for (;;)
1670 {
1671 i>>=1;
1672 if (val < tab_atan[p])
1673 {
1674 p-=i;
1675 }
1676 else
1677 {
1678 if (val < tab_atan[p+1])
1679 {
1680 return p+(p>>1);
1681 }
1682 else
1683 {
1684 p+=i;
1685 }
1686 }
1687 }
1688 }
1689
1690 static int
1691 atan2_int(int dx, int dy)
1692 {
1693 int mul=1,add=0,ret;
1694 if (! dx)
1695 {
1696 return dy < 0 ? 180 : 0;
1697 }
1698
1699 if (! dy)
1700 {
1701 return dx < 0 ? -90 : 90;
1702 }
1703
1704 if (dx < 0)
1705 {
1706 dx=-dx;
1707 mul=-1;
1708 }
1709
1710 if (dy < 0)
1711 {
1712 dy=-dy;
1713 add=180*mul;
1714 mul*=-1;
1715 }
1716
1717 while (dx > 20000 || dy > 20000)
1718 {
1719 dx/=10;
1720 dy/=10;
1721 }
1722
1723 if (dx > dy)
1724 {
1725 ret = 90 - atan2_int_lookup(dy*10000/dx);
1726 }
1727 else
1728 {
1729 ret = atan2_int_lookup(dx*10000/dy);
1730 }
1731
1732 return ret * mul + add;
1733 }
1734 #endif
1735
1736 int transform_get_angle_delta(struct coord *c1, struct coord *c2, int dir)
1737 {
1738 int dx = c2->x - c1->x;
1739 int dy = c2->y - c1->y;
1740
1741 #ifndef AVOID_FLOAT
1742
1743 //dbg(0, "**use FLOAT\n");
1744
1745 double angle;
1746 angle = atan2(dx, dy);
1747 angle *= 180 / M_PI;
1748 #else
1749
1750 //dbg(0, "++AVOID FLOAT\n");
1751
1752 int angle;
1753 angle = atan2_int(dx,dy);
1754 #endif
1755
1756 if (dir == -1)
1757 {
1758 angle = angle - 180;
1759 }
1760
1761 if (angle < 0)
1762 {
1763 angle += 360;
1764 }
1765
1766 return angle;
1767 }
1768
1769 int transform_get_angle_delta_accurate(struct coord *c1, struct coord *c2, int dir)
1770 {
1771 int dx = c2->x - c1->x;
1772 int dy = c2->y - c1->y;
1773
1774
1775 double angle;
1776 angle = atan2(dx, dy);
1777 angle *= 180 / M_PI;
1778
1779 if (dir == -1)
1780 {
1781 angle = angle - 180;
1782 }
1783
1784 if (angle < 0)
1785 {
1786 angle += 360;
1787 }
1788
1789 return angle;
1790 }
1791
1792
1793 int transform_within_border(struct transformation *this_, struct point *p, int border)
1794 {
1795 struct map_selection *ms = this_->screen_sel;
1796 while (ms)
1797 {
1798 struct point_rect *r = &ms->u.p_rect;
1799 if (p->x >= r->lu.x + border && p->x <= r->rl.x - border && p->y >= r->lu.y + border && p->y <= r->rl.y - border)
1800 return 1;
1801 ms = ms->next;
1802 }
1803 return 0;
1804 }
1805
1806 int transform_within_dist_point(struct coord *ref, struct coord *c, int dist)
1807 {
1808 if (c->x - dist > ref->x)
1809 return 0;
1810
1811 if (c->x + dist < ref->x)
1812 return 0;
1813
1814 if (c->y - dist > ref->y)
1815 return 0;
1816
1817 if (c->y + dist < ref->y)
1818 return 0;
1819
1820 if ((c->x - ref->x) * (c->x - ref->x) + (c->y - ref->y) * (c->y - ref->y) <= dist * dist)
1821 return 1;
1822
1823 return 0;
1824 }
1825
1826 int transform_within_dist_line(struct coord *ref, struct coord *c0, struct coord *c1, int dist)
1827 {
1828 int vx, vy, wx, wy;
1829 int n1, n2;
1830 struct coord lc;
1831
1832 if (c0->x < c1->x)
1833 {
1834 if (c0->x - dist > ref->x)
1835 return 0;
1836
1837 if (c1->x + dist < ref->x)
1838 return 0;
1839 }
1840 else
1841 {
1842 if (c1->x - dist > ref->x)
1843 return 0;
1844 if (c0->x + dist < ref->x)
1845 return 0;
1846 }
1847 if (c0->y < c1->y)
1848 {
1849 if (c0->y - dist > ref->y)
1850 return 0;
1851
1852 if (c1->y + dist < ref->y)
1853 return 0;
1854 }
1855 else
1856 {
1857 if (c1->y - dist > ref->y)
1858 return 0;
1859
1860 if (c0->y + dist < ref->y)
1861 return 0;
1862 }
1863 vx = c1->x - c0->x;
1864 vy = c1->y - c0->y;
1865 wx = ref->x - c0->x;
1866 wy = ref->y - c0->y;
1867
1868 n1 = vx * wx + vy * wy;
1869 if (n1 <= 0)
1870 return transform_within_dist_point(ref, c0, dist);
1871
1872 n2 = vx * vx + vy * vy;
1873
1874 if (n2 <= n1)
1875 return transform_within_dist_point(ref, c1, dist);
1876
1877 lc.x = c0->x + vx * n1 / n2;
1878 lc.y = c0->y + vy * n1 / n2;
1879
1880 return transform_within_dist_point(ref, &lc, dist);
1881 }
1882
1883 int transform_within_dist_polyline(struct coord *ref, struct coord *c, int count, int close, int dist)
1884 {
1885 int i;
1886 for (i = 0; i < count - 1; i++)
1887 {
1888 if (transform_within_dist_line(ref, c + i, c + i + 1, dist))
1889 {
1890 return 1;
1891 }
1892 }
1893
1894 if (close)
1895 return (transform_within_dist_line(ref, c, c + count - 1, dist));
1896
1897 return 0;
1898 }
1899
1900 int transform_within_dist_polygon(struct coord *ref, struct coord *c, int count, int dist)
1901 {
1902 int i, j, ci = 0;
1903 for (i = 0, j = count - 1; i < count; j = i++)
1904 {
1905 if ((((c[i].y <= ref->y) && (ref->y < c[j].y)) || ((c[j].y <= ref->y) && (ref->y < c[i].y))) && (ref->x < (c[j].x - c[i].x) * (ref->y - c[i].y) / (c[j].y - c[i].y) + c[i].x))
1906 ci = !ci;
1907 }
1908
1909 if (!ci)
1910 {
1911 if (dist)
1912 return transform_within_dist_polyline(ref, c, count, dist, 1);
1913 else
1914 return 0;
1915 }
1916
1917 return 1;
1918 }
1919
1920 int transform_within_dist_item(struct coord *ref, enum item_type type, struct coord *c, int count, int dist)
1921 {
1922 if (type < type_line)
1923 return transform_within_dist_point(ref, c, dist);
1924
1925 if (type < type_area)
1926 return transform_within_dist_polyline(ref, c, count, 0, dist);
1927
1928 return transform_within_dist_polygon(ref, c, count, dist);
1929 }
1930
1931 void transform_copy(struct transformation *src, struct transformation *dst)
1932 {
1933 memcpy(dst, src, sizeof(*src));
1934 }
1935
1936 void transform_destroy(struct transformation *t)
1937 {
1938 g_free(t);
1939 }
1940
1941 /*
1942 Note: there are many mathematically equivalent ways to express these formulas. As usual, not all of them are computationally equivalent.
1943
1944 L = latitude in radians (positive north)
1945 Lo = longitude in radians (positive east)
1946 E = easting (meters)
1947 N = northing (meters)
1948
1949 For the sphere
1950
1951 E = r Lo
1952 N = r ln [ tan (pi/4 + L/2) ]
1953
1954 where
1955
1956 r = radius of the sphere (meters)
1957 ln() is the natural logarithm
1958
1959 For the ellipsoid
1960
1961 E = a Lo
1962 N = a * ln ( tan (pi/4 + L/2) * ( (1 - e * sin (L)) / (1 + e * sin (L))) ** (e/2) )
1963
1964
1965 e
1966 -
1967 pi L 1 - e sin(L) 2
1968 = a ln( tan( ---- + ---) (--------------) )
1969 4 2 1 + e sin(L)
1970
1971
1972 where
1973
1974 a = the length of the semi-major axis of the ellipsoid (meters)
1975 e = the first eccentricity of the ellipsoid
1976
1977
1978 */
1979
1980 void transform_init(void)
1981 {
1982 #if 0
1983 if (global_transform_hash == NULL)
1984 {
1985 dbg(0,"enter\n");
1986 global_transform_hash=g_hash_table_new_full(g_int_hash, g_int_equal, NULL, g_free_func);
1987 dbg(0,"leave\n");
1988 }
1989 if (global_transform_hash2 == NULL)
1990 {
1991 dbg(0,"enter\n");
1992 global_transform_hash2=g_hash_table_new_full(g_int_hash, g_int_equal, NULL, g_free_func);
1993 dbg(0,"leave\n");
1994 }
1995 #endif
1996 }
1997

   
Visit the ZANavi Wiki