/*
* scan.c
* kirk johnson
* august 1993
*
* Copyright (C) 1989, 1990, 1993-1995, 1999 Kirk Lauritz Johnson
*
* Parts of the source code (as marked) are:
* Copyright (C) 1989, 1990, 1991 by Jim Frost
* Copyright (C) 1992 by Jamie Zawinski <jwz@lucid.com>
*
* Permission to use, copy, modify and freely distribute xearth for
* non-commercial and not-for-profit purposes is hereby granted
* without fee, provided that both the above copyright notice and this
* permission notice appear in all copies and in supporting
* documentation.
*
* Unisys Corporation holds worldwide patent rights on the Lempel Zev
* Welch (LZW) compression technique employed in the CompuServe GIF
* image file format as well as in other formats. Unisys has made it
* clear, however, that it does not require licensing or fees to be
* paid for freely distributed, non-commercial applications (such as
* xearth) that employ LZW/GIF technology. Those wishing further
* information about licensing the LZW patent should contact Unisys
* directly at (lzw_info@unisys.com) or by writing to
*
* Unisys Corporation
* Welch Licensing Department
* M/S-C1SW19
* P.O. Box 500
* Blue Bell, PA 19424
*
* The author makes no representations about the suitability of this
* software for any purpose. It is provided "as is" without express or
* implied warranty.
*
* THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
* INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS,
* IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT
* OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
* LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT,
* NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
* CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
*/
#include "xearth.h"
#include "kljcpyrt.h"
#define MAP_DATA_SCALE (30000)
#define XingTypeEntry (0)
#define XingTypeExit (1)
typedef struct
{
int type;
int cidx;
double x, y;
double angle;
} EdgeXing;
void scan_map _P((void));
void orth_scan_outline _P((void));
void orth_scan_curves _P((void));
double *orth_extract_curve _P((int, short *));
void orth_scan_along_curve _P((double *, double *, int));
void orth_find_edge_xing _P((double *, double *, double *));
void orth_handle_xings _P((void));
void orth_scan_arc _P((double, double, double, double, double, double));
void merc_scan_outline _P((void));
void merc_scan_curves _P((void));
double *merc_extract_curve _P((int, short *));
void merc_scan_along_curve _P((double *, double *, int));
double merc_find_edge_xing _P((double *, double *));
void merc_handle_xings _P((void));
void merc_scan_edge _P((EdgeXing *, EdgeXing *));
void cyl_scan_outline _P((void));
void cyl_scan_curves _P((void));
double *cyl_extract_curve _P((int, short *));
void cyl_scan_along_curve _P((double *, double *, int));
double cyl_find_edge_xing _P((double *, double *));
void cyl_handle_xings _P((void));
void cyl_scan_edge _P((EdgeXing *, EdgeXing *));
void xing_error _P((const char *, int, int, int, EdgeXing *));
void scan _P((double, double, double, double));
void get_scanbits _P((int));
static int double_comp _P((const void *, const void *));
static int scanbit_comp _P((const void *, const void *));
static int orth_edgexing_comp _P((const void *, const void *));
static int merc_edgexing_comp _P((const void *, const void *));
static int cyl_edgexing_comp _P((const void *, const void *));
ViewPosInfo view_pos_info;
ProjInfo proj_info;
int first_scan = 1;
ExtArr scanbits;
ExtArr edgexings;
int min_y, max_y;
ExtArr *scanbuf;
static int double_comp(a, b)
const void *a;
const void *b;
{
double val_a;
double val_b;
int rslt;
val_a = *((const double *) a);
val_b = *((const double *) b);
if (val_a < val_b)
rslt = -1;
else if (val_a > val_b)
rslt = 1;
else
rslt = 0;
return rslt;
}
static int scanbit_comp(a, b)
const void *a;
const void *b;
{
return (((const ScanBit *) a)->y - ((const ScanBit *) b)->y);
}
static int orth_edgexing_comp(a, b)
const void *a;
const void *b;
{
double angle_a;
double angle_b;
int rslt;
angle_a = ((const EdgeXing *) a)->angle;
angle_b = ((const EdgeXing *) b)->angle;
if (angle_a < angle_b)
rslt = -1;
else if (angle_a > angle_b)
rslt = 1;
else
rslt = 0;
return rslt;
}
static int merc_edgexing_comp(a, b)
const void *a;
const void *b;
{
double val_a;
double val_b;
int rslt;
val_a = ((const EdgeXing *) a)->angle;
val_b = ((const EdgeXing *) b)->angle;
if (val_a < val_b)
{
rslt = -1;
}
else if (val_a > val_b)
{
rslt = 1;
}
else if (val_a == 0)
{
val_a = ((const EdgeXing *) a)->y;
val_b = ((const EdgeXing *) b)->y;
if (val_a < val_b)
rslt = -1;
else if (val_a > val_b)
rslt = 1;
else
rslt = 0;
}
else if (val_a == 2)
{
val_a = ((const EdgeXing *) a)->y;
val_b = ((const EdgeXing *) b)->y;
if (val_a > val_b)
rslt = -1;
else if (val_a < val_b)
rslt = 1;
else
rslt = 0;
}
else
{
/* keep lint happy */
rslt = 0;
assert(0);
}
return rslt;
}
static int cyl_edgexing_comp(a, b)
const void *a;
const void *b;
{
double val_a;
double val_b;
int rslt;
val_a = ((const EdgeXing *) a)->angle;
val_b = ((const EdgeXing *) b)->angle;
if (val_a < val_b)
{
rslt = -1;
}
else if (val_a > val_b)
{
rslt = 1;
}
else if (val_a == 0)
{
val_a = ((const EdgeXing *) a)->y;
val_b = ((const EdgeXing *) b)->y;
if (val_a < val_b)
rslt = -1;
else if (val_a > val_b)
rslt = 1;
else
rslt = 0;
}
else if (val_a == 2)
{
val_a = ((const EdgeXing *) a)->y;
val_b = ((const EdgeXing *) b)->y;
if (val_a > val_b)
rslt = -1;
else if (val_a < val_b)
rslt = 1;
else
rslt = 0;
}
else
{
/* keep lint happy */
rslt = 0;
assert(0);
}
return rslt;
}
void scan_map()
{
int i;
ViewPosInfo *vpi;
ProjInfo *pi;
vpi = &view_pos_info;
vpi->cos_lat = cos(view_lat * (M_PI/180));
vpi->sin_lat = sin(view_lat * (M_PI/180));
vpi->cos_lon = cos(view_lon * (M_PI/180));
vpi->sin_lon = sin(view_lon * (M_PI/180));
vpi->cos_rot = cos(view_rot * (M_PI/180));
vpi->sin_rot = sin(view_rot * (M_PI/180));
pi = &proj_info;
if (proj_type == ProjTypeOrthographic)
{
pi->proj_scale = ((hght < wdth) ? hght : wdth) * (view_mag / 2) * 0.99;
}
else
{
/* proj_type is either ProjTypeMercator or ProjTypeCylindrical
*/
pi->proj_scale = (view_mag * wdth) / (2 * M_PI);
}
pi->proj_xofs = (double) wdth / 2 + shift_x;
pi->proj_yofs = (double) hght / 2 + shift_y;
pi->inv_proj_scale = 1 / pi->proj_scale;
/* the first time through, allocate scanbits and edgexings;
* on subsequent passes, simply reset them.
*/
if (first_scan)
{
scanbits = extarr_alloc(sizeof(ScanBit));
edgexings = extarr_alloc(sizeof(EdgeXing));
}
else
{
scanbits->count = 0;
edgexings->count = 0;
}
/* maybe only allocate these once and reset them on
* subsequent passes (like scanbits and edgexings)?
*/
scanbuf = (ExtArr *) malloc((unsigned) sizeof(ExtArr) * hght);
assert(scanbuf != NULL);
for (i=0; i<hght; i++)
scanbuf[i] = extarr_alloc(sizeof(double));
if (proj_type == ProjTypeOrthographic)
{
orth_scan_outline();
orth_scan_curves();
}
else if (proj_type == ProjTypeMercator)
{
merc_scan_outline();
merc_scan_curves();
}
else /* (proj_type == ProjTypeCylindrical) */
{
cyl_scan_outline();
cyl_scan_curves();
}
for (i=0; i<hght; i++)
extarr_free(scanbuf[i]);
free(scanbuf);
qsort(scanbits->body, scanbits->count, sizeof(ScanBit), scanbit_comp);
first_scan = 0;
}
void orth_scan_outline()
{
min_y = hght;
max_y = -1;
orth_scan_arc(1.0, 0.0, 0.0, 1.0, 0.0, (2*M_PI));
get_scanbits(64);
}
void orth_scan_curves()
{
int i;
int cidx;
int npts;
int val;
short *raw;
double *pos;
double *prev;
double *curr;
cidx = 0;
raw = map_data;
while (1)
{
npts = raw[0];
if (npts == 0) break;
val = raw[1];
raw += 2;
pos = orth_extract_curve(npts, raw);
prev = pos + (npts-1)*3;
curr = pos;
min_y = hght;
max_y = -1;
for (i=0; i<npts; i++)
{
orth_scan_along_curve(prev, curr, cidx);
prev = curr;
curr += 3;
}
free(pos);
if (edgexings->count > 0)
orth_handle_xings();
if (min_y <= max_y)
get_scanbits(val);
cidx += 1;
raw += 3*npts;
}
}
double *orth_extract_curve(npts, data)
int npts;
short *data;
{
int i;
int x, y, z;
double scale;
double *pos;
double *rslt;
rslt = (double *) malloc((unsigned) sizeof(double) * 3 * npts);
assert(rslt != NULL);
x = 0;
y = 0;
z = 0;
scale = 1.0 / MAP_DATA_SCALE;
pos = rslt;
for (i=0; i<npts; i++)
{
x += data[0];
y += data[1];
z += data[2];
pos[0] = x * scale;
pos[1] = y * scale;
pos[2] = z * scale;
XFORM_ROTATE(pos, view_pos_info);
data += 3;
pos += 3;
}
return rslt;
}
void orth_scan_along_curve(prev, curr, cidx)
double *prev;
double *curr;
int cidx;
{
double extra[3];
EdgeXing *xing;
if (prev[2] <= 0) /* prev not visible */
{
if (curr[2] <= 0)
return; /* neither point visible */
orth_find_edge_xing(prev, curr, extra);
/* extra[] is an edge crossing (entry point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeEntry;
xing->cidx = cidx;
xing->x = extra[0];
xing->y = extra[1];
xing->angle = atan2(extra[1], extra[0]);
prev = extra;
}
else if (curr[2] <= 0) /* curr not visible */
{
orth_find_edge_xing(prev, curr, extra);
/* extra[] is an edge crossing (exit point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeExit;
xing->cidx = cidx;
xing->x = extra[0];
xing->y = extra[1];
xing->angle = atan2(extra[1], extra[0]);
curr = extra;
}
scan(XPROJECT(prev[0]), YPROJECT(prev[1]),
XPROJECT(curr[0]), YPROJECT(curr[1]));
}
void orth_find_edge_xing(prev, curr, rslt)
double *prev;
double *curr;
double *rslt;
{
double tmp;
double r0, r1;
tmp = curr[2] / (curr[2] - prev[2]);
r0 = curr[0] - tmp * (curr[0] - prev[0]);
r1 = curr[1] - tmp * (curr[1] - prev[1]);
tmp = sqrt((r0*r0) + (r1*r1));
rslt[0] = r0 / tmp;
rslt[1] = r1 / tmp;
rslt[2] = 0;
}
void orth_handle_xings()
{
int i;
int nxings;
EdgeXing *xings;
EdgeXing *from;
EdgeXing *to;
xings = (EdgeXing *) edgexings->body;
nxings = edgexings->count;
assert((nxings % 2) == 0);
qsort(xings, (unsigned) nxings, sizeof(EdgeXing), orth_edgexing_comp);
if (xings[0].type == XingTypeExit)
{
for (i=0; i<nxings; i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
orth_scan_arc(from->x, from->y, from->angle,
to->x, to->y, to->angle);
}
}
else
{
from = &(xings[nxings-1]);
to = &(xings[0]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry) ||
(from->angle < to->angle))
xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
orth_scan_arc(from->x, from->y, from->angle,
to->x, to->y, to->angle+(2*M_PI));
for (i=1; i<(nxings-1); i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
orth_scan_arc(from->x, from->y, from->angle,
to->x, to->y, to->angle);
}
}
edgexings->count = 0;
}
void orth_scan_arc(x_0, y_0, a_0, x_1, y_1, a_1)
double x_0, y_0, a_0;
double x_1, y_1, a_1;
{
int i;
int lo, hi;
double angle, step;
double prev_x, prev_y;
double curr_x, curr_y;
double c_step, s_step;
double arc_x, arc_y;
double tmp;
assert(a_0 < a_1);
step = proj_info.inv_proj_scale * 10;
if (step > 0.05) step = 0.05;
lo = ceil(a_0 / step);
hi = floor(a_1 / step);
prev_x = XPROJECT(x_0);
prev_y = YPROJECT(y_0);
if (lo <= hi)
{
c_step = cos(step);
s_step = sin(step);
angle = lo * step;
arc_x = cos(angle);
arc_y = sin(angle);
for (i=lo; i<=hi; i++)
{
curr_x = XPROJECT(arc_x);
curr_y = YPROJECT(arc_y);
scan(prev_x, prev_y, curr_x, curr_y);
/* instead of repeatedly calling cos() and sin() to get the next
* values for arc_x and arc_y, simply rotate the existing values
*/
tmp = (c_step * arc_x) - (s_step * arc_y);
arc_y = (s_step * arc_x) + (c_step * arc_y);
arc_x = tmp;
prev_x = curr_x;
prev_y = curr_y;
}
}
curr_x = XPROJECT(x_1);
curr_y = YPROJECT(y_1);
scan(prev_x, prev_y, curr_x, curr_y);
}
void merc_scan_outline()
{
double left, right;
double top, bottom;
min_y = hght;
max_y = -1;
left = XPROJECT(-M_PI);
right = XPROJECT(M_PI);
top = YPROJECT(BigNumber);
bottom = YPROJECT(-BigNumber);
scan(right, top, left, top);
scan(left, top, left, bottom);
scan(left, bottom, right, bottom);
scan(right, bottom, right, top);
get_scanbits(64);
}
void merc_scan_curves()
{
int i;
int cidx;
int npts;
int val;
short *raw;
double *pos;
double *prev;
double *curr;
cidx = 0;
raw = map_data;
while (1)
{
npts = raw[0];
if (npts == 0) break;
val = raw[1];
raw += 2;
pos = merc_extract_curve(npts, raw);
prev = pos + (npts-1)*5;
curr = pos;
min_y = hght;
max_y = -1;
for (i=0; i<npts; i++)
{
merc_scan_along_curve(prev, curr, cidx);
prev = curr;
curr += 5;
}
free(pos);
if (edgexings->count > 0)
merc_handle_xings();
if (min_y <= max_y)
get_scanbits(val);
cidx += 1;
raw += 3*npts;
}
}
double *merc_extract_curve(npts, data)
int npts;
short *data;
{
int i;
int x, y, z;
double scale;
double *pos;
double *rslt;
rslt = (double *) malloc((unsigned) sizeof(double) * 5 * npts);
assert(rslt != NULL);
x = 0;
y = 0;
z = 0;
scale = 1.0 / MAP_DATA_SCALE;
pos = rslt;
for (i=0; i<npts; i++)
{
x += data[0];
y += data[1];
z += data[2];
pos[0] = x * scale;
pos[1] = y * scale;
pos[2] = z * scale;
XFORM_ROTATE(pos, view_pos_info);
/* apply mercator projection
*/
pos[3] = MERCATOR_X(pos[0], pos[2]);
pos[4] = MERCATOR_Y(pos[1]);
data += 3;
pos += 5;
}
return rslt;
}
void merc_scan_along_curve(prev, curr, cidx)
double *prev;
double *curr;
int cidx;
{
double px, py;
double cx, cy;
double dx;
double mx, my;
EdgeXing *xing;
px = prev[3];
cx = curr[3];
py = prev[4];
cy = curr[4];
dx = cx - px;
if (dx > 0)
{
/* curr to the right of prev
*/
if (dx > ((2*M_PI) - dx))
{
/* vertical edge crossing to the left of prev
*/
/* find exit point (left edge) */
mx = - M_PI;
my = merc_find_edge_xing(prev, curr);
/* scan from prev to exit point */
scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
/* (mx, my) is an edge crossing (exit point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeExit;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 2; /* left edge */
/* scan from entry point (right edge) to curr */
mx = M_PI;
scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
/* (mx, my) is an edge crossing (entry point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeEntry;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 0; /* right edge */
}
else
{
/* no vertical edge crossing
*/
scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
}
}
else
{
/* curr to the left of prev
*/
dx = - dx;
if (dx > ((2*M_PI) - dx))
{
/* vertical edge crossing to the right of prev
*/
/* find exit point (right edge) */
mx = M_PI;
my = merc_find_edge_xing(prev, curr);
/* scan from prev to exit point */
scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
/* (mx, my) is an edge crossing (exit point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeExit;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 0; /* right edge */
/* scan from entry point (left edge) to curr */
mx = - M_PI;
scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
/* (mx, my) is an edge crossing (entry point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeEntry;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 2; /* left edge */
}
else
{
/* no vertical edge crossing
*/
scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
}
}
}
double merc_find_edge_xing(prev, curr)
double *prev;
double *curr;
{
double ratio;
double scale;
double z1, z2;
double rslt;
if (curr[0] != 0)
{
ratio = (prev[0] / curr[0]);
z1 = prev[1] - (ratio * curr[1]);
z2 = prev[2] - (ratio * curr[2]);
}
else
{
z1 = curr[1];
z2 = curr[2];
}
scale = ((z2 > 0) ? -1 : 1) / sqrt((z1*z1) + (z2*z2));
z1 *= scale;
rslt = MERCATOR_Y(z1);
return rslt;
}
void merc_handle_xings()
{
int i;
int nxings;
EdgeXing *xings;
EdgeXing *from;
EdgeXing *to;
xings = (EdgeXing *) edgexings->body;
nxings = edgexings->count;
assert((nxings % 2) == 0);
qsort(xings, (unsigned) nxings, sizeof(EdgeXing), merc_edgexing_comp);
if (xings[0].type == XingTypeExit)
{
for (i=0; i<nxings; i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
merc_scan_edge(from, to);
}
}
else
{
from = &(xings[nxings-1]);
to = &(xings[0]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry) ||
(from->angle < to->angle))
xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
merc_scan_edge(from, to);
for (i=1; i<(nxings-1); i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
merc_scan_edge(from, to);
}
}
edgexings->count = 0;
}
void merc_scan_edge(from, to)
EdgeXing *from;
EdgeXing *to;
{
int s0, s1, s_new;
double x_0, x_1, x_new;
double y_0, y_1, y_new;
s0 = from->angle;
x_0 = XPROJECT(from->x);
y_0 = YPROJECT(from->y);
s1 = to->angle;
x_1 = XPROJECT(to->x);
y_1 = YPROJECT(to->y);
while (s0 != s1)
{
switch (s0)
{
case 0:
x_new = XPROJECT(M_PI);
y_new = YPROJECT(BigNumber);
s_new = 1;
break;
case 1:
x_new = XPROJECT(-M_PI);
y_new = YPROJECT(BigNumber);
s_new = 2;
break;
case 2:
x_new = XPROJECT(-M_PI);
y_new = YPROJECT(-BigNumber);
s_new = 3;
break;
case 3:
x_new = XPROJECT(M_PI);
y_new = YPROJECT(-BigNumber);
s_new = 0;
break;
default:
/* keep lint happy */
x_new = 0;
y_new = 0;
s_new = 0;
assert(0);
}
scan(x_0, y_0, x_new, y_new);
x_0 = x_new;
y_0 = y_new;
s0 = s_new;
}
scan(x_0, y_0, x_1, y_1);
}
void cyl_scan_outline()
{
double left, right;
double top, bottom;
min_y = hght;
max_y = -1;
left = XPROJECT(-M_PI);
right = XPROJECT(M_PI);
top = YPROJECT(BigNumber);
bottom = YPROJECT(-BigNumber);
scan(right, top, left, top);
scan(left, top, left, bottom);
scan(left, bottom, right, bottom);
scan(right, bottom, right, top);
get_scanbits(64);
}
void cyl_scan_curves()
{
int i;
int cidx;
int npts;
int val;
short *raw;
double *pos;
double *prev;
double *curr;
cidx = 0;
raw = map_data;
while (1)
{
npts = raw[0];
if (npts == 0) break;
val = raw[1];
raw += 2;
pos = cyl_extract_curve(npts, raw);
prev = pos + (npts-1)*5;
curr = pos;
min_y = hght;
max_y = -1;
for (i=0; i<npts; i++)
{
cyl_scan_along_curve(prev, curr, cidx);
prev = curr;
curr += 5;
}
free(pos);
if (edgexings->count > 0)
cyl_handle_xings();
if (min_y <= max_y)
get_scanbits(val);
cidx += 1;
raw += 3*npts;
}
}
double *cyl_extract_curve(npts, data)
int npts;
short *data;
{
int i;
int x, y, z;
double scale;
double *pos;
double *rslt;
rslt = (double *) malloc((unsigned) sizeof(double) * 5 * npts);
assert(rslt != NULL);
x = 0;
y = 0;
z = 0;
scale = 1.0 / MAP_DATA_SCALE;
pos = rslt;
for (i=0; i<npts; i++)
{
x += data[0];
y += data[1];
z += data[2];
pos[0] = x * scale;
pos[1] = y * scale;
pos[2] = z * scale;
XFORM_ROTATE(pos, view_pos_info);
/* apply cylindrical projection
*/
pos[3] = CYLINDRICAL_X(pos[0], pos[2]);
pos[4] = CYLINDRICAL_Y(pos[1]);
data += 3;
pos += 5;
}
return rslt;
}
void cyl_scan_along_curve(prev, curr, cidx)
double *prev;
double *curr;
int cidx;
{
double px, py;
double cx, cy;
double dx;
double mx, my;
EdgeXing *xing;
px = prev[3];
cx = curr[3];
py = prev[4];
cy = curr[4];
dx = cx - px;
if (dx > 0)
{
/* curr to the right of prev
*/
if (dx > ((2*M_PI) - dx))
{
/* vertical edge crossing to the left of prev
*/
/* find exit point (left edge) */
mx = - M_PI;
my = cyl_find_edge_xing(prev, curr);
/* scan from prev to exit point */
scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
/* (mx, my) is an edge crossing (exit point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeExit;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 2; /* left edge */
/* scan from entry point (right edge) to curr */
mx = M_PI;
scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
/* (mx, my) is an edge crossing (entry point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeEntry;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 0; /* right edge */
}
else
{
/* no vertical edge crossing
*/
scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
}
}
else
{
/* curr to the left of prev
*/
dx = - dx;
if (dx > ((2*M_PI) - dx))
{
/* vertical edge crossing to the right of prev
*/
/* find exit point (right edge) */
mx = M_PI;
my = cyl_find_edge_xing(prev, curr);
/* scan from prev to exit point */
scan(XPROJECT(px), YPROJECT(py), XPROJECT(mx), YPROJECT(my));
/* (mx, my) is an edge crossing (exit point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeExit;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 0; /* right edge */
/* scan from entry point (left edge) to curr */
mx = - M_PI;
scan(XPROJECT(mx), YPROJECT(my), XPROJECT(cx), YPROJECT(cy));
/* (mx, my) is an edge crossing (entry point) */
xing = (EdgeXing *) extarr_next(edgexings);
xing->type = XingTypeEntry;
xing->cidx = cidx;
xing->x = mx;
xing->y = my;
xing->angle = 2; /* left edge */
}
else
{
/* no vertical edge crossing
*/
scan(XPROJECT(px), YPROJECT(py), XPROJECT(cx), YPROJECT(cy));
}
}
}
double cyl_find_edge_xing(prev, curr)
double *prev;
double *curr;
{
double ratio;
double scale;
double z1, z2;
double rslt;
if (curr[0] != 0)
{
ratio = (prev[0] / curr[0]);
z1 = prev[1] - (ratio * curr[1]);
z2 = prev[2] - (ratio * curr[2]);
}
else
{
z1 = curr[1];
z2 = curr[2];
}
scale = ((z2 > 0) ? -1 : 1) / sqrt((z1*z1) + (z2*z2));
z1 *= scale;
rslt = CYLINDRICAL_Y(z1);
return rslt;
}
void cyl_handle_xings()
{
int i;
int nxings;
EdgeXing *xings;
EdgeXing *from;
EdgeXing *to;
xings = (EdgeXing *) edgexings->body;
nxings = edgexings->count;
assert((nxings % 2) == 0);
qsort(xings, (unsigned) nxings, sizeof(EdgeXing), cyl_edgexing_comp);
if (xings[0].type == XingTypeExit)
{
for (i=0; i<nxings; i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
cyl_scan_edge(from, to);
}
}
else
{
from = &(xings[nxings-1]);
to = &(xings[0]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry) ||
(from->angle < to->angle))
xing_error(__FILE__, __LINE__, nxings-1, nxings, xings);
cyl_scan_edge(from, to);
for (i=1; i<(nxings-1); i+=2)
{
from = &(xings[i]);
to = &(xings[i+1]);
if ((from->type != XingTypeExit) ||
(to->type != XingTypeEntry))
xing_error(__FILE__, __LINE__, i, nxings, xings);
cyl_scan_edge(from, to);
}
}
edgexings->count = 0;
}
void cyl_scan_edge(from, to)
EdgeXing *from;
EdgeXing *to;
{
int s0, s1, s_new;
double x_0, x_1, x_new;
double y_0, y_1, y_new;
s0 = from->angle;
x_0 = XPROJECT(from->x);
y_0 = YPROJECT(from->y);
s1 = to->angle;
x_1 = XPROJECT(to->x);
y_1 = YPROJECT(to->y);
while (s0 != s1)
{
switch (s0)
{
case 0:
x_new = XPROJECT(M_PI);
y_new = YPROJECT(BigNumber);
s_new = 1;
break;
case 1:
x_new = XPROJECT(-M_PI);
y_new = YPROJECT(BigNumber);
s_new = 2;
break;
case 2:
x_new = XPROJECT(-M_PI);
y_new = YPROJECT(-BigNumber);
s_new = 3;
break;
case 3:
x_new = XPROJECT(M_PI);
y_new = YPROJECT(-BigNumber);
s_new = 0;
break;
default:
/* keep lint happy */
x_new = 0;
y_new = 0;
s_new = 0;
assert(0);
}
scan(x_0, y_0, x_new, y_new);
x_0 = x_new;
y_0 = y_new;
s0 = s_new;
}
scan(x_0, y_0, x_1, y_1);
}
void xing_error(file, line, idx, nxings, xings)
const char *file;
int line;
int idx;
int nxings;
EdgeXing *xings;
{
fflush(stdout);
fprintf(stderr, "xearth %s: incorrect edgexing sequence (%s:%d)\n",
VersionString, file, line);
fprintf(stderr, " (cidx %d) (view_lat %.16f) (view_lon %.16f)\n",
xings[idx].cidx, view_lat, view_lon);
fprintf(stderr, "\n");
exit(1);
}
void scan(x_0, y_0, x_1, y_1)
double x_0, y_0;
double x_1, y_1;
{
int i;
int lo_y, hi_y;
double x_value;
double x_delta;
if (y_0 < y_1)
{
lo_y = ceil(y_0 - 0.5);
hi_y = floor(y_1 - 0.5);
if (hi_y == (y_1 - 0.5))
hi_y -= 1;
}
else
{
lo_y = ceil(y_1 - 0.5);
hi_y = floor(y_0 - 0.5);
if (hi_y == (y_0 - 0.5))
hi_y -= 1;
}
if (lo_y < 0) lo_y = 0;
if (hi_y >= hght) hi_y = hght-1;
if (lo_y > hi_y)
return; /* no scan lines crossed */
if (lo_y < min_y) min_y = lo_y;
if (hi_y > max_y) max_y = hi_y;
x_delta = (x_1 - x_0) / (y_1 - y_0);
x_value = x_0 + x_delta * ((lo_y + 0.5) - y_0);
for (i=lo_y; i<=hi_y; i++)
{
*((double *) extarr_next(scanbuf[i])) = x_value;
x_value += x_delta;
}
}
void get_scanbits(val)
int val;
{
int i, j;
int lo_x, hi_x;
unsigned nvals;
double *vals;
ScanBit *scanbit;
for (i=min_y; i<=max_y; i++)
{
vals = (double *) scanbuf[i]->body;
nvals = scanbuf[i]->count;
assert((nvals % 2) == 0);
qsort(vals, nvals, sizeof(double), double_comp);
for (j=0; j<nvals; j+=2)
{
lo_x = ceil(vals[j] - 0.5);
hi_x = floor(vals[j+1] - 0.5);
if (lo_x < 0) lo_x = 0;
if (hi_x >= wdth) hi_x = wdth-1;
if (lo_x <= hi_x)
{
scanbit = (ScanBit *) extarr_next(scanbits);
scanbit->y = i;
scanbit->lo_x = lo_x;
scanbit->hi_x = hi_x;
scanbit->val = val;
}
}
scanbuf[i]->count = 0;
}
}
syntax highlighted by Code2HTML, v. 0.9.1