1 /*
2 * render.c
3 * kirk johnson
4 * october 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 static void new_stars _P((double));
50 static void new_grid _P((int, int));
51 static void new_grid_dot _P((double *, double *));
52 static int dot_comp _P((const void *, const void *));
53 static void render_rows_setup _P((void));
54 static void render_next_row _P((s8or32 *, int));
55 static void no_shade_row _P((s8or32 *, u_char *));
56 static void compute_sun_vector _P((double *));
57 static void orth_compute_inv_x _P((double *));
58 static void orth_shade_row _P((int, s8or32 *, double *, double *, u_char *));
59 static void merc_shade_row _P((int, s8or32 *, double *, u_char *));
60 static void cyl_shade_row _P((int, s8or32 *, double *, u_char *));
61
62 static int scanbitcnt;
63 static ScanBit *scanbit;
64 static s8or32 scan_to_pix[256];
65
66 static int night_val;
67 static int day_val_base;
68 static double day_val_delta;
69
70 static ExtArr dots = NULL;
71 static int dotcnt;
72 static ScanDot *dot;
73
74
dot_comp(a,b)75 static int dot_comp(a, b)
76 const void *a;
77 const void *b;
78 {
79 return (((const ScanDot *) a)->y - ((const ScanDot *) b)->y);
80 }
81
82
render_rows_setup()83 static void render_rows_setup()
84 {
85 int i;
86
87 scanbitcnt = scanbits->count;
88 scanbit = (ScanBit *) scanbits->body;
89 dotcnt = dots->count;
90 dot = (ScanDot *) dots->body;
91
92 /* precompute table for translating between
93 * scan buffer values and pixel types
94 */
95 for (i=0; i<256; i++)
96 if (i == 0)
97 scan_to_pix[i] = PixTypeSpace;
98 else if (i > 64)
99 scan_to_pix[i] = PixTypeLand;
100 else
101 scan_to_pix[i] = PixTypeWater;
102 }
103
104
render_next_row(buf,idx)105 static void render_next_row(buf, idx)
106 s8or32 *buf;
107 int idx;
108 {
109 int i, i_lim;
110 int tmp;
111 int _scanbitcnt;
112 ScanBit *_scanbit;
113
114 xearth_bzero((char *) buf, (unsigned) (sizeof(s8or32) * wdth));
115
116 /* explicitly copy scanbitcnt and scanbit to local variables
117 * to help compilers figure out that they can be registered
118 */
119 _scanbitcnt = scanbitcnt;
120 _scanbit = scanbit;
121
122 while ((_scanbitcnt > 0) && (_scanbit->y == idx))
123 {
124 /* use i_lim to encourage compilers to register loop limit
125 */
126 i_lim = _scanbit->hi_x;
127 tmp = _scanbit->val;
128 for (i=_scanbit->lo_x; i<=i_lim; i++)
129 buf[i] += tmp;
130
131 _scanbit += 1;
132 _scanbitcnt -= 1;
133 }
134
135 /* copy changes to scanbitcnt and scanbit out to memory
136 */
137 scanbitcnt = _scanbitcnt;
138 scanbit = _scanbit;
139
140 /* use i_lim to encourage compilers to register loop limit
141 */
142 i_lim = wdth;
143 for (i=0; i<i_lim; i++)
144 buf[i] = scan_to_pix[(int) (buf[i] & 0xff)];
145
146 while ((dotcnt > 0) && (dot->y == idx))
147 {
148 tmp = dot->x;
149
150 if (dot->type == DotTypeStar)
151 {
152 if (buf[tmp] == PixTypeSpace)
153 buf[tmp] = PixTypeStar;
154 }
155 else
156 {
157 switch (buf[tmp])
158 {
159 case PixTypeLand:
160 buf[tmp] = PixTypeGridLand;
161 break;
162
163 case PixTypeWater:
164 buf[tmp] = PixTypeGridWater;
165 break;
166 }
167 }
168
169 dot += 1;
170 dotcnt -= 1;
171 }
172 }
173
174
no_shade_row(scanbuf,rslt)175 static void no_shade_row(scanbuf, rslt)
176 s8or32 *scanbuf;
177 u_char *rslt;
178 {
179 int i, i_lim;
180
181 /* use i_lim to encourage compilers to register loop limit
182 */
183 i_lim = wdth;
184 for (i=0; i<i_lim; i++)
185 {
186 switch (scanbuf[i])
187 {
188 case PixTypeSpace: /* black */
189 rslt[0] = 0;
190 rslt[1] = 0;
191 rslt[2] = 0;
192 break;
193
194 case PixTypeStar: /* white */
195 case PixTypeGridLand:
196 case PixTypeGridWater:
197 rslt[0] = 255;
198 rslt[1] = 255;
199 rslt[2] = 255;
200 break;
201
202 case PixTypeLand: /* green */
203 rslt[0] = 0;
204 rslt[1] = 255;
205 rslt[2] = 0;
206 break;
207
208 case PixTypeWater: /* blue */
209 rslt[0] = 0;
210 rslt[1] = 0;
211 rslt[2] = 255;
212 break;
213
214 default:
215 assert(0);
216 }
217
218 rslt += 3;
219 }
220 }
221
222
compute_sun_vector(rslt)223 static void compute_sun_vector(rslt)
224 double *rslt;
225 {
226 rslt[0] = sin(sun_lon * (M_PI/180)) * cos(sun_lat * (M_PI/180));
227 rslt[1] = sin(sun_lat * (M_PI/180));
228 rslt[2] = cos(sun_lon * (M_PI/180)) * cos(sun_lat * (M_PI/180));
229
230 XFORM_ROTATE(rslt, view_pos_info);
231 }
232
233
orth_compute_inv_x(inv_x)234 static void orth_compute_inv_x(inv_x)
235 double *inv_x;
236 {
237 int i, i_lim;
238
239 i_lim = wdth;
240 for (i=0; i<i_lim; i++)
241 inv_x[i] = INV_XPROJECT(i);
242 }
243
244
orth_shade_row(idx,scanbuf,sol,inv_x,rslt)245 static void orth_shade_row(idx, scanbuf, sol, inv_x, rslt)
246 int idx;
247 s8or32 *scanbuf;
248 double *sol;
249 double *inv_x;
250 u_char *rslt;
251 {
252 int i, i_lim;
253 int scanbuf_val;
254 int val;
255 double x, y, z;
256 double scale;
257 double tmp;
258 double y_sol_1;
259
260 y = INV_YPROJECT(idx);
261
262 /* save a little computation in the inner loop
263 */
264 tmp = 1 - (y*y);
265 y_sol_1 = y * sol[1];
266
267 /* use i_lim to encourage compilers to register loop limit
268 */
269 i_lim = wdth;
270 for (i=0; i<i_lim; i++)
271 {
272 scanbuf_val = scanbuf[i];
273
274 switch (scanbuf_val)
275 {
276 case PixTypeSpace: /* black */
277 rslt[0] = 0;
278 rslt[1] = 0;
279 rslt[2] = 0;
280 break;
281
282 case PixTypeStar: /* white */
283 case PixTypeGridLand:
284 case PixTypeGridWater:
285 rslt[0] = 255;
286 rslt[1] = 255;
287 rslt[2] = 255;
288 break;
289
290 case PixTypeLand: /* green, blue */
291 case PixTypeWater:
292 x = inv_x[i];
293 z = tmp - (x*x);
294 z = SQRT(z);
295 scale = (x * sol[0]) + y_sol_1 + (z * sol[2]);
296 if (scale < 0)
297 {
298 val = night_val;
299 }
300 else
301 {
302 val = day_val_base + (scale * day_val_delta);
303 if (val > 255)
304 val = 255;
305 else
306 assert(val >= 0);
307 }
308
309 if (scanbuf_val == PixTypeLand)
310 {
311 /* land (green) */
312 rslt[0] = 0;
313 rslt[1] = val;
314 rslt[2] = 0;
315 }
316 else
317 {
318 /* water (blue) */
319 rslt[0] = 0;
320 rslt[1] = 0;
321 rslt[2] = val;
322 }
323 break;
324
325 default:
326 assert(0);
327 }
328
329 rslt += 3;
330 }
331 }
332
333
merc_shade_row(idx,scanbuf,sol,rslt)334 static void merc_shade_row(idx, scanbuf, sol, rslt)
335 int idx;
336 s8or32 *scanbuf;
337 double *sol;
338 u_char *rslt;
339 {
340 int i, i_lim;
341 int scanbuf_val;
342 int val;
343 double x, y, z;
344 double sin_theta;
345 double cos_theta;
346 double scale;
347 double tmp;
348 double y_sol_1;
349
350 y = INV_YPROJECT(idx);
351 y = INV_MERCATOR_Y(y);
352
353 /* conceptually, on each iteration of the i loop, we want:
354 *
355 * x = sin(INV_XPROJECT(i)) * sqrt(1 - (y*y));
356 * z = cos(INV_XPROJECT(i)) * sqrt(1 - (y*y));
357 *
358 * computing this directly is rather expensive, however, so we only
359 * compute the first (i=0) pair of values directly; all other pairs
360 * (i>0) are obtained through successive rotations of the original
361 * pair (by inv_proj_scale radians).
362 */
363
364 /* compute initial (x, z) values
365 */
366 tmp = sqrt(1 - (y*y));
367 x = sin(INV_XPROJECT(0)) * tmp;
368 z = cos(INV_XPROJECT(0)) * tmp;
369
370 /* compute rotation coefficients used
371 * to find subsequent (x, z) values
372 */
373 tmp = proj_info.inv_proj_scale;
374 sin_theta = sin(tmp);
375 cos_theta = cos(tmp);
376
377 /* save a little computation in the inner loop
378 */
379 y_sol_1 = y * sol[1];
380
381 /* use i_lim to encourage compilers to register loop limit
382 */
383 i_lim = wdth;
384 for (i=0; i<i_lim; i++)
385 {
386 scanbuf_val = scanbuf[i];
387
388 switch (scanbuf_val)
389 {
390 case PixTypeSpace: /* black */
391 rslt[0] = 0;
392 rslt[1] = 0;
393 rslt[2] = 0;
394 break;
395
396 case PixTypeStar: /* white */
397 case PixTypeGridLand:
398 case PixTypeGridWater:
399 rslt[0] = 255;
400 rslt[1] = 255;
401 rslt[2] = 255;
402 break;
403
404 case PixTypeLand: /* green, blue */
405 case PixTypeWater:
406 scale = (x * sol[0]) + y_sol_1 + (z * sol[2]);
407 if (scale < 0)
408 {
409 val = night_val;
410 }
411 else
412 {
413 val = day_val_base + (scale * day_val_delta);
414 if (val > 255)
415 val = 255;
416 else
417 assert(val >= 0);
418 }
419
420 if (scanbuf_val == PixTypeLand)
421 {
422 /* land (green) */
423 rslt[0] = 0;
424 rslt[1] = val;
425 rslt[2] = 0;
426 }
427 else
428 {
429 /* water (blue) */
430 rslt[0] = 0;
431 rslt[1] = 0;
432 rslt[2] = val;
433 }
434 break;
435
436 default:
437 assert(0);
438 }
439
440 /* compute next (x, z) values via 2-d rotation
441 */
442 tmp = (cos_theta * z) - (sin_theta * x);
443 x = (sin_theta * z) + (cos_theta * x);
444 z = tmp;
445
446 rslt += 3;
447 }
448 }
449
450
cyl_shade_row(idx,scanbuf,sol,rslt)451 static void cyl_shade_row(idx, scanbuf, sol, rslt)
452 int idx;
453 s8or32 *scanbuf;
454 double *sol;
455 u_char *rslt;
456 {
457 int i, i_lim;
458 int scanbuf_val;
459 int val;
460 double x, y, z;
461 double sin_theta;
462 double cos_theta;
463 double scale;
464 double tmp;
465 double y_sol_1;
466
467 y = INV_YPROJECT(idx);
468 y = INV_CYLINDRICAL_Y(y);
469
470 /* conceptually, on each iteration of the i loop, we want:
471 *
472 * x = sin(INV_XPROJECT(i)) * sqrt(1 - (y*y));
473 * z = cos(INV_XPROJECT(i)) * sqrt(1 - (y*y));
474 *
475 * computing this directly is rather expensive, however, so we only
476 * compute the first (i=0) pair of values directly; all other pairs
477 * (i>0) are obtained through successive rotations of the original
478 * pair (by inv_proj_scale radians).
479 */
480
481 /* compute initial (x, z) values
482 */
483 tmp = sqrt(1 - (y*y));
484 x = sin(INV_XPROJECT(0)) * tmp;
485 z = cos(INV_XPROJECT(0)) * tmp;
486
487 /* compute rotation coefficients used
488 * to find subsequent (x, z) values
489 */
490 tmp = proj_info.inv_proj_scale;
491 sin_theta = sin(tmp);
492 cos_theta = cos(tmp);
493
494 /* save a little computation in the inner loop
495 */
496 y_sol_1 = y * sol[1];
497
498 /* use i_lim to encourage compilers to register loop limit
499 */
500 i_lim = wdth;
501 for (i=0; i<i_lim; i++)
502 {
503 scanbuf_val = scanbuf[i];
504
505 switch (scanbuf_val)
506 {
507 case PixTypeSpace: /* black */
508 rslt[0] = 0;
509 rslt[1] = 0;
510 rslt[2] = 0;
511 break;
512
513 case PixTypeStar: /* white */
514 case PixTypeGridLand:
515 case PixTypeGridWater:
516 rslt[0] = 255;
517 rslt[1] = 255;
518 rslt[2] = 255;
519 break;
520
521 case PixTypeLand: /* green, blue */
522 case PixTypeWater:
523 scale = (x * sol[0]) + y_sol_1 + (z * sol[2]);
524 if (scale < 0)
525 {
526 val = night_val;
527 }
528 else
529 {
530 val = day_val_base + (scale * day_val_delta);
531 if (val > 255)
532 val = 255;
533 else
534 assert(val >= 0);
535 }
536
537 if (scanbuf_val == PixTypeLand)
538 {
539 /* land (green) */
540 rslt[0] = 0;
541 rslt[1] = val;
542 rslt[2] = 0;
543 }
544 else
545 {
546 /* water (blue) */
547 rslt[0] = 0;
548 rslt[1] = 0;
549 rslt[2] = val;
550 }
551 break;
552
553 default:
554 assert(0);
555 }
556
557 /* compute next (x, z) values via 2-d rotation
558 */
559 tmp = (cos_theta * z) - (sin_theta * x);
560 x = (sin_theta * z) + (cos_theta * x);
561 z = tmp;
562
563 rslt += 3;
564 }
565 }
566
567
568 void render(rowfunc)
569 int (*rowfunc) _P((u_char *));
570 {
571 int i, i_lim;
572 s8or32 *scanbuf;
573 u_char *row;
574 double *inv_x;
575 double sol[3];
576 double tmp;
577
578 scanbuf = (s8or32 *) malloc((unsigned) (sizeof(s8or32) * wdth));
579 row = (u_char *) malloc((unsigned) wdth*3);
580 assert((scanbuf != NULL) && (row != NULL));
581
582 inv_x = NULL;
583 render_rows_setup();
584
585 if (do_shade)
586 {
587 /* inv_x[] only gets used with orthographic projection
588 */
589 if (proj_type == ProjTypeOrthographic)
590 {
591 inv_x = (double *) malloc((unsigned) sizeof(double) * wdth);
592 assert(inv_x != NULL);
593 orth_compute_inv_x(inv_x);
594 }
595
596 compute_sun_vector(sol);
597
598 /* precompute shading parameters
599 */
600 night_val = night * (255.99/100.0);
601 tmp = terminator / 100.0;
602 day_val_base = ((tmp * day) + ((1-tmp) * night)) * (255.99/100.0);
603 day_val_delta = (day * (255.99/100.0)) - day_val_base;
604 }
605
606 /* main render loop
607 * (use i_lim to encourage compilers to register loop limit)
608 */
609 i_lim = hght;
610 for (i=0; i<i_lim; i++)
611 {
612 render_next_row(scanbuf, i);
613
614 if (!do_shade)
615 no_shade_row(scanbuf, row);
616 else if (proj_type == ProjTypeOrthographic)
617 orth_shade_row(i, scanbuf, sol, inv_x, row);
618 else if (proj_type == ProjTypeMercator)
619 merc_shade_row(i, scanbuf, sol, row);
620 else /* (proj_type == ProjTypeCylindrical) */
621 cyl_shade_row(i, scanbuf, sol, row);
622
623 rowfunc(row);
624 }
625
626 free(scanbuf);
627 free(row);
628
629 if (inv_x != NULL) free(inv_x);
630 }
631
632
do_dots()633 void do_dots()
634 {
635 if (dots == NULL)
636 dots = extarr_alloc(sizeof(ScanDot));
637 else
638 dots->count = 0;
639
640 if (do_stars) new_stars(star_freq);
641 if (do_grid) new_grid(grid_big, grid_small);
642
643 qsort(dots->body, dots->count, sizeof(ScanDot), dot_comp);
644 }
645
646
new_stars(freq)647 static void new_stars(freq)
648 double freq;
649 {
650 int i;
651 int x, y;
652 int max_stars;
653 ScanDot *newdot;
654
655 max_stars = wdth * hght * freq;
656
657 for (i=0; i<max_stars; i++)
658 {
659 x = random() % wdth;
660 y = random() % hght;
661
662 newdot = (ScanDot *) extarr_next(dots);
663 newdot->x = x;
664 newdot->y = y;
665 newdot->type = DotTypeStar;
666
667 if ((big_stars) && (x+1 < wdth) && ((random() % 100) < big_stars))
668 {
669 newdot = (ScanDot *) extarr_next(dots);
670 newdot->x = x+1;
671 newdot->y = y;
672 newdot->type = DotTypeStar;
673 }
674 }
675 }
676
677
new_grid(big,small)678 static void new_grid(big, small)
679 int big;
680 int small;
681 {
682 int i, j;
683 int cnt;
684 double lat, lon;
685 double lat_scale, lon_scale;
686 double cs_lat[2];
687 double cs_lon[2];
688
689 /* lines of longitude
690 */
691 lon_scale = M_PI / (2 * big);
692 lat_scale = M_PI / (2 * big * small);
693 for (i=(-2*big); i<(2*big); i++)
694 {
695 lon = i * lon_scale;
696 cs_lon[0] = cos(lon);
697 cs_lon[1] = sin(lon);
698
699 for (j=(-(big*small)+1); j<(big*small); j++)
700 {
701 lat = j * lat_scale;
702 cs_lat[0] = cos(lat);
703 cs_lat[1] = sin(lat);
704
705 new_grid_dot(cs_lat, cs_lon);
706 }
707 }
708
709 /* lines of latitude
710 */
711 lat_scale = M_PI / (2 * big);
712 for (i=(1-big); i<big; i++)
713 {
714 lat = i * lat_scale;
715 cs_lat[0] = cos(lat);
716 cs_lat[1] = sin(lat);
717 cnt = 2 * ((int) ((cs_lat[0] * small) + 0.5)) * big;
718 lon_scale = M_PI / cnt;
719
720 for (j=(-cnt); j<cnt; j++)
721 {
722 lon = j * lon_scale;
723 cs_lon[0] = cos(lon);
724 cs_lon[1] = sin(lon);
725
726 new_grid_dot(cs_lat, cs_lon);
727 }
728 }
729 }
730
731
new_grid_dot(cs_lat,cs_lon)732 static void new_grid_dot(cs_lat, cs_lon)
733 double *cs_lat;
734 double *cs_lon;
735 {
736 int x, y;
737 double pos[3];
738 ScanDot *new;
739
740 pos[0] = cs_lon[1] * cs_lat[0];
741 pos[1] = cs_lat[1];
742 pos[2] = cs_lon[0] * cs_lat[0];
743
744 XFORM_ROTATE(pos, view_pos_info);
745
746 if (proj_type == ProjTypeOrthographic)
747 {
748 /* if the grid dot isn't visible, return immediately
749 */
750 if (pos[2] <= 0) return;
751 }
752 else if (proj_type == ProjTypeMercator)
753 {
754 /* apply mercator projection
755 */
756 pos[0] = MERCATOR_X(pos[0], pos[2]);
757 pos[1] = MERCATOR_Y(pos[1]);
758 }
759 else /* (proj_type == ProjTypeCylindrical) */
760 {
761 /* apply cylindrical projection
762 */
763 pos[0] = CYLINDRICAL_X(pos[0], pos[2]);
764 pos[1] = CYLINDRICAL_Y(pos[1]);
765 }
766
767 x = XPROJECT(pos[0]);
768 y = YPROJECT(pos[1]);
769
770 if ((x >= 0) && (x < wdth) && (y >= 0) && (y < hght))
771 {
772 new = (ScanDot *) extarr_next(dots);
773 new->x = x;
774 new->y = y;
775 new->type = DotTypeGrid;
776 }
777 }
778