1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
|
#ifndef QUADTILE_H
#define QUADTILE_H
typedef signed long long quadtile;
// Unseralize a quadtile from 8 bytes
quadtile buf2ll(const unsigned char *s)
{
quadtile result=0;
for(int i=0; i<8;i++)
result = (result << 8LL) | s[i];
return result;
}
// Unserialize an unsigned long from 4 bytes
unsigned long buf2l(const unsigned char *s)
{
unsigned long result=0;
for(int i=0; i<4;i++)
result = (result << 8UL) | s[i];
return result;
}
// Take two (<=4bytes) numbers, and alternate their binary digits.
quadtile mux(quadtile x, quadtile y)
{
long long t=0;
for(int i=0;i<32;i++) {
t <<= 1;
if(x & (1ULL<<31)) t |= 1;
t <<= 1;
if(y & (1ULL<<31)) t |= 1;
x <<= 1; y <<= 1;
}
return t;
}
void demux(quadtile tile, quadtile *x, quadtile *y)
{
*x = 0; *y = 0;
for(int i=0;i<31;i++) {
//binary_printf(tile); binary_printf(*x); binary_printf(*y);printf("\n");
(*y) <<= 1; (*x) <<= 1;
if(tile & (1ULL<<61)) (*x)|=1;
tile <<= 1;
if(tile & (1ULL<<61)) (*y)|=1;
tile <<= 1;
}
//binary_printf(tile); binary_printf(*x); binary_printf(*y);printf("\n");
}
/* Store a pair of doubles in the range (0.0->1.0) as integers with alternating
binary bits */
quadtile xy2q(double fx, double fy)
{
return mux((long long) (fx * (1ULL<<31)), (long long) (fy * (1ULL<<31)));
}
#ifdef NEED_QTILE_WRITE
// Below here is only needed in the preprocessor.
#include <math.h>
/* latlon to projected x/y coords */
void ll2pxy(double lat, double lon, unsigned long *x, unsigned long *y)
{
double fx,fy;
fx = (lon+180.0) / 360.0;
fy = (1.0 - log(tan(lat*M_PI/180.0) +
(1.0/cos(lat*M_PI/180.0)))/M_PI)/2.0;
*x = (unsigned long) (fx * (1UL<<31));
*y = (unsigned long) (fy * (1UL<<31));
}
/* Serialise a quadtile to 4 unsigned char bytes, MSB first */
unsigned char *ll2buf(quadtile q)
{
static unsigned char result[8];
for(int i=0; i<8;i++) {
result[7-i] = q & 255;
q >>= 8;
}
return result;
}
unsigned char *l2buf(unsigned long l)
{
static unsigned char result[4];
for(int i=0; i<4;i++) {
result[3-i] = l & 255;
l >>= 8;
}
return result;
}
/* Horrible name! Assume a line drawn between q1 and q2.
if (q1 & qmask) != (q2 & qmask) then return the quadtile q corresponding
to the first point along the line q1->q2 where (q & qmask != q1 & qmask) */
quadtile line_edge_intersect(quadtile q1, quadtile q2, quadtile qmask)
{
quadtile qmin, qmax, xmin, xmax, ymin, ymax, x1, y1, x2, y2;
if((q1 & qmask) == (q2 & qmask)) {
printf("Error - q1 and q2 are in the same qtile\n");
exit(-1);
}
//Get the bounds of q1's quadtile.
qmin = q1 & qmask;
qmax = q1 | (~qmask);
demux(qmin, &xmin, &ymin);
demux(qmax, &xmax, &ymax);
demux(q1, &x1, &y1);
demux(q2, &x2, &y2);
//y = ax + b x = (y-b)/a
//First a special case for x1==x2 (a==infinity. Doesn't work.)
if(x1==x2) {
if(y2>y1) return mux(x1, ymax+1);
else return mux(x1, ymin-1);
}
double a = ((double)(y2-y1)) / ((double)(x2-x1));
double b = y1 - a *x1;
//Find where the line crosses the bounds of this quadtile.
//FIXME this will fail at the international date line - do roads cross this?
if(x2 > x1) {
quadtile y = roundl(a*(xmax) + b);
if(y>=ymin && y<=ymax) return mux(xmax+1, roundl(a*(xmax+1) + b));
} else {
quadtile y = roundl(a*(xmin) + b);
if(y>=ymin && y<=ymax) return mux(xmin-1, roundl(a*(xmin-1) + b));
}
if(y2 > y1) {
quadtile x = roundl((ymax-b)/a);
if(x>=xmin && x<=xmax) return mux(roundl((ymax+1-b)/a), ymax+1);
} else {
quadtile x = roundl((ymin-b)/a);
if(x>=xmin && x<=xmax) return mux(roundl((ymin-b-1)/a), ymin-1);
}
printf("Error - shouldn't get here in intersect calcs\n");
//Debugging output.
printf("q1=%llx q2=%llx qmask=%llx\na=%f, b=%f\n1=%lld,%lld 2=%lld,%lld\n",
q1, q2, qmask, a, b, x1, y1, x2, y2);
printf("xmin = %lld xmax = %lld, ymin = %lld, ymax=%lld\n",
xmin, xmax, ymin, ymax);
printf("Maxy = %f, Miny=%f\n", a*(xmax+1)+b, a*(ymin-1)+b);
printf("Maxx = %f, Minx=%f\n", (ymax+1-b)/a, (ymin-1-b)/a);
exit(-1);
return 0;
}
#endif /*NEED_QTILE_WRITE*/
#endif /*QUADTILE_H*/
|
