diegoripley/tippecanoe
1
0
Fork 0
mirror of https://github.com/felt/tippecanoe.git synced 2026-10-02 16:35:40 +02:00
Code Issues Packages Projects Releases Wiki Activity

Use 128-bit arithmetic to avoid overflow

Browse Source
This commit is contained in:
Erica Fischer
2024-03-18 11:37:47 -07:00
parent cc6135dd02
commit 3b948e0f20
4 changed files with 24 additions and 36 deletions
Download Patch File Download Diff File Expand all files Collapse all files
clip.cpp
+5 -7
View File
@@ -177,28 +177,26 @@ int clip(long long *x0, long long *y0, long long *x1, long long *y1, long long x
} else {
// failed both tests, so calculate the line segment to clip
// from an outside point to an intersection with clip edge
long long x = *x0, y = *y0;
__int128 x = *x0, y = *y0;
// At least one endpoint is outside the clip rectangle; pick it.
int outcodeOut = outcode0 ? outcode0 : outcode1;
// XXX truncating division
long long shift = 1LL << (GLOBAL_DETAIL - 32);
// Now find the intersection point;
// use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0)
if (outcodeOut & TOP) { // point is above the clip rectangle
x = *x0 + ((*x1 - *x0) / shift) * ((ymax - *y0) / shift) / (*y1 - *y0) * shift * shift;
x = *x0 + ((__int128) (*x1 - *x0)) * ((ymax - *y0)) / (*y1 - *y0);
y = ymax;
} else if (outcodeOut & BOTTOM) { // point is below the clip rectangle
x = *x0 + ((*x1 - *x0) / shift) * ((ymin - *y0) / shift) / (*y1 - *y0) * shift * shift;
x = *x0 + ((__int128) (*x1 - *x0)) * ((ymin - *y0)) / (*y1 - *y0);
y = ymin;
} else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle
y = *y0 + ((*y1 - *y0) / shift) * ((xmax - *x0) / shift) / (*x1 - *x0) * shift * shift;
y = *y0 + ((__int128) (*y1 - *y0)) * ((xmax - *x0)) / (*x1 - *x0);
x = xmax;
} else if (outcodeOut & LEFT) { // point is to the left of clip rectangle
y = *y0 + ((*y1 - *y0) / shift) * ((xmin - *x0) / shift) / (*x1 - *x0) * shift * shift;
y = *y0 + ((__int128) (*y1 - *y0)) * ((xmin - *x0)) / (*x1 - *x0);
x = xmin;
}
geometry.cpp
+17 -27
View File
@@ -300,21 +300,15 @@ bool point_within_tile(long long x, long long y, int z) {
return x >= 0 && y >= 0 && x < area && y < area;
}
double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) {
long long p2x = segB_x - segA_x;
long long p2y = segB_y - segA_y;
double distance_from_line(__int128 point_x, __int128 point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) {
__int128 p2x = segB_x - segA_x;
__int128 p2y = segB_y - segA_y;
// These calculations must be made in integers instead of floating point
// to make them consistent between x86 and arm floating point implementations.
//
// In a 32-bit world, coordinates may be up to 34 bits, so their product is up to 68 bits,
// making their sum up to 69 bits. Downshift before multiplying to keep them in range.
//
// If the world is bigger than 32 bits, scale down to 32 bits.
long long shift = 1LL << (GLOBAL_DETAIL - 32);
double something = ((p2x / 4 / shift) * (p2x / 8 / shift) + (p2y / 4 / shift) * (p2y / 8 / shift)) * 32.0 * shift * shift;
double something = (p2x) * (p2x) + (p2y) * (p2y);
// likewise
double u = (0 == something) ? 0 : ((point_x - segA_x) / 4 / shift * (p2x / 8 / shift) + (point_y - segA_y) / 4 / shift * (p2y / 8 / shift)) * 32.0 * shift * shift / (something);
double u = (0 == something) ? 0 : ((point_x - segA_x) * (p2x) + (point_y - segA_y) * (p2y)) / (something);
if (u >= 1) {
u = 1;
@@ -674,8 +668,8 @@ drawvec fix_polygon(const drawvec &geom) {
// calculate centroid
// a + 1 < size() because point 0 is duplicated at the end
long long xtotal = 0;
long long ytotal = 0;
__int128 xtotal = 0;
__int128 ytotal = 0;
long long count = 0;
for (size_t a = 0; a + 1 < ring.size(); a++) {
xtotal += ring[a].x;
@@ -685,16 +679,13 @@ drawvec fix_polygon(const drawvec &geom) {
xtotal /= count;
ytotal /= count;
long long shift = 1LL << (GLOBAL_DETAIL - 32);
// figure out which point is furthest from the centroid
long long dist2 = 0;
long long furthest = 0;
__int128 dist2 = 0;
size_t furthest = 0;
for (size_t a = 0; a + 1 < ring.size(); a++) {
// division by 16 because these are z0 coordinates and we need to avoid overflow
long long xd = (ring[a].x - xtotal) / 16 / shift;
long long yd = (ring[a].y - ytotal) / 16 / shift;
long long d2 = xd * xd + yd * yd;
__int128 xd = (ring[a].x - xtotal);
__int128 yd = (ring[a].y - ytotal);
__int128 d2 = xd * xd + yd * yd;
if (d2 > dist2 || (d2 == dist2 && ring[a] < ring[furthest])) {
dist2 = d2;
furthest = a;
@@ -704,13 +695,12 @@ drawvec fix_polygon(const drawvec &geom) {
// then figure out which point is furthest from *that*,
// which will hopefully be a good origin point since it should be
// at a far edge of the shape.
long long dist2b = 0;
long long furthestb = 0;
__int128 dist2b = 0;
size_t furthestb = 0;
for (size_t a = 0; a + 1 < ring.size(); a++) {
// division by 16 because these are z0 coordinates and we need to avoid overflow
long long xd = (ring[a].x - ring[furthest].x) / 16 / shift;
long long yd = (ring[a].y - ring[furthest].y) / 16 / shift;
long long d2 = xd * xd + yd * yd;
__int128 xd = (ring[a].x - ring[furthest].x);
__int128 yd = (ring[a].y - ring[furthest].y);
__int128 d2 = xd * xd + yd * yd;
if (d2 > dist2b || (d2 == dist2b && ring[a] < ring[furthestb])) {
dist2b = d2;
furthestb = a;
geometry.hpp
+1 -1
View File
@@ -98,7 +98,7 @@ drawvec clip_lines(drawvec &geom, long long x1, long long y1, long long x2, long
drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long long y2);
void visvalingam(drawvec &ls, size_t start, size_t end, double threshold, size_t retain);
int pnpoly(const drawvec &vert, size_t start, size_t nvert, long long testx, long long testy);
double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y);
double distance_from_line(__int128 point_x, __int128 point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y);
std::string overzoom(const mvt_tile &tile, int oz, int ox, int oy, int nz, int nx, int ny,
int detail, int buffer, std::set<std::string> const &keep, bool do_compress,
main.cpp
+1 -1
View File
@@ -2431,7 +2431,7 @@ std::pair<int, metadata> read_input(std::vector<source> &sources, char *fname, i
double total_tile_count = 0;
for (int i = 1; i <= maxzoom; i++) {
double tile_count = ceil(area_sum / ((1LL << (GLOBAL_DETAIL - i)) * (1LL << (GLOBAL_DETAIL - i))));
double tile_count = ceil(area_sum / ((__int128) (1LL << (GLOBAL_DETAIL - i)) * (1LL << (GLOBAL_DETAIL - i))));
total_tile_count += tile_count;
// 2M tiles is an arbitrary limit, chosen to make tiling jobs