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

Contents of /navit/navit/transform.c

Parent Directory Parent Directory | Revision Log Revision Log


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

   
Visit the ZANavi Wiki