1 /*
2 * scan.c
3 * kirk johnson
4 * august 1993
5 *
6 * Copyright (C) 1989, 1990, 1993-1995, 1999 Kirk Lauritz Johnson
7 *
8 * Parts of the source code (as marked) are:
9 * Copyright (C) 1989, 1990, 1991 by Jim Frost
10 * Copyright (C) 1992 by Jamie Zawinski <jwz@lucid.com>
11 *
12 * Permission to use, copy, modify and freely distribute xearth for
13 * non-commercial and not-for-profit purposes is hereby granted
14 * without fee, provided that both the above copyright notice and this
15 * permission notice appear in all copies and in supporting
16 * documentation.
17 *
18 * Unisys Corporation holds worldwide patent rights on the Lempel Zev
19 * Welch (LZW) compression technique employed in the CompuServe GIF
20 * image file format as well as in other formats. Unisys has made it
21 * clear, however, that it does not require licensing or fees to be
22 * paid for freely distributed, non-commercial applications (such as
23 * xearth) that employ LZW/GIF technology. Those wishing further
24 * information about licensing the LZW patent should contact Unisys
25 * directly at (lzw_info@unisys.com) or by writing to
26 *
27 * Unisys Corporation
28 * Welch Licensing Department
29 * M/S-C1SW19
30 * P.O. Box 500
31 * Blue Bell, PA 19424
32 *
33 * The author makes no representations about the suitability of this
34 * software for any purpose. It is provided "as is" without express or
35 * implied warranty.
36 *
37 * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
38 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
39 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT
40 * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
41 * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
42 * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
43 * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
44 */
45
46 #include "xearth.h"
47 #include "kljcpyrt.h"
48
49 #define MAP_DATA_SCALE (30000)
50
51 #define XingTypeEntry (0)
52 #define XingTypeExit (1)
53
54 typedef struct
55 {
56 int type;
57 int cidx;
58 double x, y;
59 double angle;
60 } EdgeXing;
61
62 void scan_map _P((void));
63 void orth_scan_outline _P((void));
64 void orth_scan_curves _P((void));
65 double *orth_extract_curve _P((int, short *));
66 void orth_scan_along_curve _P((double *, double *, int));
67 void orth_find_edge_xing _P((double *, double *, double *));
68 void orth_handle_xings _P((void));
69 void orth_scan_arc _P((double, double, double, double, double, double));
70 void merc_scan_outline _P((void));
71 void merc_scan_curves _P((void));
72 double *merc_extract_curve _P((int, short *));
73 void merc_scan_along_curve _P((double *, double *, int));
74 double merc_find_edge_xing _P((double *, double *));
75 void merc_handle_xings _P((void));
76 void merc_scan_edge _P((EdgeXing *, EdgeXing *));
77 void cyl_scan_outline _P((void));
78 void cyl_scan_curves _P((void));
79 double *cyl_extract_curve _P((int, short *));
80 void cyl_scan_along_curve _P((double *, double *, int));
81 double cyl_find_edge_xing _P((double *, double *));
82 void cyl_handle_xings _P((void));
83 void cyl_scan_edge _P((EdgeXing *, EdgeXing *));
84 void xing_error _P((const char *, int, int, int, EdgeXing *));
85 void scan _P((double, double, double, double));
86 void get_scanbits _P((int));
87
88 static int double_comp _P((const void *, const void *));
89 static int scanbit_comp _P((const void *, const void *));
90 static int orth_edgexing_comp _P((const void *, const void *));
91 static int merc_edgexing_comp _P((const void *, const void *));
92 static int cyl_edgexing_comp _P((const void *, const void *));
93
94 ViewPosInfo view_pos_info;
95 ProjInfo proj_info;
96
97 int first_scan = 1;
98 ExtArr scanbits;
99 ExtArr edgexings;
100 int min_y, max_y;
101 ExtArr *scanbuf;
102
103
double_comp(a,b)104 static int double_comp(a, b)
105 const void *a;
106 const void *b;
107 {
108 double val_a;
109 double val_b;
110 int rslt;
111
112 val_a = *((const double *) a);
113 val_b = *((const double *) b);
114
115 if (val_a < val_b)
116 rslt = -1;
117 else if (val_a > val_b)
118 rslt = 1;
119 else
120 rslt = 0;
121
122 return rslt;
123 }
124
125
scanbit_comp(a,b)126 static int scanbit_comp(a, b)
127 const void *a;
128 const void *b;
129 {
130 return (((const ScanBit *) a)->y - ((const ScanBit *) b)->y);
131 }
132
133
orth_edgexing_comp(a,b)134 static int orth_edgexing_comp(a, b)
135 const void *a;
136 const void *b;
137 {
138 double angle_a;
139 double angle_b;
140 int rslt;
141
142 angle_a = ((const EdgeXing *) a)->angle;
143 angle_b = ((const EdgeXing *) b)->angle;
144
145 if (angle_a < angle_b)
146 rslt = -1;
147 else if (angle_a > angle_b)
148 rslt = 1;
149 else
150 rslt = 0;
151
152 return rslt;
153 }
154
155
merc_edgexing_comp(a,b)156 static int merc_edgexing_comp(a, b)
157 const void *a;
158 const void *b;
159 {
160 double val_a;
161 double val_b;
162 int rslt;
163
164 val_a = ((const EdgeXing *) a)->angle;
165 val_b = ((const EdgeXing *) b)->angle;
166
167 if (val_a < val_b)
168 {
169 rslt = -1;
170 }
171 else if (val_a > val_b)
172 {
173 rslt = 1;
174 }
175 else if (val_a == 0)
176 {
177 val_a = ((const EdgeXing *) a)->y;
178 val_b = ((const EdgeXing *) b)->y;
179
180 if (val_a < val_b)
181 rslt = -1;
182 else if (val_a > val_b)
183 rslt = 1;
184 else
185 rslt = 0;
186 }
187 else if (val_a == 2)
188 {
189 val_a = ((const EdgeXing *) a)->y;
190 val_b = ((const EdgeXing *) b)->y;
191
192 if (val_a > val_b)
193 rslt = -1;
194 else if (val_a < val_b)
195 rslt = 1;
196 else
197 rslt = 0;
198 }
199 else
200 {
201 /* keep lint happy */
202 rslt = 0;
203 assert(0);
204 }
205
206 return rslt;
207 }
208
209
cyl_edgexing_comp(a,b)210 static int cyl_edgexing_comp(a, b)
211 const void *a;
212 const void *b;
213 {
214 double val_a;
215 double val_b;
216 int rslt;
217
218 val_a = ((const EdgeXing *) a)->angle;
219 val_b = ((const EdgeXing *) b)->angle;
220
221 if (val_a < val_b)
222 {
223 rslt = -1;
224 }
225 else if (val_a > val_b)
226 {
227 rslt = 1;
228 }
229 else if (val_a == 0)
230 {
231 val_a = ((const EdgeXing *) a)->y;
232 val_b = ((const EdgeXing *) b)->y;
233
234 if (val_a < val_b)
235 rslt = -1;
236 else if (val_a > val_b)
237 rslt = 1;
238 else
239 rslt = 0;
240 }
241 else if (val_a == 2)
242 {
243 val_a = ((const EdgeXing *) a)->y;
244 val_b = ((const EdgeXing *) b)->y;
245
246 if (val_a > val_b)
247 rslt = -1;
248 else if (val_a < val_b)
249 rslt = 1;
250 else
251 rslt = 0;
252 }
253 else
254 {
255 /* keep lint happy */
256 rslt = 0;
257 assert(0);
258 }
259
260 return rslt;
261 }
262
263
scan_map()264 void scan_map()
265 {
266 int i;
267 ViewPosInfo *vpi;
268 ProjInfo *pi;
269
270 vpi = &view_pos_info;
271 vpi->cos_lat = cos(view_lat * (M_PI/180));
272 vpi->sin_lat = sin(view_lat * (M_PI/180));
273 vpi->cos_lon = cos(view_lon * (M_PI/180));
274 vpi->sin_lon = sin(view_lon * (M_PI/180));
275 vpi->cos_rot = cos(view_rot * (M_PI/180));
276 vpi->sin_rot = sin(view_rot * (M_PI/180));
277
278 pi = &proj_info;
279 if (proj_type == ProjTypeOrthographic)
280 {
281 pi->proj_scale = ((hght < wdth) ? hght : wdth) * (view_mag / 2) * 0.99;
282 }
283 else
284 {
285 /* proj_type is either ProjTypeMercator or ProjTypeCylindrical
286 */
287 pi->proj_scale = (view_mag * wdth) / (2 * M_PI);
288 }
289
290 pi->proj_xofs = (double) wdth / 2 + shift_x;
291 pi->proj_yofs = (double) hght / 2 + shift_y;
292 pi->inv_proj_scale = 1 / pi->proj_scale;
293
294 /* the first time through, allocate scanbits and edgexings;
295 * on subsequent passes, simply reset them.
296 */
297 if (first_scan)
298 {
299 scanbits = extarr_alloc(sizeof(ScanBit));
300 edgexings = extarr_alloc(sizeof(EdgeXing));
301 }
302 else
303 {
304 scanbits->count = 0;
305 edgexings->count = 0;
306 }
307
308 /* maybe only allocate these once and reset them on
309 * subsequent passes (like scanbits and edgexings)?
310 */
311 scanbuf = (ExtArr *) malloc((unsigned) sizeof(ExtArr) * hght);
312 assert(scanbuf != NULL);
313 for (i=0; i<hght; i++)
314 scanbuf[i] = extarr_alloc(sizeof(double));
315
316 if (proj_type == ProjTypeOrthographic)
317 {
318 orth_scan_outline();
319 orth_scan_curves();
320 }
321 else if (proj_type == ProjTypeMercator)
322 {
323 merc_scan_outline();
324 merc_scan_curves();
325 }
326 else /* (proj_type == ProjTypeCylindrical) */
327 {
328 cyl_scan_outline();
329 cyl_scan_curves();
330 }
331
332 for (i=0; i<hght; i++)
333 extarr_free(scanbuf[i]);
334 free(scanbuf);
335
336 qsort(scanbits->body, scanbits->count, sizeof(ScanBit), scanbit_comp);
337
338 first_scan = 0;
339 }
340
341
orth_scan_outline()342 void orth_scan_outline()
343 {
344 min_y = hght;
345 max_y = -1;
346
347 orth_scan_arc(1.0, 0.0, 0.0, 1.0, 0.0, (2*M_PI));
348
349 get_scanbits(64);
350 }
351
352
orth_scan_curves()353 void orth_scan_curves()
354 {
355 int i;
356 int cidx;
357 int npts;
358 int val;
359 short *raw;
360 double *pos;
361 double *prev;
362 double *curr;
363
364 cidx = 0;
365 raw = map_data;
366 while (1)
367 {
368 npts = raw[0];
369 if (npts == 0) break;
370 val = raw[1];
371 raw += 2;
372
373 pos = orth_extract_curve(npts, raw);
374 prev = pos + (npts-1)*3;
375 curr = pos;
376 min_y = hght;
377 max_y = -1;
378
379 for (i=0; i<npts; i++)
380 {
381 orth_scan_along_curve(prev, curr, cidx);
382 prev = curr;
383 curr += 3;
384 }
385
386 free(pos);
387 if (edgexings->count > 0)
388 orth_handle_xings();
389 if (min_y <= max_y)
390 get_scanbits(val);
391
392 cidx += 1;
393 raw += 3*npts;
394 }
395 }
396
397
orth_extract_curve(npts,data)398 double *orth_extract_curve(npts, data)
399 int npts;
400 short *data;
401 {
402 int i;
403 int x, y, z;
404 double scale;
405 double *pos;
406 double *rslt;
407
408 rslt = (double *) malloc((unsigned) sizeof(double) * 3 * npts);
409 assert(rslt != NULL);
410
411 x = 0;
412 y = 0;
413 z = 0;
414 scale = 1.0 / MAP_DATA_SCALE;
415 pos = rslt;
416
417 for (i=0; i<npts; i++)
418 {
419 x += data[0];
420 y += data[1];
421 z += data[2];
422
423 pos[0] = x * scale;
424 pos[1] = y * scale;
425 pos[2] = z * scale;
426
427 XFORM_ROTATE(pos, view_pos_info);
428
429 data += 3;
430 pos += 3;
431 }
432
433 return rslt;
434 }
435
436
orth_scan_along_curve(prev,curr,cidx)437 void orth_scan_along_curve(prev, curr, cidx)
438 double *prev;
439 double *curr;
440 int cidx;
441 {
442 double extra[3];
443 EdgeXing *xing;
444
445 if (prev[2] <= 0) /* prev not visible */
446 {
447 if (curr[2] <= 0)
448 return; /* neither point visible */
449
450 orth_find_edge_xing(prev, curr, extra);
451
452 /* extra[] is an edge crossing (entry point) */
453 xing = (EdgeXing *) extarr_next(edgexings);
454 xing->type = XingTypeEntry;
455 xing->cidx = cidx;
456 xing->x = extra[0];
457 xing->y = extra[1];
458 xing->angle = atan2(extra[1], extra[0]);
459
460 prev = extra;
461 }
462 else if (curr[2] <= 0) /* curr not visible */
463 {
464 orth_find_edge_xing(prev, curr, extra);
465
466 /* extra[] is an edge crossing (exit point) */
467 xing = (EdgeXing *) extarr_next(edgexings);
468 xing->type = XingTypeExit;
469 xing->cidx = cidx;
470 xing->x = extra[0];
471 xing->y = extra[1];
472 xing->angle = atan2(extra[1], extra[0]);
473
474 curr = extra;
475 }
476
477 scan(XPROJECT(prev[0]), YPROJECT(prev[1]),
478 XPROJECT(curr[0]), YPROJECT(curr[1]));
479 }
480
481
orth_find_edge_xing(prev,curr,rslt)482 void orth_find_edge_xing(prev, curr, rslt)
483 double *prev;
484 double *curr;
485 double *rslt;
486 {
487 double tmp;
488 double r0, r1;
489
490 tmp = curr[2] / (curr[2] - prev[2]);
491 r0 = curr[0] - tmp * (curr[0] - prev[0]);
492 r1 = curr[1] - tmp * (curr[1] - prev[1]);
493
494 tmp = sqrt((r0*r0) + (r1*r1));
495 rslt[0] = r0 / tmp;
496 rslt[1] = r1 / tmp;
497 rslt[2] = 0;
498 }
499
500
orth_handle_xings()501 void orth_handle_xings()
502 {
503 int i;
504 int nxings;
505 EdgeXing *xings;
506 EdgeXing *from;
507 EdgeXing *to;
508
509 xings = (EdgeXing *) edgexings->body;
510 nxings = edgexings->count;
511
512 assert((nxings % 2) == 0);
513 qsort(xings, (unsigned) nxings, sizeof(EdgeXing), orth_edgexing_comp);
514
515 if (xings[0].type == XingTypeExit)
516 {
517 for (i=0; i<nxings; i+=2)
518 {
519 from = &(xings[i]);
520 to = &(xings[i+1]);
521
522 if ((from->type != XingTypeExit) ||
523 (to->type != XingTypeEntry))
524 xing_error(__FILE__, __LINE__, i, nxings, xings);
525
526 orth_scan_arc(from->x, from->y, from->angle,
527 to->x, to->y, to->angle);
528 }
529 }
530 else
531 {
532 from = &(xings[nxings-1]);
533 to = &(xings[0]);
534
535 if ((from->type != XingTypeExit) ||
536 (to->type != XingTypeEntry) ||
537 (from->angle < to->angle))
538 xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
539
540 orth_scan_arc(from->x, from->y, from->angle,
541 to->x, to->y, to->angle+(2*M_PI));
542
543 for (i=1; i<(nxings-1); i+=2)
544 {
545 from = &(xings[i]);
546 to = &(xings[i+1]);
547
548 if ((from->type != XingTypeExit) ||
549 (to->type != XingTypeEntry))
550 xing_error(__FILE__, __LINE__, i, nxings, xings);
551
552 orth_scan_arc(from->x, from->y, from->angle,
553 to->x, to->y, to->angle);
554 }
555 }
556
557 edgexings->count = 0;
558 }
559
560
orth_scan_arc(x_0,y_0,a_0,x_1,y_1,a_1)561 void orth_scan_arc(x_0, y_0, a_0, x_1, y_1, a_1)
562 double x_0, y_0, a_0;
563 double x_1, y_1, a_1;
564 {
565 int i;
566 int lo, hi;
567 double angle, step;
568 double prev_x, prev_y;
569 double curr_x, curr_y;
570 double c_step, s_step;
571 double arc_x, arc_y;
572 double tmp;
573
574 assert(a_0 < a_1);
575
576 step = proj_info.inv_proj_scale * 10;
577 if (step > 0.05) step = 0.05;
578 lo = ceil(a_0 / step);
579 hi = floor(a_1 / step);
580
581 prev_x = XPROJECT(x_0);
582 prev_y = YPROJECT(y_0);
583
584 if (lo <= hi)
585 {
586 c_step = cos(step);
587 s_step = sin(step);
588
589 angle = lo * step;
590 arc_x = cos(angle);
591 arc_y = sin(angle);
592
593 for (i=lo; i<=hi; i++)
594 {
595 curr_x = XPROJECT(arc_x);
596 curr_y = YPROJECT(arc_y);
597 scan(prev_x, prev_y, curr_x, curr_y);
598
599 /* instead of repeatedly calling cos() and sin() to get the next
600 * values for arc_x and arc_y, simply rotate the existing values
601 */
602 tmp = (c_step * arc_x) - (s_step * arc_y);
603 arc_y = (s_step * arc_x) + (c_step * arc_y);
604 arc_x = tmp;
605
606 prev_x = curr_x;
607 prev_y = curr_y;
608 }
609 }
610
611 curr_x = XPROJECT(x_1);
612 curr_y = YPROJECT(y_1);
613 scan(prev_x, prev_y, curr_x, curr_y);
614 }
615
616
merc_scan_outline()617 void merc_scan_outline()
618 {
619 double left, right;
620 double top, bottom;
621
622 min_y = hght;
623 max_y = -1;
624
625 left = XPROJECT(-M_PI);
626 right = XPROJECT(M_PI);
627 top = YPROJECT(BigNumber);
628 bottom = YPROJECT(-BigNumber);
629
630 scan(right, top, left, top);
631 scan(left, top, left, bottom);
632 scan(left, bottom, right, bottom);
633 scan(right, bottom, right, top);
634
635 get_scanbits(64);
636 }
637
638
merc_scan_curves()639 void merc_scan_curves()
640 {
641 int i;
642 int cidx;
643 int npts;
644 int val;
645 short *raw;
646 double *pos;
647 double *prev;
648 double *curr;
649
650 cidx = 0;
651 raw = map_data;
652 while (1)
653 {
654 npts = raw[0];
655 if (npts == 0) break;
656 val = raw[1];
657 raw += 2;
658
659 pos = merc_extract_curve(npts, raw);
660 prev = pos + (npts-1)*5;
661 curr = pos;
662 min_y = hght;
663 max_y = -1;
664
665 for (i=0; i<npts; i++)
666 {
667 merc_scan_along_curve(prev, curr, cidx);
668 prev = curr;
669 curr += 5;
670 }
671
672 free(pos);
673 if (edgexings->count > 0)
674 merc_handle_xings();
675 if (min_y <= max_y)
676 get_scanbits(val);
677
678 cidx += 1;
679 raw += 3*npts;
680 }
681 }
682
683
merc_extract_curve(npts,data)684 double *merc_extract_curve(npts, data)
685 int npts;
686 short *data;
687 {
688 int i;
689 int x, y, z;
690 double scale;
691 double *pos;
692 double *rslt;
693
694 rslt = (double *) malloc((unsigned) sizeof(double) * 5 * npts);
695 assert(rslt != NULL);
696
697 x = 0;
698 y = 0;
699 z = 0;
700 scale = 1.0 / MAP_DATA_SCALE;
701 pos = rslt;
702
703 for (i=0; i<npts; i++)
704 {
705 x += data[0];
706 y += data[1];
707 z += data[2];
708
709 pos[0] = x * scale;
710 pos[1] = y * scale;
711 pos[2] = z * scale;
712
713 XFORM_ROTATE(pos, view_pos_info);
714
715 /* apply mercator projection
716 */
717 pos[3] = MERCATOR_X(pos[0], pos[2]);
718 pos[4] = MERCATOR_Y(pos[1]);
719
720 data += 3;
721 pos += 5;
722 }
723
724 return rslt;
725 }
726
727
merc_scan_along_curve(prev,curr,cidx)728 void merc_scan_along_curve(prev, curr, cidx)
729 double *prev;
730 double *curr;
731 int cidx;
732 {
733 double px, py;
734 double cx, cy;
735 double dx;
736 double mx, my;
737 EdgeXing *xing;
738
739 px = prev[3];
740 cx = curr[3];
741 py = prev[4];
742 cy = curr[4];
743 dx = cx - px;
744
745 if (dx > 0)
746 {
747 /* curr to the right of prev
748 */
749
750 if (dx > ((2*M_PI) - dx))
751 {
752 /* vertical edge crossing to the left of prev
753 */
754
755 /* find exit point (left edge) */
756 mx = - M_PI;
757 my = merc_find_edge_xing(prev, curr);
758
759 /* scan from prev to exit point */
760 scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
761
762 /* (mx, my) is an edge crossing (exit point) */
763 xing = (EdgeXing *) extarr_next(edgexings);
764 xing->type = XingTypeExit;
765 xing->cidx = cidx;
766 xing->x = mx;
767 xing->y = my;
768 xing->angle = 2; /* left edge */
769
770 /* scan from entry point (right edge) to curr */
771 mx = M_PI;
772 scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
773
774 /* (mx, my) is an edge crossing (entry point) */
775 xing = (EdgeXing *) extarr_next(edgexings);
776 xing->type = XingTypeEntry;
777 xing->cidx = cidx;
778 xing->x = mx;
779 xing->y = my;
780 xing->angle = 0; /* right edge */
781 }
782 else
783 {
784 /* no vertical edge crossing
785 */
786 scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
787 }
788 }
789 else
790 {
791 /* curr to the left of prev
792 */
793 dx = - dx;
794
795 if (dx > ((2*M_PI) - dx))
796 {
797 /* vertical edge crossing to the right of prev
798 */
799
800 /* find exit point (right edge) */
801 mx = M_PI;
802 my = merc_find_edge_xing(prev, curr);
803
804 /* scan from prev to exit point */
805 scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
806
807 /* (mx, my) is an edge crossing (exit point) */
808 xing = (EdgeXing *) extarr_next(edgexings);
809 xing->type = XingTypeExit;
810 xing->cidx = cidx;
811 xing->x = mx;
812 xing->y = my;
813 xing->angle = 0; /* right edge */
814
815 /* scan from entry point (left edge) to curr */
816 mx = - M_PI;
817 scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
818
819 /* (mx, my) is an edge crossing (entry point) */
820 xing = (EdgeXing *) extarr_next(edgexings);
821 xing->type = XingTypeEntry;
822 xing->cidx = cidx;
823 xing->x = mx;
824 xing->y = my;
825 xing->angle = 2; /* left edge */
826 }
827 else
828 {
829 /* no vertical edge crossing
830 */
831 scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
832 }
833 }
834 }
835
836
merc_find_edge_xing(prev,curr)837 double merc_find_edge_xing(prev, curr)
838 double *prev;
839 double *curr;
840 {
841 double ratio;
842 double scale;
843 double z1, z2;
844 double rslt;
845
846 if (curr[0] != 0)
847 {
848 ratio = (prev[0] / curr[0]);
849 z1 = prev[1] - (ratio * curr[1]);
850 z2 = prev[2] - (ratio * curr[2]);
851 }
852 else
853 {
854 z1 = curr[1];
855 z2 = curr[2];
856 }
857
858 scale = ((z2 > 0) ? -1 : 1) / sqrt((z1*z1) + (z2*z2));
859 z1 *= scale;
860
861 rslt = MERCATOR_Y(z1);
862
863 return rslt;
864 }
865
866
merc_handle_xings()867 void merc_handle_xings()
868 {
869 int i;
870 int nxings;
871 EdgeXing *xings;
872 EdgeXing *from;
873 EdgeXing *to;
874
875 xings = (EdgeXing *) edgexings->body;
876 nxings = edgexings->count;
877
878 assert((nxings % 2) == 0);
879 qsort(xings, (unsigned) nxings, sizeof(EdgeXing), merc_edgexing_comp);
880
881 if (xings[0].type == XingTypeExit)
882 {
883 for (i=0; i<nxings; i+=2)
884 {
885 from = &(xings[i]);
886 to = &(xings[i+1]);
887
888 if ((from->type != XingTypeExit) ||
889 (to->type != XingTypeEntry))
890 xing_error(__FILE__, __LINE__, i, nxings, xings);
891
892 merc_scan_edge(from, to);
893 }
894 }
895 else
896 {
897 from = &(xings[nxings-1]);
898 to = &(xings[0]);
899
900 if ((from->type != XingTypeExit) ||
901 (to->type != XingTypeEntry) ||
902 (from->angle < to->angle))
903 xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
904
905 merc_scan_edge(from, to);
906
907 for (i=1; i<(nxings-1); i+=2)
908 {
909 from = &(xings[i]);
910 to = &(xings[i+1]);
911
912 if ((from->type != XingTypeExit) ||
913 (to->type != XingTypeEntry))
914 xing_error(__FILE__, __LINE__, i, nxings, xings);
915
916 merc_scan_edge(from, to);
917 }
918 }
919
920 edgexings->count = 0;
921 }
922
923
merc_scan_edge(from,to)924 void merc_scan_edge(from, to)
925 EdgeXing *from;
926 EdgeXing *to;
927 {
928 int s0, s1, s_new;
929 double x_0, x_1, x_new;
930 double y_0, y_1, y_new;
931
932 s0 = from->angle;
933 x_0 = XPROJECT(from->x);
934 y_0 = YPROJECT(from->y);
935
936 s1 = to->angle;
937 x_1 = XPROJECT(to->x);
938 y_1 = YPROJECT(to->y);
939
940 while (s0 != s1)
941 {
942 switch (s0)
943 {
944 case 0:
945 x_new = XPROJECT(M_PI);
946 y_new = YPROJECT(BigNumber);
947 s_new = 1;
948 break;
949
950 case 1:
951 x_new = XPROJECT(-M_PI);
952 y_new = YPROJECT(BigNumber);
953 s_new = 2;
954 break;
955
956 case 2:
957 x_new = XPROJECT(-M_PI);
958 y_new = YPROJECT(-BigNumber);
959 s_new = 3;
960 break;
961
962 case 3:
963 x_new = XPROJECT(M_PI);
964 y_new = YPROJECT(-BigNumber);
965 s_new = 0;
966 break;
967
968 default:
969 /* keep lint happy */
970 x_new = 0;
971 y_new = 0;
972 s_new = 0;
973 assert(0);
974 }
975
976 scan(x_0, y_0, x_new, y_new);
977 x_0 = x_new;
978 y_0 = y_new;
979 s0 = s_new;
980 }
981
982 scan(x_0, y_0, x_1, y_1);
983 }
984
985
cyl_scan_outline()986 void cyl_scan_outline()
987 {
988 double left, right;
989 double top, bottom;
990
991 min_y = hght;
992 max_y = -1;
993
994 left = XPROJECT(-M_PI);
995 right = XPROJECT(M_PI);
996 top = YPROJECT(BigNumber);
997 bottom = YPROJECT(-BigNumber);
998
999 scan(right, top, left, top);
1000 scan(left, top, left, bottom);
1001 scan(left, bottom, right, bottom);
1002 scan(right, bottom, right, top);
1003
1004 get_scanbits(64);
1005 }
1006
1007
cyl_scan_curves()1008 void cyl_scan_curves()
1009 {
1010 int i;
1011 int cidx;
1012 int npts;
1013 int val;
1014 short *raw;
1015 double *pos;
1016 double *prev;
1017 double *curr;
1018
1019 cidx = 0;
1020 raw = map_data;
1021 while (1)
1022 {
1023 npts = raw[0];
1024 if (npts == 0) break;
1025 val = raw[1];
1026 raw += 2;
1027
1028 pos = cyl_extract_curve(npts, raw);
1029 prev = pos + (npts-1)*5;
1030 curr = pos;
1031 min_y = hght;
1032 max_y = -1;
1033
1034 for (i=0; i<npts; i++)
1035 {
1036 cyl_scan_along_curve(prev, curr, cidx);
1037 prev = curr;
1038 curr += 5;
1039 }
1040
1041 free(pos);
1042 if (edgexings->count > 0)
1043 cyl_handle_xings();
1044 if (min_y <= max_y)
1045 get_scanbits(val);
1046
1047 cidx += 1;
1048 raw += 3*npts;
1049 }
1050 }
1051
1052
cyl_extract_curve(npts,data)1053 double *cyl_extract_curve(npts, data)
1054 int npts;
1055 short *data;
1056 {
1057 int i;
1058 int x, y, z;
1059 double scale;
1060 double *pos;
1061 double *rslt;
1062
1063 rslt = (double *) malloc((unsigned) sizeof(double) * 5 * npts);
1064 assert(rslt != NULL);
1065
1066 x = 0;
1067 y = 0;
1068 z = 0;
1069 scale = 1.0 / MAP_DATA_SCALE;
1070 pos = rslt;
1071
1072 for (i=0; i<npts; i++)
1073 {
1074 x += data[0];
1075 y += data[1];
1076 z += data[2];
1077
1078 pos[0] = x * scale;
1079 pos[1] = y * scale;
1080 pos[2] = z * scale;
1081
1082 XFORM_ROTATE(pos, view_pos_info);
1083
1084 /* apply cylindrical projection
1085 */
1086 pos[3] = CYLINDRICAL_X(pos[0], pos[2]);
1087 pos[4] = CYLINDRICAL_Y(pos[1]);
1088
1089 data += 3;
1090 pos += 5;
1091 }
1092
1093 return rslt;
1094 }
1095
1096
cyl_scan_along_curve(prev,curr,cidx)1097 void cyl_scan_along_curve(prev, curr, cidx)
1098 double *prev;
1099 double *curr;
1100 int cidx;
1101 {
1102 double px, py;
1103 double cx, cy;
1104 double dx;
1105 double mx, my;
1106 EdgeXing *xing;
1107
1108 px = prev[3];
1109 cx = curr[3];
1110 py = prev[4];
1111 cy = curr[4];
1112 dx = cx - px;
1113
1114 if (dx > 0)
1115 {
1116 /* curr to the right of prev
1117 */
1118
1119 if (dx > ((2*M_PI) - dx))
1120 {
1121 /* vertical edge crossing to the left of prev
1122 */
1123
1124 /* find exit point (left edge) */
1125 mx = - M_PI;
1126 my = cyl_find_edge_xing(prev, curr);
1127
1128 /* scan from prev to exit point */
1129 scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
1130
1131 /* (mx, my) is an edge crossing (exit point) */
1132 xing = (EdgeXing *) extarr_next(edgexings);
1133 xing->type = XingTypeExit;
1134 xing->cidx = cidx;
1135 xing->x = mx;
1136 xing->y = my;
1137 xing->angle = 2; /* left edge */
1138
1139 /* scan from entry point (right edge) to curr */
1140 mx = M_PI;
1141 scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
1142
1143 /* (mx, my) is an edge crossing (entry point) */
1144 xing = (EdgeXing *) extarr_next(edgexings);
1145 xing->type = XingTypeEntry;
1146 xing->cidx = cidx;
1147 xing->x = mx;
1148 xing->y = my;
1149 xing->angle = 0; /* right edge */
1150 }
1151 else
1152 {
1153 /* no vertical edge crossing
1154 */
1155 scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
1156 }
1157 }
1158 else
1159 {
1160 /* curr to the left of prev
1161 */
1162 dx = - dx;
1163
1164 if (dx > ((2*M_PI) - dx))
1165 {
1166 /* vertical edge crossing to the right of prev
1167 */
1168
1169 /* find exit point (right edge) */
1170 mx = M_PI;
1171 my = cyl_find_edge_xing(prev, curr);
1172
1173 /* scan from prev to exit point */
1174 scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
1175
1176 /* (mx, my) is an edge crossing (exit point) */
1177 xing = (EdgeXing *) extarr_next(edgexings);
1178 xing->type = XingTypeExit;
1179 xing->cidx = cidx;
1180 xing->x = mx;
1181 xing->y = my;
1182 xing->angle = 0; /* right edge */
1183
1184 /* scan from entry point (left edge) to curr */
1185 mx = - M_PI;
1186 scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
1187
1188 /* (mx, my) is an edge crossing (entry point) */
1189 xing = (EdgeXing *) extarr_next(edgexings);
1190 xing->type = XingTypeEntry;
1191 xing->cidx = cidx;
1192 xing->x = mx;
1193 xing->y = my;
1194 xing->angle = 2; /* left edge */
1195 }
1196 else
1197 {
1198 /* no vertical edge crossing
1199 */
1200 scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
1201 }
1202 }
1203 }
1204
1205
cyl_find_edge_xing(prev,curr)1206 double cyl_find_edge_xing(prev, curr)
1207 double *prev;
1208 double *curr;
1209 {
1210 double ratio;
1211 double scale;
1212 double z1, z2;
1213 double rslt;
1214
1215 if (curr[0] != 0)
1216 {
1217 ratio = (prev[0] / curr[0]);
1218 z1 = prev[1] - (ratio * curr[1]);
1219 z2 = prev[2] - (ratio * curr[2]);
1220 }
1221 else
1222 {
1223 z1 = curr[1];
1224 z2 = curr[2];
1225 }
1226
1227 scale = ((z2 > 0) ? -1 : 1) / sqrt((z1*z1) + (z2*z2));
1228 z1 *= scale;
1229
1230 rslt = CYLINDRICAL_Y(z1);
1231
1232 return rslt;
1233 }
1234
1235
cyl_handle_xings()1236 void cyl_handle_xings()
1237 {
1238 int i;
1239 int nxings;
1240 EdgeXing *xings;
1241 EdgeXing *from;
1242 EdgeXing *to;
1243
1244 xings = (EdgeXing *) edgexings->body;
1245 nxings = edgexings->count;
1246
1247 assert((nxings % 2) == 0);
1248 qsort(xings, (unsigned) nxings, sizeof(EdgeXing), cyl_edgexing_comp);
1249
1250 if (xings[0].type == XingTypeExit)
1251 {
1252 for (i=0; i<nxings; i+=2)
1253 {
1254 from = &(xings[i]);
1255 to = &(xings[i+1]);
1256
1257 if ((from->type != XingTypeExit) ||
1258 (to->type != XingTypeEntry))
1259 xing_error(__FILE__, __LINE__, i, nxings, xings);
1260
1261 cyl_scan_edge(from, to);
1262 }
1263 }
1264 else
1265 {
1266 from = &(xings[nxings-1]);
1267 to = &(xings[0]);
1268
1269 if ((from->type != XingTypeExit) ||
1270 (to->type != XingTypeEntry) ||
1271 (from->angle < to->angle))
1272 xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
1273
1274 cyl_scan_edge(from, to);
1275
1276 for (i=1; i<(nxings-1); i+=2)
1277 {
1278 from = &(xings[i]);
1279 to = &(xings[i+1]);
1280
1281 if ((from->type != XingTypeExit) ||
1282 (to->type != XingTypeEntry))
1283 xing_error(__FILE__, __LINE__, i, nxings, xings);
1284
1285 cyl_scan_edge(from, to);
1286 }
1287 }
1288
1289 edgexings->count = 0;
1290 }
1291
1292
cyl_scan_edge(from,to)1293 void cyl_scan_edge(from, to)
1294 EdgeXing *from;
1295 EdgeXing *to;
1296 {
1297 int s0, s1, s_new;
1298 double x_0, x_1, x_new;
1299 double y_0, y_1, y_new;
1300
1301 s0 = from->angle;
1302 x_0 = XPROJECT(from->x);
1303 y_0 = YPROJECT(from->y);
1304
1305 s1 = to->angle;
1306 x_1 = XPROJECT(to->x);
1307 y_1 = YPROJECT(to->y);
1308
1309 while (s0 != s1)
1310 {
1311 switch (s0)
1312 {
1313 case 0:
1314 x_new = XPROJECT(M_PI);
1315 y_new = YPROJECT(BigNumber);
1316 s_new = 1;
1317 break;
1318
1319 case 1:
1320 x_new = XPROJECT(-M_PI);
1321 y_new = YPROJECT(BigNumber);
1322 s_new = 2;
1323 break;
1324
1325 case 2:
1326 x_new = XPROJECT(-M_PI);
1327 y_new = YPROJECT(-BigNumber);
1328 s_new = 3;
1329 break;
1330
1331 case 3:
1332 x_new = XPROJECT(M_PI);
1333 y_new = YPROJECT(-BigNumber);
1334 s_new = 0;
1335 break;
1336
1337 default:
1338 /* keep lint happy */
1339 x_new = 0;
1340 y_new = 0;
1341 s_new = 0;
1342 assert(0);
1343 }
1344
1345 scan(x_0, y_0, x_new, y_new);
1346 x_0 = x_new;
1347 y_0 = y_new;
1348 s0 = s_new;
1349 }
1350
1351 scan(x_0, y_0, x_1, y_1);
1352 }
1353
1354
xing_error(file,line,idx,nxings,xings)1355 void xing_error(file, line, idx, nxings, xings)
1356 const char *file;
1357 int line;
1358 int idx;
1359 int nxings;
1360 EdgeXing *xings;
1361 {
1362 fflush(stdout);
1363 fprintf(stderr, "xearth %s: incorrect edgexing sequence (%s:%d)\n",
1364 VersionString, file, line);
1365 fprintf(stderr, " (cidx %d) (view_lat %.16f) (view_lon %.16f)\n",
1366 xings[idx].cidx, view_lat, view_lon);
1367 fprintf(stderr, "\n");
1368 exit(1);
1369 }
1370
1371
scan(x_0,y_0,x_1,y_1)1372 void scan(x_0, y_0, x_1, y_1)
1373 double x_0, y_0;
1374 double x_1, y_1;
1375 {
1376 int i;
1377 int lo_y, hi_y;
1378 double x_value;
1379 double x_delta;
1380
1381 if (y_0 < y_1)
1382 {
1383 lo_y = ceil(y_0 - 0.5);
1384 hi_y = floor(y_1 - 0.5);
1385
1386 if (hi_y == (y_1 - 0.5))
1387 hi_y -= 1;
1388 }
1389 else
1390 {
1391 lo_y = ceil(y_1 - 0.5);
1392 hi_y = floor(y_0 - 0.5);
1393
1394 if (hi_y == (y_0 - 0.5))
1395 hi_y -= 1;
1396 }
1397
1398 if (lo_y < 0) lo_y = 0;
1399 if (hi_y >= hght) hi_y = hght-1;
1400
1401 if (lo_y > hi_y)
1402 return; /* no scan lines crossed */
1403
1404 if (lo_y < min_y) min_y = lo_y;
1405 if (hi_y > max_y) max_y = hi_y;
1406
1407 x_delta = (x_1 - x_0) / (y_1 - y_0);
1408 x_value = x_0 + x_delta * ((lo_y + 0.5) - y_0);
1409
1410 for (i=lo_y; i<=hi_y; i++)
1411 {
1412 *((double *) extarr_next(scanbuf[i])) = x_value;
1413 x_value += x_delta;
1414 }
1415 }
1416
1417
get_scanbits(val)1418 void get_scanbits(val)
1419 int val;
1420 {
1421 int i, j;
1422 int lo_x, hi_x;
1423 unsigned nvals;
1424 double *vals;
1425 ScanBit *scanbit;
1426
1427 for (i=min_y; i<=max_y; i++)
1428 {
1429 vals = (double *) scanbuf[i]->body;
1430 nvals = scanbuf[i]->count;
1431 assert((nvals % 2) == 0);
1432 qsort(vals, nvals, sizeof(double), double_comp);
1433
1434 for (j=0; j<nvals; j+=2)
1435 {
1436 lo_x = ceil(vals[j] - 0.5);
1437 hi_x = floor(vals[j+1] - 0.5);
1438
1439 if (lo_x < 0) lo_x = 0;
1440 if (hi_x >= wdth) hi_x = wdth-1;
1441
1442 if (lo_x <= hi_x)
1443 {
1444 scanbit = (ScanBit *) extarr_next(scanbits);
1445 scanbit->y = i;
1446 scanbit->lo_x = lo_x;
1447 scanbit->hi_x = hi_x;
1448 scanbit->val = val;
1449 }
1450 }
1451
1452 scanbuf[i]->count = 0;
1453 }
1454 }
1455