1 /********************************************************************** 2 * 3 * PostGIS - Spatial Types for PostgreSQL 4 * http://postgis.net 5 * 6 * PostGIS is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 2 of the License, or 9 * (at your option) any later version. 10 * 11 * PostGIS is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with PostGIS. If not, see <http://www.gnu.org/licenses/>. 18 * 19 ********************************************************************** 20 * 21 * Copyright (C) 2011-2021 Sandro Santilli <strk@kbt.io> 22 * Copyright (C) 2011 Paul Ramsey <pramsey@cleverelephant.ca> 23 * Copyright (C) 2007-2008 Mark Cave-Ayland 24 * Copyright (C) 2001-2006 Refractions Research Inc. 25 * 26 **********************************************************************/ 27 28 29 #ifndef _LIBLWGEOM_INTERNAL_H 30 #define _LIBLWGEOM_INTERNAL_H 1 31 32 #include "../postgis_config.h" 33 34 #include "lwgeom_log.h" 35 36 #include <assert.h> 37 #include <stdarg.h> 38 #include <stdint.h> 39 #include <stdio.h> 40 #include <string.h> 41 #include <stdlib.h> 42 #include <float.h> 43 #include <math.h> 44 45 #if HAVE_IEEEFP_H 46 #include <ieeefp.h> 47 #endif 48 49 #include "liblwgeom.h" 50 51 /** 52 * Floating point comparators. 53 */ 54 #define FP_TOLERANCE 1e-12 55 #define FP_IS_ZERO(A) (fabs(A) <= FP_TOLERANCE) 56 #define FP_MAX(A, B) (((A) > (B)) ? (A) : (B)) 57 #define FP_MIN(A, B) (((A) < (B)) ? (A) : (B)) 58 #define FP_ABS(a) ((a) < (0) ? -(a) : (a)) 59 #define FP_EQUALS(A, B) (fabs((A)-(B)) <= FP_TOLERANCE) 60 #define FP_NEQUALS(A, B) (fabs((A)-(B)) > FP_TOLERANCE) 61 #define FP_LT(A, B) (((A) + FP_TOLERANCE) < (B)) 62 #define FP_LTEQ(A, B) (((A) - FP_TOLERANCE) <= (B)) 63 #define FP_GT(A, B) (((A) - FP_TOLERANCE) > (B)) 64 #define FP_GTEQ(A, B) (((A) + FP_TOLERANCE) >= (B)) 65 #define FP_CONTAINS_TOP(A, X, B) (FP_LT(A, X) && FP_LTEQ(X, B)) 66 #define FP_CONTAINS_BOTTOM(A, X, B) (FP_LTEQ(A, X) && FP_LT(X, B)) 67 #define FP_CONTAINS_INCL(A, X, B) (FP_LTEQ(A, X) && FP_LTEQ(X, B)) 68 #define FP_CONTAINS_EXCL(A, X, B) (FP_LT(A, X) && FP_LT(X, B)) 69 #define FP_CONTAINS(A, X, B) FP_CONTAINS_EXCL(A, X, B) 70 71 #define STR_EQUALS(A, B) strcmp((A), (B)) == 0 72 #define STR_IEQUALS(A, B) (strcasecmp((A), (B)) == 0) 73 #define STR_ISTARTS(A, B) (strncasecmp((A), (B), strlen((B))) == 0) 74 75 76 /* 77 * this will change to NaN when I figure out how to 78 * get NaN in a platform-independent way 79 */ 80 #define NO_VALUE 0.0 81 #define NO_Z_VALUE NO_VALUE 82 #define NO_M_VALUE NO_VALUE 83 84 85 /** 86 * Well-Known Text (WKT) Output Variant Types 87 */ 88 #define WKT_NO_TYPE 0x08 /* Internal use only */ 89 #define WKT_NO_PARENS 0x10 /* Internal use only */ 90 #define WKT_IS_CHILD 0x20 /* Internal use only */ 91 92 /** 93 * Well-Known Binary (WKB) Output Variant Types 94 */ 95 96 #define WKB_DOUBLE_SIZE 8 /* Internal use only */ 97 #define WKB_INT_SIZE 4 /* Internal use only */ 98 #define WKB_BYTE_SIZE 1 /* Internal use only */ 99 100 /** 101 * Well-Known Binary (WKB) Geometry Types 102 */ 103 #define WKB_POINT_TYPE 1 104 #define WKB_LINESTRING_TYPE 2 105 #define WKB_POLYGON_TYPE 3 106 #define WKB_MULTIPOINT_TYPE 4 107 #define WKB_MULTILINESTRING_TYPE 5 108 #define WKB_MULTIPOLYGON_TYPE 6 109 #define WKB_GEOMETRYCOLLECTION_TYPE 7 110 #define WKB_CIRCULARSTRING_TYPE 8 111 #define WKB_COMPOUNDCURVE_TYPE 9 112 #define WKB_CURVEPOLYGON_TYPE 10 113 #define WKB_MULTICURVE_TYPE 11 114 #define WKB_MULTISURFACE_TYPE 12 115 #define WKB_CURVE_TYPE 13 /* from ISO draft, not sure is real */ 116 #define WKB_SURFACE_TYPE 14 /* from ISO draft, not sure is real */ 117 #define WKB_POLYHEDRALSURFACE_TYPE 15 118 #define WKB_TIN_TYPE 16 119 #define WKB_TRIANGLE_TYPE 17 120 121 122 /** 123 * Macro that returns: 124 * -1 if n < 0, 125 * 1 if n > 0, 126 * 0 if n == 0 127 */ 128 #define SIGNUM(n) (((n) > 0) - ((n) < 0)) 129 130 /** 131 * Tolerance used to determine equality. 132 */ 133 #define EPSILON_SQLMM 1e-8 134 135 /* 136 * Export functions 137 */ 138 139 /* Any (absolute) values outside this range will be printed in scientific notation */ 140 #define OUT_MIN_DOUBLE 1E-8 141 #define OUT_MAX_DOUBLE 1E15 142 #define OUT_DEFAULT_DECIMAL_DIGITS 15 143 144 /* 17 digits are sufficient for round-tripping 145 * Then we might add up to 8 (from OUT_MIN_DOUBLE) max leading zeroes (or 2 digits for "e+") */ 146 #define OUT_MAX_DIGITS 17 + 8 147 148 /* Limit for the max amount of characters that a double can use, including dot and sign */ 149 /* */ 150 #define OUT_MAX_BYTES_DOUBLE (1 /* Sign */ + 2 /* 0.x */ + OUT_MAX_DIGITS) 151 #define OUT_DOUBLE_BUFFER_SIZE OUT_MAX_BYTES_DOUBLE + 1 /* +1 including NULL */ 152 153 /** 154 * Constants for point-in-polygon return values 155 */ 156 #define LW_INSIDE 1 157 #define LW_BOUNDARY 0 158 #define LW_OUTSIDE -1 159 160 /* 161 * Internal prototypes 162 */ 163 164 /* Machine endianness */ 165 #define XDR 0 /* big endian */ 166 #define NDR 1 /* little endian */ 167 168 169 /* 170 * Force dims 171 */ 172 LWGEOM* lwgeom_force_dims(const LWGEOM *lwgeom, int hasz, int hasm, double zval, double mval); 173 LWPOINT* lwpoint_force_dims(const LWPOINT *lwpoint, int hasz, int hasm, double zval, double mval); 174 LWLINE* lwline_force_dims(const LWLINE *lwline, int hasz, int hasm, double zval, double mval); 175 LWPOLY* lwpoly_force_dims(const LWPOLY *lwpoly, int hasz, int hasm, double zval, double mval); 176 LWCOLLECTION* lwcollection_force_dims(const LWCOLLECTION *lwcol, int hasz, int hasm, double zval, double mval); 177 POINTARRAY* ptarray_force_dims(const POINTARRAY *pa, int hasz, int hasm, double zval, double mval); 178 179 /** 180 * Swap ordinate values o1 and o2 on a given POINTARRAY 181 * 182 * Ordinates semantic is: 0=x 1=y 2=z 3=m 183 */ 184 void ptarray_swap_ordinates(POINTARRAY *pa, LWORD o1, LWORD o2); 185 186 /* 187 * Is Empty? 188 */ 189 int lwpoly_is_empty(const LWPOLY *poly); 190 int lwcollection_is_empty(const LWCOLLECTION *col); 191 int lwcircstring_is_empty(const LWCIRCSTRING *circ); 192 int lwtriangle_is_empty(const LWTRIANGLE *triangle); 193 int lwline_is_empty(const LWLINE *line); 194 int lwpoint_is_empty(const LWPOINT *point); 195 196 /* 197 * Number of vertices? 198 */ 199 uint32_t lwline_count_vertices(LWLINE *line); 200 uint32_t lwpoly_count_vertices(LWPOLY *poly); 201 uint32_t lwcollection_count_vertices(LWCOLLECTION *col); 202 203 /* 204 * DP simplification 205 */ 206 207 /** 208 * @param minpts minimum number of points to retain, if possible. 209 */ 210 void ptarray_simplify_in_place(POINTARRAY *pa, double tolerance, uint32_t minpts); 211 212 /* 213 * The possible ways a pair of segments can interact. Returned by lw_segment_intersects 214 */ 215 enum CG_SEGMENT_INTERSECTION_TYPE { 216 SEG_ERROR = -1, 217 SEG_NO_INTERSECTION = 0, 218 SEG_COLINEAR = 1, 219 SEG_CROSS_LEFT = 2, 220 SEG_CROSS_RIGHT = 3, 221 SEG_TOUCH_LEFT = 4, 222 SEG_TOUCH_RIGHT = 5 223 }; 224 225 /* 226 * Do the segments intersect? How? 227 */ 228 int lw_segment_intersects(const POINT2D *p1, const POINT2D *p2, const POINT2D *q1, const POINT2D *q2); 229 230 /* 231 * Get/Set an enumeratoed ordinate. (x,y,z,m) 232 */ 233 double lwpoint_get_ordinate(const POINT4D *p, char ordinate); 234 void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value); 235 236 /* 237 * Generate an interpolated coordinate p given an interpolation value and ordinate to apply it to 238 */ 239 int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value); 240 241 242 /* 243 * Geohash 244 */ 245 int lwgeom_geohash_precision(GBOX bbox, GBOX *bounds); 246 lwvarlena_t *geohash_point(double longitude, double latitude, int precision); 247 void decode_geohash_bbox(char *geohash, double *lat, double *lon, int precision); 248 249 /* 250 * Point comparisons 251 */ 252 int p4d_same(const POINT4D *p1, const POINT4D *p2); 253 int p3d_same(const POINT3D *p1, const POINT3D *p2); 254 int p2d_same(const POINT2D *p1, const POINT2D *p2); 255 256 /* 257 * Area calculations 258 */ 259 double lwpoly_area(const LWPOLY *poly); 260 double lwcurvepoly_area(const LWCURVEPOLY *curvepoly); 261 double lwtriangle_area(const LWTRIANGLE *triangle); 262 263 /** 264 * Pull a #GBOX from the header of a #GSERIALIZED, if one is available. If 265 * it is not, return LW_FAILURE. 266 */ 267 int gserialized_read_gbox_p(const GSERIALIZED *g, GBOX *gbox); 268 269 /* 270 * Populate a bounding box *without* allocating an LWGEOM. Useful for some performance 271 * purposes. Use only if gserialized_read_gbox_p failed 272 */ 273 int gserialized_peek_gbox_p(const GSERIALIZED *g, GBOX *gbox); 274 275 /** 276 * Calculate required memory segment to contain a serialized form of the LWGEOM. 277 * Primarily used internally by the serialization code. Exposed to allow the cunit 278 * tests to exercise it. 279 */ 280 size_t gserialized_from_lwgeom_size(const LWGEOM *geom); 281 282 /* 283 * Length calculations 284 */ 285 double lwcompound_length(const LWCOMPOUND *comp); 286 double lwcompound_length_2d(const LWCOMPOUND *comp); 287 double lwline_length(const LWLINE *line); 288 double lwline_length_2d(const LWLINE *line); 289 double lwcircstring_length(const LWCIRCSTRING *circ); 290 double lwcircstring_length_2d(const LWCIRCSTRING *circ); 291 double lwpoly_perimeter(const LWPOLY *poly); 292 double lwpoly_perimeter_2d(const LWPOLY *poly); 293 double lwcurvepoly_perimeter(const LWCURVEPOLY *poly); 294 double lwcurvepoly_perimeter_2d(const LWCURVEPOLY *poly); 295 double lwtriangle_perimeter(const LWTRIANGLE *triangle); 296 double lwtriangle_perimeter_2d(const LWTRIANGLE *triangle); 297 298 /* 299 * Segmentization 300 */ 301 LWPOLY *lwcurvepoly_stroke(const LWCURVEPOLY *curvepoly, uint32_t perQuad); 302 303 /* 304 * Affine 305 */ 306 void ptarray_affine(POINTARRAY *pa, const AFFINE *affine); 307 void affine_invert(AFFINE *affine); 308 309 /* 310 * Scale 311 */ 312 void ptarray_scale(POINTARRAY *pa, const POINT4D *factor); 313 314 /* 315 * Scroll 316 */ 317 int ptarray_scroll_in_place(POINTARRAY *pa, const POINT4D *newbase); 318 319 /* 320 * PointArray 321 */ 322 int ptarray_has_z(const POINTARRAY *pa); 323 int ptarray_has_m(const POINTARRAY *pa); 324 double ptarray_signed_area(const POINTARRAY *pa); 325 326 /* 327 * Length 328 */ 329 double ptarray_length(const POINTARRAY *pts); 330 double ptarray_arc_length_2d(const POINTARRAY *pts); 331 332 /* 333 * Clone support 334 */ 335 LWPOINT *lwpoint_clone(const LWPOINT *lwgeom); 336 LWLINE *lwline_clone(const LWLINE *lwgeom); 337 LWPOLY *lwpoly_clone(const LWPOLY *lwgeom); 338 LWTRIANGLE *lwtriangle_clone(const LWTRIANGLE *lwgeom); 339 LWCOLLECTION *lwcollection_clone(const LWCOLLECTION *lwgeom); 340 LWCIRCSTRING *lwcircstring_clone(const LWCIRCSTRING *curve); 341 POINTARRAY *ptarray_clone(const POINTARRAY *ptarray); 342 LWLINE *lwline_clone_deep(const LWLINE *lwgeom); 343 LWPOLY *lwpoly_clone_deep(const LWPOLY *lwgeom); 344 LWCOLLECTION *lwcollection_clone_deep(const LWCOLLECTION *lwgeom); 345 GBOX *gbox_clone(const GBOX *gbox); 346 347 /* 348 * Clockwise 349 */ 350 void lwpoly_force_clockwise(LWPOLY *poly); 351 void lwtriangle_force_clockwise(LWTRIANGLE *triangle); 352 int lwpoly_is_clockwise(LWPOLY *poly); 353 int lwtriangle_is_clockwise(LWTRIANGLE *triangle); 354 int ptarray_isccw(const POINTARRAY *pa); 355 356 /* 357 * Same 358 */ 359 char ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2); 360 char lwpoint_same(const LWPOINT *p1, const LWPOINT *p2); 361 char lwline_same(const LWLINE *p1, const LWLINE *p2); 362 char lwpoly_same(const LWPOLY *p1, const LWPOLY *p2); 363 char lwtriangle_same(const LWTRIANGLE *p1, const LWTRIANGLE *p2); 364 char lwcollection_same(const LWCOLLECTION *p1, const LWCOLLECTION *p2); 365 char lwcircstring_same(const LWCIRCSTRING *p1, const LWCIRCSTRING *p2); 366 367 /* 368 * Shift 369 */ 370 void ptarray_longitude_shift(POINTARRAY *pa); 371 372 /* 373 * Support for in place modification of point arrays, fast 374 * function to move coordinate values around 375 */ 376 void ptarray_copy_point(POINTARRAY *pa, uint32_t from, uint32_t to); 377 378 /* 379 * Reverse 380 */ 381 void ptarray_reverse_in_place(POINTARRAY *pa); 382 383 /* 384 * Startpoint 385 */ 386 int lwpoly_startpoint(const LWPOLY* lwpoly, POINT4D* pt); 387 int ptarray_startpoint(const POINTARRAY* pa, POINT4D* pt); 388 int lwcollection_startpoint(const LWCOLLECTION* col, POINT4D* pt); 389 390 /* 391 * Write into *ret the coordinates of the closest point on 392 * segment A-B to the reference input point R 393 */ 394 void closest_point_on_segment(const POINT4D *R, const POINT4D *A, const POINT4D *B, POINT4D *ret); 395 396 /* 397 * Repeated points 398 */ 399 POINTARRAY *ptarray_remove_repeated_points(const POINTARRAY *in, double tolerance); 400 LWGEOM* lwline_remove_repeated_points(const LWLINE *in, double tolerance); 401 void ptarray_remove_repeated_points_in_place(POINTARRAY *pa, double tolerance, uint32_t min_points); 402 403 /* 404 * Closure test 405 */ 406 int lwline_is_closed(const LWLINE *line); 407 int lwpoly_is_closed(const LWPOLY *poly); 408 int lwcircstring_is_closed(const LWCIRCSTRING *curve); 409 int lwcompound_is_closed(const LWCOMPOUND *curve); 410 int lwpsurface_is_closed(const LWPSURFACE *psurface); 411 int lwtin_is_closed(const LWTIN *tin); 412 413 /** 414 * Snap to grid 415 */ 416 void ptarray_grid_in_place(POINTARRAY *pa, const gridspec *grid); 417 418 /* 419 * What side of the line formed by p1 and p2 does q fall? 420 * Returns -1 for left and 1 for right and 0 for co-linearity 421 */ 422 int lw_segment_side(const POINT2D *p1, const POINT2D *p2, const POINT2D *q); 423 int lw_arc_side(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, const POINT2D *Q); 424 int lw_arc_calculate_gbox_cartesian_2d(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3, GBOX *gbox); 425 double lw_arc_center(const POINT2D *p1, const POINT2D *p2, const POINT2D *p3, POINT2D *result); 426 int lw_pt_in_seg(const POINT2D *P, const POINT2D *A1, const POINT2D *A2); 427 int lw_pt_in_arc(const POINT2D *P, const POINT2D *A1, const POINT2D *A2, const POINT2D *A3); 428 int lw_arc_is_pt(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3); 429 double lw_seg_length(const POINT2D *A1, const POINT2D *A2); 430 double lw_arc_length(const POINT2D *A1, const POINT2D *A2, const POINT2D *A3); 431 int pt_in_ring_2d(const POINT2D *p, const POINTARRAY *ring); 432 int ptarray_contains_point(const POINTARRAY *pa, const POINT2D *pt); 433 int ptarrayarc_contains_point(const POINTARRAY *pa, const POINT2D *pt); 434 int ptarray_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number); 435 int ptarrayarc_contains_point_partial(const POINTARRAY *pa, const POINT2D *pt, int check_closed, int *winding_number); 436 int lwcompound_contains_point(const LWCOMPOUND *comp, const POINT2D *pt); 437 int lwgeom_contains_point(const LWGEOM *geom, const POINT2D *pt); 438 439 /** 440 * Split a line by a point and push components to the provided multiline. 441 * If the point doesn't split the line, push nothing to the container. 442 * Returns 0 if the point is off the line. 443 * Returns 1 if the point is on the line boundary (endpoints). 444 * Return 2 if the point is on the interior of the line (only case in which 445 * a split happens). 446 * 447 * NOTE: the components pushed to the output vector have their SRID stripped 448 */ 449 int lwline_split_by_point_to(const LWLINE* ln, const LWPOINT* pt, LWMLINE* to); 450 451 /** Ensure the collection can hold at least up to ngeoms geometries */ 452 void lwcollection_reserve(LWCOLLECTION *col, uint32_t ngeoms); 453 454 /** Check if subtype is allowed in collectiontype */ 455 int lwcollection_allows_subtype(int collectiontype, int subtype); 456 457 /** GBOX utility functions to figure out coverage/location on the globe */ 458 double gbox_angular_height(const GBOX* gbox); 459 double gbox_angular_width(const GBOX* gbox); 460 int gbox_centroid(const GBOX* gbox, POINT2D* out); 461 462 /* Utilities */ 463 int lwprint_double(double d, int maxdd, char *buf); 464 extern uint8_t MULTITYPE[NUMTYPES]; 465 466 extern lwinterrupt_callback *_lwgeom_interrupt_callback; 467 extern int _lwgeom_interrupt_requested; 468 #define LW_ON_INTERRUPT(x) { \ 469 if ( _lwgeom_interrupt_callback ) { \ 470 (*_lwgeom_interrupt_callback)(); \ 471 } \ 472 if ( _lwgeom_interrupt_requested ) { \ 473 _lwgeom_interrupt_requested = 0; \ 474 lwnotice("liblwgeom code interrupted"); \ 475 x; \ 476 } \ 477 } 478 479 int ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox); 480 int gbox_contains_point2d(const GBOX *g, const POINT2D *p); 481 int lwpoly_contains_point(const LWPOLY *poly, const POINT2D *pt); 482 POINT4D* lwmpoint_extract_points_4d(const LWMPOINT* g, uint32_t* npoints, int* input_empty); 483 char* lwstrdup(const char* a); 484 485 #endif /* _LIBLWGEOM_INTERNAL_H */ 486