1 /*-------------------------------------------------------------------------
2 *
3 * geo_ops.c
4 * 2D geometric operations
5 *
6 * Portions Copyright (c) 1996-2018, PostgreSQL Global Development Group
7 * Portions Copyright (c) 1994, Regents of the University of California
8 *
9 *
10 * IDENTIFICATION
11 * src/backend/utils/adt/geo_ops.c
12 *
13 *-------------------------------------------------------------------------
14 */
15 #include "postgres.h"
16
17 #include <math.h>
18 #include <limits.h>
19 #include <float.h>
20 #include <ctype.h>
21
22 #include "libpq/pqformat.h"
23 #include "miscadmin.h"
24 #include "utils/builtins.h"
25 #include "utils/geo_decls.h"
26
27 #ifndef M_PI
28 #define M_PI 3.14159265358979323846
29 #endif
30
31
32 /*
33 * Internal routines
34 */
35
36 enum path_delim
37 {
38 PATH_NONE, PATH_OPEN, PATH_CLOSED
39 };
40
41 static int point_inside(Point *p, int npts, Point *plist);
42 static int lseg_crossing(double x, double y, double px, double py);
43 static BOX *box_construct(double x1, double x2, double y1, double y2);
44 static BOX *box_fill(BOX *result, double x1, double x2, double y1, double y2);
45 static bool box_ov(BOX *box1, BOX *box2);
46 static double box_ht(BOX *box);
47 static double box_wd(BOX *box);
48 static double circle_ar(CIRCLE *circle);
49 static CIRCLE *circle_copy(CIRCLE *circle);
50 static LINE *line_construct_pm(Point *pt, double m);
51 static void line_construct_pts(LINE *line, Point *pt1, Point *pt2);
52 static bool lseg_intersect_internal(LSEG *l1, LSEG *l2);
53 static double lseg_dt(LSEG *l1, LSEG *l2);
54 static bool on_ps_internal(Point *pt, LSEG *lseg);
55 static void make_bound_box(POLYGON *poly);
56 static bool plist_same(int npts, Point *p1, Point *p2);
57 static Point *point_construct(double x, double y);
58 static Point *point_copy(Point *pt);
59 static double single_decode(char *num, char **endptr_p,
60 const char *type_name, const char *orig_string);
61 static void single_encode(float8 x, StringInfo str);
62 static void pair_decode(char *str, double *x, double *y, char **endptr_p,
63 const char *type_name, const char *orig_string);
64 static void pair_encode(float8 x, float8 y, StringInfo str);
65 static int pair_count(char *s, char delim);
66 static void path_decode(char *str, bool opentype, int npts, Point *p,
67 bool *isopen, char **endptr_p,
68 const char *type_name, const char *orig_string);
69 static char *path_encode(enum path_delim path_delim, int npts, Point *pt);
70 static void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2);
71 static double box_ar(BOX *box);
72 static void box_cn(Point *center, BOX *box);
73 static Point *interpt_sl(LSEG *lseg, LINE *line);
74 static bool has_interpt_sl(LSEG *lseg, LINE *line);
75 static double dist_pl_internal(Point *pt, LINE *line);
76 static double dist_ps_internal(Point *pt, LSEG *lseg);
77 static Point *line_interpt_internal(LINE *l1, LINE *l2);
78 static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start);
79 static Point *lseg_interpt_internal(LSEG *l1, LSEG *l2);
80 static double dist_ppoly_internal(Point *pt, POLYGON *poly);
81
82
83 /*
84 * Delimiters for input and output strings.
85 * LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively.
86 * LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints.
87 */
88
89 #define LDELIM '('
90 #define RDELIM ')'
91 #define DELIM ','
92 #define LDELIM_EP '['
93 #define RDELIM_EP ']'
94 #define LDELIM_C '<'
95 #define RDELIM_C '>'
96
97
98 /*
99 * Geometric data types are composed of points.
100 * This code tries to support a common format throughout the data types,
101 * to allow for more predictable usage and data type conversion.
102 * The fundamental unit is the point. Other units are line segments,
103 * open paths, boxes, closed paths, and polygons (which should be considered
104 * non-intersecting closed paths).
105 *
106 * Data representation is as follows:
107 * point: (x,y)
108 * line segment: [(x1,y1),(x2,y2)]
109 * box: (x1,y1),(x2,y2)
110 * open path: [(x1,y1),...,(xn,yn)]
111 * closed path: ((x1,y1),...,(xn,yn))
112 * polygon: ((x1,y1),...,(xn,yn))
113 *
114 * For boxes, the points are opposite corners with the first point at the top right.
115 * For closed paths and polygons, the points should be reordered to allow
116 * fast and correct equality comparisons.
117 *
118 * XXX perhaps points in complex shapes should be reordered internally
119 * to allow faster internal operations, but should keep track of input order
120 * and restore that order for text output - tgl 97/01/16
121 */
122
123 static double
single_decode(char * num,char ** endptr_p,const char * type_name,const char * orig_string)124 single_decode(char *num, char **endptr_p,
125 const char *type_name, const char *orig_string)
126 {
127 return float8in_internal(num, endptr_p, type_name, orig_string);
128 } /* single_decode() */
129
130 static void
single_encode(float8 x,StringInfo str)131 single_encode(float8 x, StringInfo str)
132 {
133 char *xstr = float8out_internal(x);
134
135 appendStringInfoString(str, xstr);
136 pfree(xstr);
137 } /* single_encode() */
138
139 static void
pair_decode(char * str,double * x,double * y,char ** endptr_p,const char * type_name,const char * orig_string)140 pair_decode(char *str, double *x, double *y, char **endptr_p,
141 const char *type_name, const char *orig_string)
142 {
143 bool has_delim;
144
145 while (isspace((unsigned char) *str))
146 str++;
147 if ((has_delim = (*str == LDELIM)))
148 str++;
149
150 *x = float8in_internal(str, &str, type_name, orig_string);
151
152 if (*str++ != DELIM)
153 ereport(ERROR,
154 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
155 errmsg("invalid input syntax for type %s: \"%s\"",
156 type_name, orig_string)));
157
158 *y = float8in_internal(str, &str, type_name, orig_string);
159
160 if (has_delim)
161 {
162 if (*str++ != RDELIM)
163 ereport(ERROR,
164 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
165 errmsg("invalid input syntax for type %s: \"%s\"",
166 type_name, orig_string)));
167 while (isspace((unsigned char) *str))
168 str++;
169 }
170
171 /* report stopping point if wanted, else complain if not end of string */
172 if (endptr_p)
173 *endptr_p = str;
174 else if (*str != '\0')
175 ereport(ERROR,
176 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
177 errmsg("invalid input syntax for type %s: \"%s\"",
178 type_name, orig_string)));
179 }
180
181 static void
pair_encode(float8 x,float8 y,StringInfo str)182 pair_encode(float8 x, float8 y, StringInfo str)
183 {
184 char *xstr = float8out_internal(x);
185 char *ystr = float8out_internal(y);
186
187 appendStringInfo(str, "%s,%s", xstr, ystr);
188 pfree(xstr);
189 pfree(ystr);
190 }
191
192 static void
path_decode(char * str,bool opentype,int npts,Point * p,bool * isopen,char ** endptr_p,const char * type_name,const char * orig_string)193 path_decode(char *str, bool opentype, int npts, Point *p,
194 bool *isopen, char **endptr_p,
195 const char *type_name, const char *orig_string)
196 {
197 int depth = 0;
198 char *cp;
199 int i;
200
201 while (isspace((unsigned char) *str))
202 str++;
203 if ((*isopen = (*str == LDELIM_EP)))
204 {
205 /* no open delimiter allowed? */
206 if (!opentype)
207 ereport(ERROR,
208 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
209 errmsg("invalid input syntax for type %s: \"%s\"",
210 type_name, orig_string)));
211 depth++;
212 str++;
213 }
214 else if (*str == LDELIM)
215 {
216 cp = (str + 1);
217 while (isspace((unsigned char) *cp))
218 cp++;
219 if (*cp == LDELIM)
220 {
221 depth++;
222 str = cp;
223 }
224 else if (strrchr(str, LDELIM) == str)
225 {
226 depth++;
227 str = cp;
228 }
229 }
230
231 for (i = 0; i < npts; i++)
232 {
233 pair_decode(str, &(p->x), &(p->y), &str, type_name, orig_string);
234 if (*str == DELIM)
235 str++;
236 p++;
237 }
238
239 while (isspace((unsigned char) *str))
240 str++;
241 while (depth > 0)
242 {
243 if ((*str == RDELIM)
244 || ((*str == RDELIM_EP) && (*isopen) && (depth == 1)))
245 {
246 depth--;
247 str++;
248 while (isspace((unsigned char) *str))
249 str++;
250 }
251 else
252 ereport(ERROR,
253 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
254 errmsg("invalid input syntax for type %s: \"%s\"",
255 type_name, orig_string)));
256 }
257
258 /* report stopping point if wanted, else complain if not end of string */
259 if (endptr_p)
260 *endptr_p = str;
261 else if (*str != '\0')
262 ereport(ERROR,
263 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
264 errmsg("invalid input syntax for type %s: \"%s\"",
265 type_name, orig_string)));
266 } /* path_decode() */
267
268 static char *
path_encode(enum path_delim path_delim,int npts,Point * pt)269 path_encode(enum path_delim path_delim, int npts, Point *pt)
270 {
271 StringInfoData str;
272 int i;
273
274 initStringInfo(&str);
275
276 switch (path_delim)
277 {
278 case PATH_CLOSED:
279 appendStringInfoChar(&str, LDELIM);
280 break;
281 case PATH_OPEN:
282 appendStringInfoChar(&str, LDELIM_EP);
283 break;
284 case PATH_NONE:
285 break;
286 }
287
288 for (i = 0; i < npts; i++)
289 {
290 if (i > 0)
291 appendStringInfoChar(&str, DELIM);
292 appendStringInfoChar(&str, LDELIM);
293 pair_encode(pt->x, pt->y, &str);
294 appendStringInfoChar(&str, RDELIM);
295 pt++;
296 }
297
298 switch (path_delim)
299 {
300 case PATH_CLOSED:
301 appendStringInfoChar(&str, RDELIM);
302 break;
303 case PATH_OPEN:
304 appendStringInfoChar(&str, RDELIM_EP);
305 break;
306 case PATH_NONE:
307 break;
308 }
309
310 return str.data;
311 } /* path_encode() */
312
313 /*-------------------------------------------------------------
314 * pair_count - count the number of points
315 * allow the following notation:
316 * '((1,2),(3,4))'
317 * '(1,3,2,4)'
318 * require an odd number of delim characters in the string
319 *-------------------------------------------------------------*/
320 static int
pair_count(char * s,char delim)321 pair_count(char *s, char delim)
322 {
323 int ndelim = 0;
324
325 while ((s = strchr(s, delim)) != NULL)
326 {
327 ndelim++;
328 s++;
329 }
330 return (ndelim % 2) ? ((ndelim + 1) / 2) : -1;
331 }
332
333
334 /***********************************************************************
335 **
336 ** Routines for two-dimensional boxes.
337 **
338 ***********************************************************************/
339
340 /*----------------------------------------------------------
341 * Formatting and conversion routines.
342 *---------------------------------------------------------*/
343
344 /* box_in - convert a string to internal form.
345 *
346 * External format: (two corners of box)
347 * "(f8, f8), (f8, f8)"
348 * also supports the older style "(f8, f8, f8, f8)"
349 */
350 Datum
box_in(PG_FUNCTION_ARGS)351 box_in(PG_FUNCTION_ARGS)
352 {
353 char *str = PG_GETARG_CSTRING(0);
354 BOX *box = (BOX *) palloc(sizeof(BOX));
355 bool isopen;
356 double x,
357 y;
358
359 path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str);
360
361 /* reorder corners if necessary... */
362 if (box->high.x < box->low.x)
363 {
364 x = box->high.x;
365 box->high.x = box->low.x;
366 box->low.x = x;
367 }
368 if (box->high.y < box->low.y)
369 {
370 y = box->high.y;
371 box->high.y = box->low.y;
372 box->low.y = y;
373 }
374
375 PG_RETURN_BOX_P(box);
376 }
377
378 /* box_out - convert a box to external form.
379 */
380 Datum
box_out(PG_FUNCTION_ARGS)381 box_out(PG_FUNCTION_ARGS)
382 {
383 BOX *box = PG_GETARG_BOX_P(0);
384
385 PG_RETURN_CSTRING(path_encode(PATH_NONE, 2, &(box->high)));
386 }
387
388 /*
389 * box_recv - converts external binary format to box
390 */
391 Datum
box_recv(PG_FUNCTION_ARGS)392 box_recv(PG_FUNCTION_ARGS)
393 {
394 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
395 BOX *box;
396 double x,
397 y;
398
399 box = (BOX *) palloc(sizeof(BOX));
400
401 box->high.x = pq_getmsgfloat8(buf);
402 box->high.y = pq_getmsgfloat8(buf);
403 box->low.x = pq_getmsgfloat8(buf);
404 box->low.y = pq_getmsgfloat8(buf);
405
406 /* reorder corners if necessary... */
407 if (box->high.x < box->low.x)
408 {
409 x = box->high.x;
410 box->high.x = box->low.x;
411 box->low.x = x;
412 }
413 if (box->high.y < box->low.y)
414 {
415 y = box->high.y;
416 box->high.y = box->low.y;
417 box->low.y = y;
418 }
419
420 PG_RETURN_BOX_P(box);
421 }
422
423 /*
424 * box_send - converts box to binary format
425 */
426 Datum
box_send(PG_FUNCTION_ARGS)427 box_send(PG_FUNCTION_ARGS)
428 {
429 BOX *box = PG_GETARG_BOX_P(0);
430 StringInfoData buf;
431
432 pq_begintypsend(&buf);
433 pq_sendfloat8(&buf, box->high.x);
434 pq_sendfloat8(&buf, box->high.y);
435 pq_sendfloat8(&buf, box->low.x);
436 pq_sendfloat8(&buf, box->low.y);
437 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
438 }
439
440
441 /* box_construct - fill in a new box.
442 */
443 static BOX *
box_construct(double x1,double x2,double y1,double y2)444 box_construct(double x1, double x2, double y1, double y2)
445 {
446 BOX *result = (BOX *) palloc(sizeof(BOX));
447
448 return box_fill(result, x1, x2, y1, y2);
449 }
450
451
452 /* box_fill - fill in a given box struct
453 */
454 static BOX *
box_fill(BOX * result,double x1,double x2,double y1,double y2)455 box_fill(BOX *result, double x1, double x2, double y1, double y2)
456 {
457 if (x1 > x2)
458 {
459 result->high.x = x1;
460 result->low.x = x2;
461 }
462 else
463 {
464 result->high.x = x2;
465 result->low.x = x1;
466 }
467 if (y1 > y2)
468 {
469 result->high.y = y1;
470 result->low.y = y2;
471 }
472 else
473 {
474 result->high.y = y2;
475 result->low.y = y1;
476 }
477
478 return result;
479 }
480
481
482 /* box_copy - copy a box
483 */
484 BOX *
box_copy(BOX * box)485 box_copy(BOX *box)
486 {
487 BOX *result = (BOX *) palloc(sizeof(BOX));
488
489 memcpy((char *) result, (char *) box, sizeof(BOX));
490
491 return result;
492 }
493
494
495 /*----------------------------------------------------------
496 * Relational operators for BOXes.
497 * <, >, <=, >=, and == are based on box area.
498 *---------------------------------------------------------*/
499
500 /* box_same - are two boxes identical?
501 */
502 Datum
box_same(PG_FUNCTION_ARGS)503 box_same(PG_FUNCTION_ARGS)
504 {
505 BOX *box1 = PG_GETARG_BOX_P(0);
506 BOX *box2 = PG_GETARG_BOX_P(1);
507
508 PG_RETURN_BOOL(FPeq(box1->high.x, box2->high.x) &&
509 FPeq(box1->low.x, box2->low.x) &&
510 FPeq(box1->high.y, box2->high.y) &&
511 FPeq(box1->low.y, box2->low.y));
512 }
513
514 /* box_overlap - does box1 overlap box2?
515 */
516 Datum
box_overlap(PG_FUNCTION_ARGS)517 box_overlap(PG_FUNCTION_ARGS)
518 {
519 BOX *box1 = PG_GETARG_BOX_P(0);
520 BOX *box2 = PG_GETARG_BOX_P(1);
521
522 PG_RETURN_BOOL(box_ov(box1, box2));
523 }
524
525 static bool
box_ov(BOX * box1,BOX * box2)526 box_ov(BOX *box1, BOX *box2)
527 {
528 return (FPle(box1->low.x, box2->high.x) &&
529 FPle(box2->low.x, box1->high.x) &&
530 FPle(box1->low.y, box2->high.y) &&
531 FPle(box2->low.y, box1->high.y));
532 }
533
534 /* box_left - is box1 strictly left of box2?
535 */
536 Datum
box_left(PG_FUNCTION_ARGS)537 box_left(PG_FUNCTION_ARGS)
538 {
539 BOX *box1 = PG_GETARG_BOX_P(0);
540 BOX *box2 = PG_GETARG_BOX_P(1);
541
542 PG_RETURN_BOOL(FPlt(box1->high.x, box2->low.x));
543 }
544
545 /* box_overleft - is the right edge of box1 at or left of
546 * the right edge of box2?
547 *
548 * This is "less than or equal" for the end of a time range,
549 * when time ranges are stored as rectangles.
550 */
551 Datum
box_overleft(PG_FUNCTION_ARGS)552 box_overleft(PG_FUNCTION_ARGS)
553 {
554 BOX *box1 = PG_GETARG_BOX_P(0);
555 BOX *box2 = PG_GETARG_BOX_P(1);
556
557 PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x));
558 }
559
560 /* box_right - is box1 strictly right of box2?
561 */
562 Datum
box_right(PG_FUNCTION_ARGS)563 box_right(PG_FUNCTION_ARGS)
564 {
565 BOX *box1 = PG_GETARG_BOX_P(0);
566 BOX *box2 = PG_GETARG_BOX_P(1);
567
568 PG_RETURN_BOOL(FPgt(box1->low.x, box2->high.x));
569 }
570
571 /* box_overright - is the left edge of box1 at or right of
572 * the left edge of box2?
573 *
574 * This is "greater than or equal" for time ranges, when time ranges
575 * are stored as rectangles.
576 */
577 Datum
box_overright(PG_FUNCTION_ARGS)578 box_overright(PG_FUNCTION_ARGS)
579 {
580 BOX *box1 = PG_GETARG_BOX_P(0);
581 BOX *box2 = PG_GETARG_BOX_P(1);
582
583 PG_RETURN_BOOL(FPge(box1->low.x, box2->low.x));
584 }
585
586 /* box_below - is box1 strictly below box2?
587 */
588 Datum
box_below(PG_FUNCTION_ARGS)589 box_below(PG_FUNCTION_ARGS)
590 {
591 BOX *box1 = PG_GETARG_BOX_P(0);
592 BOX *box2 = PG_GETARG_BOX_P(1);
593
594 PG_RETURN_BOOL(FPlt(box1->high.y, box2->low.y));
595 }
596
597 /* box_overbelow - is the upper edge of box1 at or below
598 * the upper edge of box2?
599 */
600 Datum
box_overbelow(PG_FUNCTION_ARGS)601 box_overbelow(PG_FUNCTION_ARGS)
602 {
603 BOX *box1 = PG_GETARG_BOX_P(0);
604 BOX *box2 = PG_GETARG_BOX_P(1);
605
606 PG_RETURN_BOOL(FPle(box1->high.y, box2->high.y));
607 }
608
609 /* box_above - is box1 strictly above box2?
610 */
611 Datum
box_above(PG_FUNCTION_ARGS)612 box_above(PG_FUNCTION_ARGS)
613 {
614 BOX *box1 = PG_GETARG_BOX_P(0);
615 BOX *box2 = PG_GETARG_BOX_P(1);
616
617 PG_RETURN_BOOL(FPgt(box1->low.y, box2->high.y));
618 }
619
620 /* box_overabove - is the lower edge of box1 at or above
621 * the lower edge of box2?
622 */
623 Datum
box_overabove(PG_FUNCTION_ARGS)624 box_overabove(PG_FUNCTION_ARGS)
625 {
626 BOX *box1 = PG_GETARG_BOX_P(0);
627 BOX *box2 = PG_GETARG_BOX_P(1);
628
629 PG_RETURN_BOOL(FPge(box1->low.y, box2->low.y));
630 }
631
632 /* box_contained - is box1 contained by box2?
633 */
634 Datum
box_contained(PG_FUNCTION_ARGS)635 box_contained(PG_FUNCTION_ARGS)
636 {
637 BOX *box1 = PG_GETARG_BOX_P(0);
638 BOX *box2 = PG_GETARG_BOX_P(1);
639
640 PG_RETURN_BOOL(FPle(box1->high.x, box2->high.x) &&
641 FPge(box1->low.x, box2->low.x) &&
642 FPle(box1->high.y, box2->high.y) &&
643 FPge(box1->low.y, box2->low.y));
644 }
645
646 /* box_contain - does box1 contain box2?
647 */
648 Datum
box_contain(PG_FUNCTION_ARGS)649 box_contain(PG_FUNCTION_ARGS)
650 {
651 BOX *box1 = PG_GETARG_BOX_P(0);
652 BOX *box2 = PG_GETARG_BOX_P(1);
653
654 PG_RETURN_BOOL(FPge(box1->high.x, box2->high.x) &&
655 FPle(box1->low.x, box2->low.x) &&
656 FPge(box1->high.y, box2->high.y) &&
657 FPle(box1->low.y, box2->low.y));
658 }
659
660
661 /* box_positionop -
662 * is box1 entirely {above,below} box2?
663 *
664 * box_below_eq and box_above_eq are obsolete versions that (probably
665 * erroneously) accept the equal-boundaries case. Since these are not
666 * in sync with the box_left and box_right code, they are deprecated and
667 * not supported in the PG 8.1 rtree operator class extension.
668 */
669 Datum
box_below_eq(PG_FUNCTION_ARGS)670 box_below_eq(PG_FUNCTION_ARGS)
671 {
672 BOX *box1 = PG_GETARG_BOX_P(0);
673 BOX *box2 = PG_GETARG_BOX_P(1);
674
675 PG_RETURN_BOOL(FPle(box1->high.y, box2->low.y));
676 }
677
678 Datum
box_above_eq(PG_FUNCTION_ARGS)679 box_above_eq(PG_FUNCTION_ARGS)
680 {
681 BOX *box1 = PG_GETARG_BOX_P(0);
682 BOX *box2 = PG_GETARG_BOX_P(1);
683
684 PG_RETURN_BOOL(FPge(box1->low.y, box2->high.y));
685 }
686
687
688 /* box_relop - is area(box1) relop area(box2), within
689 * our accuracy constraint?
690 */
691 Datum
box_lt(PG_FUNCTION_ARGS)692 box_lt(PG_FUNCTION_ARGS)
693 {
694 BOX *box1 = PG_GETARG_BOX_P(0);
695 BOX *box2 = PG_GETARG_BOX_P(1);
696
697 PG_RETURN_BOOL(FPlt(box_ar(box1), box_ar(box2)));
698 }
699
700 Datum
box_gt(PG_FUNCTION_ARGS)701 box_gt(PG_FUNCTION_ARGS)
702 {
703 BOX *box1 = PG_GETARG_BOX_P(0);
704 BOX *box2 = PG_GETARG_BOX_P(1);
705
706 PG_RETURN_BOOL(FPgt(box_ar(box1), box_ar(box2)));
707 }
708
709 Datum
box_eq(PG_FUNCTION_ARGS)710 box_eq(PG_FUNCTION_ARGS)
711 {
712 BOX *box1 = PG_GETARG_BOX_P(0);
713 BOX *box2 = PG_GETARG_BOX_P(1);
714
715 PG_RETURN_BOOL(FPeq(box_ar(box1), box_ar(box2)));
716 }
717
718 Datum
box_le(PG_FUNCTION_ARGS)719 box_le(PG_FUNCTION_ARGS)
720 {
721 BOX *box1 = PG_GETARG_BOX_P(0);
722 BOX *box2 = PG_GETARG_BOX_P(1);
723
724 PG_RETURN_BOOL(FPle(box_ar(box1), box_ar(box2)));
725 }
726
727 Datum
box_ge(PG_FUNCTION_ARGS)728 box_ge(PG_FUNCTION_ARGS)
729 {
730 BOX *box1 = PG_GETARG_BOX_P(0);
731 BOX *box2 = PG_GETARG_BOX_P(1);
732
733 PG_RETURN_BOOL(FPge(box_ar(box1), box_ar(box2)));
734 }
735
736
737 /*----------------------------------------------------------
738 * "Arithmetic" operators on boxes.
739 *---------------------------------------------------------*/
740
741 /* box_area - returns the area of the box.
742 */
743 Datum
box_area(PG_FUNCTION_ARGS)744 box_area(PG_FUNCTION_ARGS)
745 {
746 BOX *box = PG_GETARG_BOX_P(0);
747
748 PG_RETURN_FLOAT8(box_ar(box));
749 }
750
751
752 /* box_width - returns the width of the box
753 * (horizontal magnitude).
754 */
755 Datum
box_width(PG_FUNCTION_ARGS)756 box_width(PG_FUNCTION_ARGS)
757 {
758 BOX *box = PG_GETARG_BOX_P(0);
759
760 PG_RETURN_FLOAT8(box->high.x - box->low.x);
761 }
762
763
764 /* box_height - returns the height of the box
765 * (vertical magnitude).
766 */
767 Datum
box_height(PG_FUNCTION_ARGS)768 box_height(PG_FUNCTION_ARGS)
769 {
770 BOX *box = PG_GETARG_BOX_P(0);
771
772 PG_RETURN_FLOAT8(box->high.y - box->low.y);
773 }
774
775
776 /* box_distance - returns the distance between the
777 * center points of two boxes.
778 */
779 Datum
box_distance(PG_FUNCTION_ARGS)780 box_distance(PG_FUNCTION_ARGS)
781 {
782 BOX *box1 = PG_GETARG_BOX_P(0);
783 BOX *box2 = PG_GETARG_BOX_P(1);
784 Point a,
785 b;
786
787 box_cn(&a, box1);
788 box_cn(&b, box2);
789
790 PG_RETURN_FLOAT8(HYPOT(a.x - b.x, a.y - b.y));
791 }
792
793
794 /* box_center - returns the center point of the box.
795 */
796 Datum
box_center(PG_FUNCTION_ARGS)797 box_center(PG_FUNCTION_ARGS)
798 {
799 BOX *box = PG_GETARG_BOX_P(0);
800 Point *result = (Point *) palloc(sizeof(Point));
801
802 box_cn(result, box);
803
804 PG_RETURN_POINT_P(result);
805 }
806
807
808 /* box_ar - returns the area of the box.
809 */
810 static double
box_ar(BOX * box)811 box_ar(BOX *box)
812 {
813 return box_wd(box) * box_ht(box);
814 }
815
816
817 /* box_cn - stores the centerpoint of the box into *center.
818 */
819 static void
box_cn(Point * center,BOX * box)820 box_cn(Point *center, BOX *box)
821 {
822 center->x = (box->high.x + box->low.x) / 2.0;
823 center->y = (box->high.y + box->low.y) / 2.0;
824 }
825
826
827 /* box_wd - returns the width (length) of the box
828 * (horizontal magnitude).
829 */
830 static double
box_wd(BOX * box)831 box_wd(BOX *box)
832 {
833 return box->high.x - box->low.x;
834 }
835
836
837 /* box_ht - returns the height of the box
838 * (vertical magnitude).
839 */
840 static double
box_ht(BOX * box)841 box_ht(BOX *box)
842 {
843 return box->high.y - box->low.y;
844 }
845
846
847 /*----------------------------------------------------------
848 * Funky operations.
849 *---------------------------------------------------------*/
850
851 /* box_intersect -
852 * returns the overlapping portion of two boxes,
853 * or NULL if they do not intersect.
854 */
855 Datum
box_intersect(PG_FUNCTION_ARGS)856 box_intersect(PG_FUNCTION_ARGS)
857 {
858 BOX *box1 = PG_GETARG_BOX_P(0);
859 BOX *box2 = PG_GETARG_BOX_P(1);
860 BOX *result;
861
862 if (!box_ov(box1, box2))
863 PG_RETURN_NULL();
864
865 result = (BOX *) palloc(sizeof(BOX));
866
867 result->high.x = Min(box1->high.x, box2->high.x);
868 result->low.x = Max(box1->low.x, box2->low.x);
869 result->high.y = Min(box1->high.y, box2->high.y);
870 result->low.y = Max(box1->low.y, box2->low.y);
871
872 PG_RETURN_BOX_P(result);
873 }
874
875
876 /* box_diagonal -
877 * returns a line segment which happens to be the
878 * positive-slope diagonal of "box".
879 */
880 Datum
box_diagonal(PG_FUNCTION_ARGS)881 box_diagonal(PG_FUNCTION_ARGS)
882 {
883 BOX *box = PG_GETARG_BOX_P(0);
884 LSEG *result = (LSEG *) palloc(sizeof(LSEG));
885
886 statlseg_construct(result, &box->high, &box->low);
887
888 PG_RETURN_LSEG_P(result);
889 }
890
891 /***********************************************************************
892 **
893 ** Routines for 2D lines.
894 **
895 ***********************************************************************/
896
897 static bool
line_decode(char * s,const char * str,LINE * line)898 line_decode(char *s, const char *str, LINE *line)
899 {
900 /* s was already advanced over leading '{' */
901 line->A = single_decode(s, &s, "line", str);
902 if (*s++ != DELIM)
903 return false;
904 line->B = single_decode(s, &s, "line", str);
905 if (*s++ != DELIM)
906 return false;
907 line->C = single_decode(s, &s, "line", str);
908 if (*s++ != '}')
909 return false;
910 while (isspace((unsigned char) *s))
911 s++;
912 if (*s != '\0')
913 return false;
914 return true;
915 }
916
917 Datum
line_in(PG_FUNCTION_ARGS)918 line_in(PG_FUNCTION_ARGS)
919 {
920 char *str = PG_GETARG_CSTRING(0);
921 LINE *line = (LINE *) palloc(sizeof(LINE));
922 LSEG lseg;
923 bool isopen;
924 char *s;
925
926 s = str;
927 while (isspace((unsigned char) *s))
928 s++;
929 if (*s == '{')
930 {
931 if (!line_decode(s + 1, str, line))
932 ereport(ERROR,
933 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
934 errmsg("invalid input syntax for type %s: \"%s\"",
935 "line", str)));
936 if (FPzero(line->A) && FPzero(line->B))
937 ereport(ERROR,
938 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
939 errmsg("invalid line specification: A and B cannot both be zero")));
940 }
941 else
942 {
943 path_decode(s, true, 2, &(lseg.p[0]), &isopen, NULL, "line", str);
944 if (FPeq(lseg.p[0].x, lseg.p[1].x) && FPeq(lseg.p[0].y, lseg.p[1].y))
945 ereport(ERROR,
946 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
947 errmsg("invalid line specification: must be two distinct points")));
948 line_construct_pts(line, &lseg.p[0], &lseg.p[1]);
949 }
950
951 PG_RETURN_LINE_P(line);
952 }
953
954
955 Datum
line_out(PG_FUNCTION_ARGS)956 line_out(PG_FUNCTION_ARGS)
957 {
958 LINE *line = PG_GETARG_LINE_P(0);
959 char *astr = float8out_internal(line->A);
960 char *bstr = float8out_internal(line->B);
961 char *cstr = float8out_internal(line->C);
962
963 PG_RETURN_CSTRING(psprintf("{%s,%s,%s}", astr, bstr, cstr));
964 }
965
966 /*
967 * line_recv - converts external binary format to line
968 */
969 Datum
line_recv(PG_FUNCTION_ARGS)970 line_recv(PG_FUNCTION_ARGS)
971 {
972 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
973 LINE *line;
974
975 line = (LINE *) palloc(sizeof(LINE));
976
977 line->A = pq_getmsgfloat8(buf);
978 line->B = pq_getmsgfloat8(buf);
979 line->C = pq_getmsgfloat8(buf);
980
981 PG_RETURN_LINE_P(line);
982 }
983
984 /*
985 * line_send - converts line to binary format
986 */
987 Datum
line_send(PG_FUNCTION_ARGS)988 line_send(PG_FUNCTION_ARGS)
989 {
990 LINE *line = PG_GETARG_LINE_P(0);
991 StringInfoData buf;
992
993 pq_begintypsend(&buf);
994 pq_sendfloat8(&buf, line->A);
995 pq_sendfloat8(&buf, line->B);
996 pq_sendfloat8(&buf, line->C);
997 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
998 }
999
1000
1001 /*----------------------------------------------------------
1002 * Conversion routines from one line formula to internal.
1003 * Internal form: Ax+By+C=0
1004 *---------------------------------------------------------*/
1005
1006 /* line_construct_pm()
1007 * point-slope
1008 */
1009 static LINE *
line_construct_pm(Point * pt,double m)1010 line_construct_pm(Point *pt, double m)
1011 {
1012 LINE *result = (LINE *) palloc(sizeof(LINE));
1013
1014 if (m == DBL_MAX)
1015 {
1016 /* vertical - use "x = C" */
1017 result->A = -1;
1018 result->B = 0;
1019 result->C = pt->x;
1020 }
1021 else
1022 {
1023 /* use "mx - y + yinter = 0" */
1024 result->A = m;
1025 result->B = -1.0;
1026 result->C = pt->y - m * pt->x;
1027 }
1028
1029 return result;
1030 }
1031
1032 /*
1033 * Fill already-allocated LINE struct from two points on the line
1034 */
1035 static void
line_construct_pts(LINE * line,Point * pt1,Point * pt2)1036 line_construct_pts(LINE *line, Point *pt1, Point *pt2)
1037 {
1038 if (FPeq(pt1->x, pt2->x))
1039 { /* vertical */
1040 /* use "x = C" */
1041 line->A = -1;
1042 line->B = 0;
1043 line->C = pt1->x;
1044 #ifdef GEODEBUG
1045 printf("line_construct_pts- line is vertical\n");
1046 #endif
1047 }
1048 else if (FPeq(pt1->y, pt2->y))
1049 { /* horizontal */
1050 /* use "y = C" */
1051 line->A = 0;
1052 line->B = -1;
1053 line->C = pt1->y;
1054 #ifdef GEODEBUG
1055 printf("line_construct_pts- line is horizontal\n");
1056 #endif
1057 }
1058 else
1059 {
1060 /* use "mx - y + yinter = 0" */
1061 line->A = (pt2->y - pt1->y) / (pt2->x - pt1->x);
1062 line->B = -1.0;
1063 line->C = pt1->y - line->A * pt1->x;
1064 /* on some platforms, the preceding expression tends to produce -0 */
1065 if (line->C == 0.0)
1066 line->C = 0.0;
1067 #ifdef GEODEBUG
1068 printf("line_construct_pts- line is neither vertical nor horizontal (diffs x=%.*g, y=%.*g\n",
1069 DBL_DIG, (pt2->x - pt1->x), DBL_DIG, (pt2->y - pt1->y));
1070 #endif
1071 }
1072 }
1073
1074 /* line_construct_pp()
1075 * two points
1076 */
1077 Datum
line_construct_pp(PG_FUNCTION_ARGS)1078 line_construct_pp(PG_FUNCTION_ARGS)
1079 {
1080 Point *pt1 = PG_GETARG_POINT_P(0);
1081 Point *pt2 = PG_GETARG_POINT_P(1);
1082 LINE *result = (LINE *) palloc(sizeof(LINE));
1083
1084 line_construct_pts(result, pt1, pt2);
1085 PG_RETURN_LINE_P(result);
1086 }
1087
1088
1089 /*----------------------------------------------------------
1090 * Relative position routines.
1091 *---------------------------------------------------------*/
1092
1093 Datum
line_intersect(PG_FUNCTION_ARGS)1094 line_intersect(PG_FUNCTION_ARGS)
1095 {
1096 LINE *l1 = PG_GETARG_LINE_P(0);
1097 LINE *l2 = PG_GETARG_LINE_P(1);
1098
1099 PG_RETURN_BOOL(!DatumGetBool(DirectFunctionCall2(line_parallel,
1100 LinePGetDatum(l1),
1101 LinePGetDatum(l2))));
1102 }
1103
1104 Datum
line_parallel(PG_FUNCTION_ARGS)1105 line_parallel(PG_FUNCTION_ARGS)
1106 {
1107 LINE *l1 = PG_GETARG_LINE_P(0);
1108 LINE *l2 = PG_GETARG_LINE_P(1);
1109
1110 if (FPzero(l1->B))
1111 PG_RETURN_BOOL(FPzero(l2->B));
1112
1113 PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B)));
1114 }
1115
1116 Datum
line_perp(PG_FUNCTION_ARGS)1117 line_perp(PG_FUNCTION_ARGS)
1118 {
1119 LINE *l1 = PG_GETARG_LINE_P(0);
1120 LINE *l2 = PG_GETARG_LINE_P(1);
1121
1122 if (FPzero(l1->A))
1123 PG_RETURN_BOOL(FPzero(l2->B));
1124 else if (FPzero(l1->B))
1125 PG_RETURN_BOOL(FPzero(l2->A));
1126
1127 PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0));
1128 }
1129
1130 Datum
line_vertical(PG_FUNCTION_ARGS)1131 line_vertical(PG_FUNCTION_ARGS)
1132 {
1133 LINE *line = PG_GETARG_LINE_P(0);
1134
1135 PG_RETURN_BOOL(FPzero(line->B));
1136 }
1137
1138 Datum
line_horizontal(PG_FUNCTION_ARGS)1139 line_horizontal(PG_FUNCTION_ARGS)
1140 {
1141 LINE *line = PG_GETARG_LINE_P(0);
1142
1143 PG_RETURN_BOOL(FPzero(line->A));
1144 }
1145
1146 Datum
line_eq(PG_FUNCTION_ARGS)1147 line_eq(PG_FUNCTION_ARGS)
1148 {
1149 LINE *l1 = PG_GETARG_LINE_P(0);
1150 LINE *l2 = PG_GETARG_LINE_P(1);
1151 double k;
1152
1153 if (!FPzero(l2->A))
1154 k = l1->A / l2->A;
1155 else if (!FPzero(l2->B))
1156 k = l1->B / l2->B;
1157 else if (!FPzero(l2->C))
1158 k = l1->C / l2->C;
1159 else
1160 k = 1.0;
1161
1162 PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) &&
1163 FPeq(l1->B, k * l2->B) &&
1164 FPeq(l1->C, k * l2->C));
1165 }
1166
1167
1168 /*----------------------------------------------------------
1169 * Line arithmetic routines.
1170 *---------------------------------------------------------*/
1171
1172 /* line_distance()
1173 * Distance between two lines.
1174 */
1175 Datum
line_distance(PG_FUNCTION_ARGS)1176 line_distance(PG_FUNCTION_ARGS)
1177 {
1178 LINE *l1 = PG_GETARG_LINE_P(0);
1179 LINE *l2 = PG_GETARG_LINE_P(1);
1180 float8 result;
1181 Point *tmp;
1182
1183 if (!DatumGetBool(DirectFunctionCall2(line_parallel,
1184 LinePGetDatum(l1),
1185 LinePGetDatum(l2))))
1186 PG_RETURN_FLOAT8(0.0);
1187 if (FPzero(l1->B)) /* vertical? */
1188 PG_RETURN_FLOAT8(fabs(l1->C - l2->C));
1189 tmp = point_construct(0.0, l1->C);
1190 result = dist_pl_internal(tmp, l2);
1191 PG_RETURN_FLOAT8(result);
1192 }
1193
1194 /* line_interpt()
1195 * Point where two lines l1, l2 intersect (if any)
1196 */
1197 Datum
line_interpt(PG_FUNCTION_ARGS)1198 line_interpt(PG_FUNCTION_ARGS)
1199 {
1200 LINE *l1 = PG_GETARG_LINE_P(0);
1201 LINE *l2 = PG_GETARG_LINE_P(1);
1202 Point *result;
1203
1204 result = line_interpt_internal(l1, l2);
1205
1206 if (result == NULL)
1207 PG_RETURN_NULL();
1208 PG_RETURN_POINT_P(result);
1209 }
1210
1211 /*
1212 * Internal version of line_interpt
1213 *
1214 * returns a NULL pointer if no intersection point
1215 */
1216 static Point *
line_interpt_internal(LINE * l1,LINE * l2)1217 line_interpt_internal(LINE *l1, LINE *l2)
1218 {
1219 Point *result;
1220 double x,
1221 y;
1222
1223 /*
1224 * NOTE: if the lines are identical then we will find they are parallel
1225 * and report "no intersection". This is a little weird, but since
1226 * there's no *unique* intersection, maybe it's appropriate behavior.
1227 */
1228 if (DatumGetBool(DirectFunctionCall2(line_parallel,
1229 LinePGetDatum(l1),
1230 LinePGetDatum(l2))))
1231 return NULL;
1232
1233 if (FPzero(l1->B)) /* l1 vertical? */
1234 {
1235 x = l1->C;
1236 y = (l2->A * x + l2->C);
1237 }
1238 else if (FPzero(l2->B)) /* l2 vertical? */
1239 {
1240 x = l2->C;
1241 y = (l1->A * x + l1->C);
1242 }
1243 else
1244 {
1245 x = (l1->C - l2->C) / (l2->A - l1->A);
1246 y = (l1->A * x + l1->C);
1247 }
1248 result = point_construct(x, y);
1249
1250 #ifdef GEODEBUG
1251 printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n",
1252 DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C);
1253 printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y);
1254 #endif
1255
1256 return result;
1257 }
1258
1259
1260 /***********************************************************************
1261 **
1262 ** Routines for 2D paths (sequences of line segments, also
1263 ** called `polylines').
1264 **
1265 ** This is not a general package for geometric paths,
1266 ** which of course include polygons; the emphasis here
1267 ** is on (for example) usefulness in wire layout.
1268 **
1269 ***********************************************************************/
1270
1271 /*----------------------------------------------------------
1272 * String to path / path to string conversion.
1273 * External format:
1274 * "((xcoord, ycoord),... )"
1275 * "[(xcoord, ycoord),... ]"
1276 * "(xcoord, ycoord),... "
1277 * "[xcoord, ycoord,... ]"
1278 * Also support older format:
1279 * "(closed, npts, xcoord, ycoord,... )"
1280 *---------------------------------------------------------*/
1281
1282 Datum
path_area(PG_FUNCTION_ARGS)1283 path_area(PG_FUNCTION_ARGS)
1284 {
1285 PATH *path = PG_GETARG_PATH_P(0);
1286 double area = 0.0;
1287 int i,
1288 j;
1289
1290 if (!path->closed)
1291 PG_RETURN_NULL();
1292
1293 for (i = 0; i < path->npts; i++)
1294 {
1295 j = (i + 1) % path->npts;
1296 area += path->p[i].x * path->p[j].y;
1297 area -= path->p[i].y * path->p[j].x;
1298 }
1299
1300 area *= 0.5;
1301 PG_RETURN_FLOAT8(area < 0.0 ? -area : area);
1302 }
1303
1304
1305 Datum
path_in(PG_FUNCTION_ARGS)1306 path_in(PG_FUNCTION_ARGS)
1307 {
1308 char *str = PG_GETARG_CSTRING(0);
1309 PATH *path;
1310 bool isopen;
1311 char *s;
1312 int npts;
1313 int size;
1314 int base_size;
1315 int depth = 0;
1316
1317 if ((npts = pair_count(str, ',')) <= 0)
1318 ereport(ERROR,
1319 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1320 errmsg("invalid input syntax for type %s: \"%s\"",
1321 "path", str)));
1322
1323 s = str;
1324 while (isspace((unsigned char) *s))
1325 s++;
1326
1327 /* skip single leading paren */
1328 if ((*s == LDELIM) && (strrchr(s, LDELIM) == s))
1329 {
1330 s++;
1331 depth++;
1332 }
1333
1334 base_size = sizeof(path->p[0]) * npts;
1335 size = offsetof(PATH, p) + base_size;
1336
1337 /* Check for integer overflow */
1338 if (base_size / npts != sizeof(path->p[0]) || size <= base_size)
1339 ereport(ERROR,
1340 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1341 errmsg("too many points requested")));
1342
1343 path = (PATH *) palloc(size);
1344
1345 SET_VARSIZE(path, size);
1346 path->npts = npts;
1347
1348 path_decode(s, true, npts, &(path->p[0]), &isopen, &s, "path", str);
1349
1350 if (depth >= 1)
1351 {
1352 if (*s++ != RDELIM)
1353 ereport(ERROR,
1354 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1355 errmsg("invalid input syntax for type %s: \"%s\"",
1356 "path", str)));
1357 while (isspace((unsigned char) *s))
1358 s++;
1359 }
1360 if (*s != '\0')
1361 ereport(ERROR,
1362 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1363 errmsg("invalid input syntax for type %s: \"%s\"",
1364 "path", str)));
1365
1366 path->closed = (!isopen);
1367 /* prevent instability in unused pad bytes */
1368 path->dummy = 0;
1369
1370 PG_RETURN_PATH_P(path);
1371 }
1372
1373
1374 Datum
path_out(PG_FUNCTION_ARGS)1375 path_out(PG_FUNCTION_ARGS)
1376 {
1377 PATH *path = PG_GETARG_PATH_P(0);
1378
1379 PG_RETURN_CSTRING(path_encode(path->closed ? PATH_CLOSED : PATH_OPEN, path->npts, path->p));
1380 }
1381
1382 /*
1383 * path_recv - converts external binary format to path
1384 *
1385 * External representation is closed flag (a boolean byte), int32 number
1386 * of points, and the points.
1387 */
1388 Datum
path_recv(PG_FUNCTION_ARGS)1389 path_recv(PG_FUNCTION_ARGS)
1390 {
1391 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1392 PATH *path;
1393 int closed;
1394 int32 npts;
1395 int32 i;
1396 int size;
1397
1398 closed = pq_getmsgbyte(buf);
1399 npts = pq_getmsgint(buf, sizeof(int32));
1400 if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(PATH, p)) / sizeof(Point)))
1401 ereport(ERROR,
1402 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
1403 errmsg("invalid number of points in external \"path\" value")));
1404
1405 size = offsetof(PATH, p) + sizeof(path->p[0]) * npts;
1406 path = (PATH *) palloc(size);
1407
1408 SET_VARSIZE(path, size);
1409 path->npts = npts;
1410 path->closed = (closed ? 1 : 0);
1411 /* prevent instability in unused pad bytes */
1412 path->dummy = 0;
1413
1414 for (i = 0; i < npts; i++)
1415 {
1416 path->p[i].x = pq_getmsgfloat8(buf);
1417 path->p[i].y = pq_getmsgfloat8(buf);
1418 }
1419
1420 PG_RETURN_PATH_P(path);
1421 }
1422
1423 /*
1424 * path_send - converts path to binary format
1425 */
1426 Datum
path_send(PG_FUNCTION_ARGS)1427 path_send(PG_FUNCTION_ARGS)
1428 {
1429 PATH *path = PG_GETARG_PATH_P(0);
1430 StringInfoData buf;
1431 int32 i;
1432
1433 pq_begintypsend(&buf);
1434 pq_sendbyte(&buf, path->closed ? 1 : 0);
1435 pq_sendint32(&buf, path->npts);
1436 for (i = 0; i < path->npts; i++)
1437 {
1438 pq_sendfloat8(&buf, path->p[i].x);
1439 pq_sendfloat8(&buf, path->p[i].y);
1440 }
1441 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1442 }
1443
1444
1445 /*----------------------------------------------------------
1446 * Relational operators.
1447 * These are based on the path cardinality,
1448 * as stupid as that sounds.
1449 *
1450 * Better relops and access methods coming soon.
1451 *---------------------------------------------------------*/
1452
1453 Datum
path_n_lt(PG_FUNCTION_ARGS)1454 path_n_lt(PG_FUNCTION_ARGS)
1455 {
1456 PATH *p1 = PG_GETARG_PATH_P(0);
1457 PATH *p2 = PG_GETARG_PATH_P(1);
1458
1459 PG_RETURN_BOOL(p1->npts < p2->npts);
1460 }
1461
1462 Datum
path_n_gt(PG_FUNCTION_ARGS)1463 path_n_gt(PG_FUNCTION_ARGS)
1464 {
1465 PATH *p1 = PG_GETARG_PATH_P(0);
1466 PATH *p2 = PG_GETARG_PATH_P(1);
1467
1468 PG_RETURN_BOOL(p1->npts > p2->npts);
1469 }
1470
1471 Datum
path_n_eq(PG_FUNCTION_ARGS)1472 path_n_eq(PG_FUNCTION_ARGS)
1473 {
1474 PATH *p1 = PG_GETARG_PATH_P(0);
1475 PATH *p2 = PG_GETARG_PATH_P(1);
1476
1477 PG_RETURN_BOOL(p1->npts == p2->npts);
1478 }
1479
1480 Datum
path_n_le(PG_FUNCTION_ARGS)1481 path_n_le(PG_FUNCTION_ARGS)
1482 {
1483 PATH *p1 = PG_GETARG_PATH_P(0);
1484 PATH *p2 = PG_GETARG_PATH_P(1);
1485
1486 PG_RETURN_BOOL(p1->npts <= p2->npts);
1487 }
1488
1489 Datum
path_n_ge(PG_FUNCTION_ARGS)1490 path_n_ge(PG_FUNCTION_ARGS)
1491 {
1492 PATH *p1 = PG_GETARG_PATH_P(0);
1493 PATH *p2 = PG_GETARG_PATH_P(1);
1494
1495 PG_RETURN_BOOL(p1->npts >= p2->npts);
1496 }
1497
1498 /*----------------------------------------------------------
1499 * Conversion operators.
1500 *---------------------------------------------------------*/
1501
1502 Datum
path_isclosed(PG_FUNCTION_ARGS)1503 path_isclosed(PG_FUNCTION_ARGS)
1504 {
1505 PATH *path = PG_GETARG_PATH_P(0);
1506
1507 PG_RETURN_BOOL(path->closed);
1508 }
1509
1510 Datum
path_isopen(PG_FUNCTION_ARGS)1511 path_isopen(PG_FUNCTION_ARGS)
1512 {
1513 PATH *path = PG_GETARG_PATH_P(0);
1514
1515 PG_RETURN_BOOL(!path->closed);
1516 }
1517
1518 Datum
path_npoints(PG_FUNCTION_ARGS)1519 path_npoints(PG_FUNCTION_ARGS)
1520 {
1521 PATH *path = PG_GETARG_PATH_P(0);
1522
1523 PG_RETURN_INT32(path->npts);
1524 }
1525
1526
1527 Datum
path_close(PG_FUNCTION_ARGS)1528 path_close(PG_FUNCTION_ARGS)
1529 {
1530 PATH *path = PG_GETARG_PATH_P_COPY(0);
1531
1532 path->closed = true;
1533
1534 PG_RETURN_PATH_P(path);
1535 }
1536
1537 Datum
path_open(PG_FUNCTION_ARGS)1538 path_open(PG_FUNCTION_ARGS)
1539 {
1540 PATH *path = PG_GETARG_PATH_P_COPY(0);
1541
1542 path->closed = false;
1543
1544 PG_RETURN_PATH_P(path);
1545 }
1546
1547
1548 /* path_inter -
1549 * Does p1 intersect p2 at any point?
1550 * Use bounding boxes for a quick (O(n)) check, then do a
1551 * O(n^2) iterative edge check.
1552 */
1553 Datum
path_inter(PG_FUNCTION_ARGS)1554 path_inter(PG_FUNCTION_ARGS)
1555 {
1556 PATH *p1 = PG_GETARG_PATH_P(0);
1557 PATH *p2 = PG_GETARG_PATH_P(1);
1558 BOX b1,
1559 b2;
1560 int i,
1561 j;
1562 LSEG seg1,
1563 seg2;
1564
1565 if (p1->npts <= 0 || p2->npts <= 0)
1566 PG_RETURN_BOOL(false);
1567
1568 b1.high.x = b1.low.x = p1->p[0].x;
1569 b1.high.y = b1.low.y = p1->p[0].y;
1570 for (i = 1; i < p1->npts; i++)
1571 {
1572 b1.high.x = Max(p1->p[i].x, b1.high.x);
1573 b1.high.y = Max(p1->p[i].y, b1.high.y);
1574 b1.low.x = Min(p1->p[i].x, b1.low.x);
1575 b1.low.y = Min(p1->p[i].y, b1.low.y);
1576 }
1577 b2.high.x = b2.low.x = p2->p[0].x;
1578 b2.high.y = b2.low.y = p2->p[0].y;
1579 for (i = 1; i < p2->npts; i++)
1580 {
1581 b2.high.x = Max(p2->p[i].x, b2.high.x);
1582 b2.high.y = Max(p2->p[i].y, b2.high.y);
1583 b2.low.x = Min(p2->p[i].x, b2.low.x);
1584 b2.low.y = Min(p2->p[i].y, b2.low.y);
1585 }
1586 if (!box_ov(&b1, &b2))
1587 PG_RETURN_BOOL(false);
1588
1589 /* pairwise check lseg intersections */
1590 for (i = 0; i < p1->npts; i++)
1591 {
1592 int iprev;
1593
1594 if (i > 0)
1595 iprev = i - 1;
1596 else
1597 {
1598 if (!p1->closed)
1599 continue;
1600 iprev = p1->npts - 1; /* include the closure segment */
1601 }
1602
1603 for (j = 0; j < p2->npts; j++)
1604 {
1605 int jprev;
1606
1607 if (j > 0)
1608 jprev = j - 1;
1609 else
1610 {
1611 if (!p2->closed)
1612 continue;
1613 jprev = p2->npts - 1; /* include the closure segment */
1614 }
1615
1616 statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1617 statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1618 if (lseg_intersect_internal(&seg1, &seg2))
1619 PG_RETURN_BOOL(true);
1620 }
1621 }
1622
1623 /* if we dropped through, no two segs intersected */
1624 PG_RETURN_BOOL(false);
1625 }
1626
1627 /* path_distance()
1628 * This essentially does a cartesian product of the lsegs in the
1629 * two paths, and finds the min distance between any two lsegs
1630 */
1631 Datum
path_distance(PG_FUNCTION_ARGS)1632 path_distance(PG_FUNCTION_ARGS)
1633 {
1634 PATH *p1 = PG_GETARG_PATH_P(0);
1635 PATH *p2 = PG_GETARG_PATH_P(1);
1636 float8 min = 0.0; /* initialize to keep compiler quiet */
1637 bool have_min = false;
1638 float8 tmp;
1639 int i,
1640 j;
1641 LSEG seg1,
1642 seg2;
1643
1644 for (i = 0; i < p1->npts; i++)
1645 {
1646 int iprev;
1647
1648 if (i > 0)
1649 iprev = i - 1;
1650 else
1651 {
1652 if (!p1->closed)
1653 continue;
1654 iprev = p1->npts - 1; /* include the closure segment */
1655 }
1656
1657 for (j = 0; j < p2->npts; j++)
1658 {
1659 int jprev;
1660
1661 if (j > 0)
1662 jprev = j - 1;
1663 else
1664 {
1665 if (!p2->closed)
1666 continue;
1667 jprev = p2->npts - 1; /* include the closure segment */
1668 }
1669
1670 statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1671 statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1672
1673 tmp = DatumGetFloat8(DirectFunctionCall2(lseg_distance,
1674 LsegPGetDatum(&seg1),
1675 LsegPGetDatum(&seg2)));
1676 if (!have_min || tmp < min)
1677 {
1678 min = tmp;
1679 have_min = true;
1680 }
1681 }
1682 }
1683
1684 if (!have_min)
1685 PG_RETURN_NULL();
1686
1687 PG_RETURN_FLOAT8(min);
1688 }
1689
1690
1691 /*----------------------------------------------------------
1692 * "Arithmetic" operations.
1693 *---------------------------------------------------------*/
1694
1695 Datum
path_length(PG_FUNCTION_ARGS)1696 path_length(PG_FUNCTION_ARGS)
1697 {
1698 PATH *path = PG_GETARG_PATH_P(0);
1699 float8 result = 0.0;
1700 int i;
1701
1702 for (i = 0; i < path->npts; i++)
1703 {
1704 int iprev;
1705
1706 if (i > 0)
1707 iprev = i - 1;
1708 else
1709 {
1710 if (!path->closed)
1711 continue;
1712 iprev = path->npts - 1; /* include the closure segment */
1713 }
1714
1715 result += point_dt(&path->p[iprev], &path->p[i]);
1716 }
1717
1718 PG_RETURN_FLOAT8(result);
1719 }
1720
1721 /***********************************************************************
1722 **
1723 ** Routines for 2D points.
1724 **
1725 ***********************************************************************/
1726
1727 /*----------------------------------------------------------
1728 * String to point, point to string conversion.
1729 * External format:
1730 * "(x,y)"
1731 * "x,y"
1732 *---------------------------------------------------------*/
1733
1734 Datum
point_in(PG_FUNCTION_ARGS)1735 point_in(PG_FUNCTION_ARGS)
1736 {
1737 char *str = PG_GETARG_CSTRING(0);
1738 Point *point = (Point *) palloc(sizeof(Point));
1739
1740 pair_decode(str, &point->x, &point->y, NULL, "point", str);
1741 PG_RETURN_POINT_P(point);
1742 }
1743
1744 Datum
point_out(PG_FUNCTION_ARGS)1745 point_out(PG_FUNCTION_ARGS)
1746 {
1747 Point *pt = PG_GETARG_POINT_P(0);
1748
1749 PG_RETURN_CSTRING(path_encode(PATH_NONE, 1, pt));
1750 }
1751
1752 /*
1753 * point_recv - converts external binary format to point
1754 */
1755 Datum
point_recv(PG_FUNCTION_ARGS)1756 point_recv(PG_FUNCTION_ARGS)
1757 {
1758 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1759 Point *point;
1760
1761 point = (Point *) palloc(sizeof(Point));
1762 point->x = pq_getmsgfloat8(buf);
1763 point->y = pq_getmsgfloat8(buf);
1764 PG_RETURN_POINT_P(point);
1765 }
1766
1767 /*
1768 * point_send - converts point to binary format
1769 */
1770 Datum
point_send(PG_FUNCTION_ARGS)1771 point_send(PG_FUNCTION_ARGS)
1772 {
1773 Point *pt = PG_GETARG_POINT_P(0);
1774 StringInfoData buf;
1775
1776 pq_begintypsend(&buf);
1777 pq_sendfloat8(&buf, pt->x);
1778 pq_sendfloat8(&buf, pt->y);
1779 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
1780 }
1781
1782
1783 static Point *
point_construct(double x,double y)1784 point_construct(double x, double y)
1785 {
1786 Point *result = (Point *) palloc(sizeof(Point));
1787
1788 result->x = x;
1789 result->y = y;
1790 return result;
1791 }
1792
1793
1794 static Point *
point_copy(Point * pt)1795 point_copy(Point *pt)
1796 {
1797 Point *result;
1798
1799 if (!PointerIsValid(pt))
1800 return NULL;
1801
1802 result = (Point *) palloc(sizeof(Point));
1803
1804 result->x = pt->x;
1805 result->y = pt->y;
1806 return result;
1807 }
1808
1809
1810 /*----------------------------------------------------------
1811 * Relational operators for Points.
1812 * Since we do have a sense of coordinates being
1813 * "equal" to a given accuracy (point_vert, point_horiz),
1814 * the other ops must preserve that sense. This means
1815 * that results may, strictly speaking, be a lie (unless
1816 * EPSILON = 0.0).
1817 *---------------------------------------------------------*/
1818
1819 Datum
point_left(PG_FUNCTION_ARGS)1820 point_left(PG_FUNCTION_ARGS)
1821 {
1822 Point *pt1 = PG_GETARG_POINT_P(0);
1823 Point *pt2 = PG_GETARG_POINT_P(1);
1824
1825 PG_RETURN_BOOL(FPlt(pt1->x, pt2->x));
1826 }
1827
1828 Datum
point_right(PG_FUNCTION_ARGS)1829 point_right(PG_FUNCTION_ARGS)
1830 {
1831 Point *pt1 = PG_GETARG_POINT_P(0);
1832 Point *pt2 = PG_GETARG_POINT_P(1);
1833
1834 PG_RETURN_BOOL(FPgt(pt1->x, pt2->x));
1835 }
1836
1837 Datum
point_above(PG_FUNCTION_ARGS)1838 point_above(PG_FUNCTION_ARGS)
1839 {
1840 Point *pt1 = PG_GETARG_POINT_P(0);
1841 Point *pt2 = PG_GETARG_POINT_P(1);
1842
1843 PG_RETURN_BOOL(FPgt(pt1->y, pt2->y));
1844 }
1845
1846 Datum
point_below(PG_FUNCTION_ARGS)1847 point_below(PG_FUNCTION_ARGS)
1848 {
1849 Point *pt1 = PG_GETARG_POINT_P(0);
1850 Point *pt2 = PG_GETARG_POINT_P(1);
1851
1852 PG_RETURN_BOOL(FPlt(pt1->y, pt2->y));
1853 }
1854
1855 Datum
point_vert(PG_FUNCTION_ARGS)1856 point_vert(PG_FUNCTION_ARGS)
1857 {
1858 Point *pt1 = PG_GETARG_POINT_P(0);
1859 Point *pt2 = PG_GETARG_POINT_P(1);
1860
1861 PG_RETURN_BOOL(FPeq(pt1->x, pt2->x));
1862 }
1863
1864 Datum
point_horiz(PG_FUNCTION_ARGS)1865 point_horiz(PG_FUNCTION_ARGS)
1866 {
1867 Point *pt1 = PG_GETARG_POINT_P(0);
1868 Point *pt2 = PG_GETARG_POINT_P(1);
1869
1870 PG_RETURN_BOOL(FPeq(pt1->y, pt2->y));
1871 }
1872
1873 Datum
point_eq(PG_FUNCTION_ARGS)1874 point_eq(PG_FUNCTION_ARGS)
1875 {
1876 Point *pt1 = PG_GETARG_POINT_P(0);
1877 Point *pt2 = PG_GETARG_POINT_P(1);
1878
1879 PG_RETURN_BOOL(FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y));
1880 }
1881
1882 Datum
point_ne(PG_FUNCTION_ARGS)1883 point_ne(PG_FUNCTION_ARGS)
1884 {
1885 Point *pt1 = PG_GETARG_POINT_P(0);
1886 Point *pt2 = PG_GETARG_POINT_P(1);
1887
1888 PG_RETURN_BOOL(FPne(pt1->x, pt2->x) || FPne(pt1->y, pt2->y));
1889 }
1890
1891 /*----------------------------------------------------------
1892 * "Arithmetic" operators on points.
1893 *---------------------------------------------------------*/
1894
1895 Datum
point_distance(PG_FUNCTION_ARGS)1896 point_distance(PG_FUNCTION_ARGS)
1897 {
1898 Point *pt1 = PG_GETARG_POINT_P(0);
1899 Point *pt2 = PG_GETARG_POINT_P(1);
1900
1901 PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
1902 }
1903
1904 double
point_dt(Point * pt1,Point * pt2)1905 point_dt(Point *pt1, Point *pt2)
1906 {
1907 #ifdef GEODEBUG
1908 printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
1909 pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
1910 #endif
1911 return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
1912 }
1913
1914 Datum
point_slope(PG_FUNCTION_ARGS)1915 point_slope(PG_FUNCTION_ARGS)
1916 {
1917 Point *pt1 = PG_GETARG_POINT_P(0);
1918 Point *pt2 = PG_GETARG_POINT_P(1);
1919
1920 PG_RETURN_FLOAT8(point_sl(pt1, pt2));
1921 }
1922
1923
1924 double
point_sl(Point * pt1,Point * pt2)1925 point_sl(Point *pt1, Point *pt2)
1926 {
1927 return (FPeq(pt1->x, pt2->x)
1928 ? (double) DBL_MAX
1929 : (pt1->y - pt2->y) / (pt1->x - pt2->x));
1930 }
1931
1932
1933 /***********************************************************************
1934 **
1935 ** Routines for 2D line segments.
1936 **
1937 ***********************************************************************/
1938
1939 /*----------------------------------------------------------
1940 * String to lseg, lseg to string conversion.
1941 * External forms: "[(x1, y1), (x2, y2)]"
1942 * "(x1, y1), (x2, y2)"
1943 * "x1, y1, x2, y2"
1944 * closed form ok "((x1, y1), (x2, y2))"
1945 * (old form) "(x1, y1, x2, y2)"
1946 *---------------------------------------------------------*/
1947
1948 Datum
lseg_in(PG_FUNCTION_ARGS)1949 lseg_in(PG_FUNCTION_ARGS)
1950 {
1951 char *str = PG_GETARG_CSTRING(0);
1952 LSEG *lseg = (LSEG *) palloc(sizeof(LSEG));
1953 bool isopen;
1954
1955 path_decode(str, true, 2, &(lseg->p[0]), &isopen, NULL, "lseg", str);
1956 PG_RETURN_LSEG_P(lseg);
1957 }
1958
1959
1960 Datum
lseg_out(PG_FUNCTION_ARGS)1961 lseg_out(PG_FUNCTION_ARGS)
1962 {
1963 LSEG *ls = PG_GETARG_LSEG_P(0);
1964
1965 PG_RETURN_CSTRING(path_encode(PATH_OPEN, 2, (Point *) &(ls->p[0])));
1966 }
1967
1968 /*
1969 * lseg_recv - converts external binary format to lseg
1970 */
1971 Datum
lseg_recv(PG_FUNCTION_ARGS)1972 lseg_recv(PG_FUNCTION_ARGS)
1973 {
1974 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
1975 LSEG *lseg;
1976
1977 lseg = (LSEG *) palloc(sizeof(LSEG));
1978
1979 lseg->p[0].x = pq_getmsgfloat8(buf);
1980 lseg->p[0].y = pq_getmsgfloat8(buf);
1981 lseg->p[1].x = pq_getmsgfloat8(buf);
1982 lseg->p[1].y = pq_getmsgfloat8(buf);
1983
1984 PG_RETURN_LSEG_P(lseg);
1985 }
1986
1987 /*
1988 * lseg_send - converts lseg to binary format
1989 */
1990 Datum
lseg_send(PG_FUNCTION_ARGS)1991 lseg_send(PG_FUNCTION_ARGS)
1992 {
1993 LSEG *ls = PG_GETARG_LSEG_P(0);
1994 StringInfoData buf;
1995
1996 pq_begintypsend(&buf);
1997 pq_sendfloat8(&buf, ls->p[0].x);
1998 pq_sendfloat8(&buf, ls->p[0].y);
1999 pq_sendfloat8(&buf, ls->p[1].x);
2000 pq_sendfloat8(&buf, ls->p[1].y);
2001 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
2002 }
2003
2004
2005 /* lseg_construct -
2006 * form a LSEG from two Points.
2007 */
2008 Datum
lseg_construct(PG_FUNCTION_ARGS)2009 lseg_construct(PG_FUNCTION_ARGS)
2010 {
2011 Point *pt1 = PG_GETARG_POINT_P(0);
2012 Point *pt2 = PG_GETARG_POINT_P(1);
2013 LSEG *result = (LSEG *) palloc(sizeof(LSEG));
2014
2015 result->p[0].x = pt1->x;
2016 result->p[0].y = pt1->y;
2017 result->p[1].x = pt2->x;
2018 result->p[1].y = pt2->y;
2019
2020 PG_RETURN_LSEG_P(result);
2021 }
2022
2023 /* like lseg_construct, but assume space already allocated */
2024 static void
statlseg_construct(LSEG * lseg,Point * pt1,Point * pt2)2025 statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2)
2026 {
2027 lseg->p[0].x = pt1->x;
2028 lseg->p[0].y = pt1->y;
2029 lseg->p[1].x = pt2->x;
2030 lseg->p[1].y = pt2->y;
2031 }
2032
2033 Datum
lseg_length(PG_FUNCTION_ARGS)2034 lseg_length(PG_FUNCTION_ARGS)
2035 {
2036 LSEG *lseg = PG_GETARG_LSEG_P(0);
2037
2038 PG_RETURN_FLOAT8(point_dt(&lseg->p[0], &lseg->p[1]));
2039 }
2040
2041 /*----------------------------------------------------------
2042 * Relative position routines.
2043 *---------------------------------------------------------*/
2044
2045 /*
2046 ** find intersection of the two lines, and see if it falls on
2047 ** both segments.
2048 */
2049 Datum
lseg_intersect(PG_FUNCTION_ARGS)2050 lseg_intersect(PG_FUNCTION_ARGS)
2051 {
2052 LSEG *l1 = PG_GETARG_LSEG_P(0);
2053 LSEG *l2 = PG_GETARG_LSEG_P(1);
2054
2055 PG_RETURN_BOOL(lseg_intersect_internal(l1, l2));
2056 }
2057
2058 static bool
lseg_intersect_internal(LSEG * l1,LSEG * l2)2059 lseg_intersect_internal(LSEG *l1, LSEG *l2)
2060 {
2061 LINE ln;
2062 Point *interpt;
2063 bool retval;
2064
2065 line_construct_pts(&ln, &l2->p[0], &l2->p[1]);
2066 interpt = interpt_sl(l1, &ln);
2067
2068 if (interpt != NULL && on_ps_internal(interpt, l2))
2069 retval = true; /* interpt on l1 and l2 */
2070 else
2071 retval = false;
2072 return retval;
2073 }
2074
2075 Datum
lseg_parallel(PG_FUNCTION_ARGS)2076 lseg_parallel(PG_FUNCTION_ARGS)
2077 {
2078 LSEG *l1 = PG_GETARG_LSEG_P(0);
2079 LSEG *l2 = PG_GETARG_LSEG_P(1);
2080
2081 PG_RETURN_BOOL(FPeq(point_sl(&l1->p[0], &l1->p[1]),
2082 point_sl(&l2->p[0], &l2->p[1])));
2083 }
2084
2085 /* lseg_perp()
2086 * Determine if two line segments are perpendicular.
2087 *
2088 * This code did not get the correct answer for
2089 * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
2090 * So, modified it to check explicitly for slope of vertical line
2091 * returned by point_sl() and the results seem better.
2092 * - thomas 1998-01-31
2093 */
2094 Datum
lseg_perp(PG_FUNCTION_ARGS)2095 lseg_perp(PG_FUNCTION_ARGS)
2096 {
2097 LSEG *l1 = PG_GETARG_LSEG_P(0);
2098 LSEG *l2 = PG_GETARG_LSEG_P(1);
2099 double m1,
2100 m2;
2101
2102 m1 = point_sl(&(l1->p[0]), &(l1->p[1]));
2103 m2 = point_sl(&(l2->p[0]), &(l2->p[1]));
2104
2105 #ifdef GEODEBUG
2106 printf("lseg_perp- slopes are %g and %g\n", m1, m2);
2107 #endif
2108 if (FPzero(m1))
2109 PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
2110 else if (FPzero(m2))
2111 PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
2112
2113 PG_RETURN_BOOL(FPeq(m1 / m2, -1.0));
2114 }
2115
2116 Datum
lseg_vertical(PG_FUNCTION_ARGS)2117 lseg_vertical(PG_FUNCTION_ARGS)
2118 {
2119 LSEG *lseg = PG_GETARG_LSEG_P(0);
2120
2121 PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
2122 }
2123
2124 Datum
lseg_horizontal(PG_FUNCTION_ARGS)2125 lseg_horizontal(PG_FUNCTION_ARGS)
2126 {
2127 LSEG *lseg = PG_GETARG_LSEG_P(0);
2128
2129 PG_RETURN_BOOL(FPeq(lseg->p[0].y, lseg->p[1].y));
2130 }
2131
2132
2133 Datum
lseg_eq(PG_FUNCTION_ARGS)2134 lseg_eq(PG_FUNCTION_ARGS)
2135 {
2136 LSEG *l1 = PG_GETARG_LSEG_P(0);
2137 LSEG *l2 = PG_GETARG_LSEG_P(1);
2138
2139 PG_RETURN_BOOL(FPeq(l1->p[0].x, l2->p[0].x) &&
2140 FPeq(l1->p[0].y, l2->p[0].y) &&
2141 FPeq(l1->p[1].x, l2->p[1].x) &&
2142 FPeq(l1->p[1].y, l2->p[1].y));
2143 }
2144
2145 Datum
lseg_ne(PG_FUNCTION_ARGS)2146 lseg_ne(PG_FUNCTION_ARGS)
2147 {
2148 LSEG *l1 = PG_GETARG_LSEG_P(0);
2149 LSEG *l2 = PG_GETARG_LSEG_P(1);
2150
2151 PG_RETURN_BOOL(!FPeq(l1->p[0].x, l2->p[0].x) ||
2152 !FPeq(l1->p[0].y, l2->p[0].y) ||
2153 !FPeq(l1->p[1].x, l2->p[1].x) ||
2154 !FPeq(l1->p[1].y, l2->p[1].y));
2155 }
2156
2157 Datum
lseg_lt(PG_FUNCTION_ARGS)2158 lseg_lt(PG_FUNCTION_ARGS)
2159 {
2160 LSEG *l1 = PG_GETARG_LSEG_P(0);
2161 LSEG *l2 = PG_GETARG_LSEG_P(1);
2162
2163 PG_RETURN_BOOL(FPlt(point_dt(&l1->p[0], &l1->p[1]),
2164 point_dt(&l2->p[0], &l2->p[1])));
2165 }
2166
2167 Datum
lseg_le(PG_FUNCTION_ARGS)2168 lseg_le(PG_FUNCTION_ARGS)
2169 {
2170 LSEG *l1 = PG_GETARG_LSEG_P(0);
2171 LSEG *l2 = PG_GETARG_LSEG_P(1);
2172
2173 PG_RETURN_BOOL(FPle(point_dt(&l1->p[0], &l1->p[1]),
2174 point_dt(&l2->p[0], &l2->p[1])));
2175 }
2176
2177 Datum
lseg_gt(PG_FUNCTION_ARGS)2178 lseg_gt(PG_FUNCTION_ARGS)
2179 {
2180 LSEG *l1 = PG_GETARG_LSEG_P(0);
2181 LSEG *l2 = PG_GETARG_LSEG_P(1);
2182
2183 PG_RETURN_BOOL(FPgt(point_dt(&l1->p[0], &l1->p[1]),
2184 point_dt(&l2->p[0], &l2->p[1])));
2185 }
2186
2187 Datum
lseg_ge(PG_FUNCTION_ARGS)2188 lseg_ge(PG_FUNCTION_ARGS)
2189 {
2190 LSEG *l1 = PG_GETARG_LSEG_P(0);
2191 LSEG *l2 = PG_GETARG_LSEG_P(1);
2192
2193 PG_RETURN_BOOL(FPge(point_dt(&l1->p[0], &l1->p[1]),
2194 point_dt(&l2->p[0], &l2->p[1])));
2195 }
2196
2197
2198 /*----------------------------------------------------------
2199 * Line arithmetic routines.
2200 *---------------------------------------------------------*/
2201
2202 /* lseg_distance -
2203 * If two segments don't intersect, then the closest
2204 * point will be from one of the endpoints to the other
2205 * segment.
2206 */
2207 Datum
lseg_distance(PG_FUNCTION_ARGS)2208 lseg_distance(PG_FUNCTION_ARGS)
2209 {
2210 LSEG *l1 = PG_GETARG_LSEG_P(0);
2211 LSEG *l2 = PG_GETARG_LSEG_P(1);
2212
2213 PG_RETURN_FLOAT8(lseg_dt(l1, l2));
2214 }
2215
2216 /* lseg_dt()
2217 * Distance between two line segments.
2218 * Must check both sets of endpoints to ensure minimum distance is found.
2219 * - thomas 1998-02-01
2220 */
2221 static double
lseg_dt(LSEG * l1,LSEG * l2)2222 lseg_dt(LSEG *l1, LSEG *l2)
2223 {
2224 double result,
2225 d;
2226
2227 if (lseg_intersect_internal(l1, l2))
2228 return 0.0;
2229
2230 d = dist_ps_internal(&l1->p[0], l2);
2231 result = d;
2232 d = dist_ps_internal(&l1->p[1], l2);
2233 result = Min(result, d);
2234 d = dist_ps_internal(&l2->p[0], l1);
2235 result = Min(result, d);
2236 d = dist_ps_internal(&l2->p[1], l1);
2237 result = Min(result, d);
2238
2239 return result;
2240 }
2241
2242
2243 Datum
lseg_center(PG_FUNCTION_ARGS)2244 lseg_center(PG_FUNCTION_ARGS)
2245 {
2246 LSEG *lseg = PG_GETARG_LSEG_P(0);
2247 Point *result;
2248
2249 result = (Point *) palloc(sizeof(Point));
2250
2251 result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0;
2252 result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0;
2253
2254 PG_RETURN_POINT_P(result);
2255 }
2256
2257 static Point *
lseg_interpt_internal(LSEG * l1,LSEG * l2)2258 lseg_interpt_internal(LSEG *l1, LSEG *l2)
2259 {
2260 Point *result;
2261 LINE tmp1,
2262 tmp2;
2263
2264 /*
2265 * Find the intersection of the appropriate lines, if any.
2266 */
2267 line_construct_pts(&tmp1, &l1->p[0], &l1->p[1]);
2268 line_construct_pts(&tmp2, &l2->p[0], &l2->p[1]);
2269 result = line_interpt_internal(&tmp1, &tmp2);
2270 if (!PointerIsValid(result))
2271 return NULL;
2272
2273 /*
2274 * If the line intersection point isn't within l1 (or equivalently l2),
2275 * there is no valid segment intersection point at all.
2276 */
2277 if (!on_ps_internal(result, l1) ||
2278 !on_ps_internal(result, l2))
2279 {
2280 pfree(result);
2281 return NULL;
2282 }
2283
2284 /*
2285 * If there is an intersection, then check explicitly for matching
2286 * endpoints since there may be rounding effects with annoying lsb
2287 * residue. - tgl 1997-07-09
2288 */
2289 if ((FPeq(l1->p[0].x, l2->p[0].x) && FPeq(l1->p[0].y, l2->p[0].y)) ||
2290 (FPeq(l1->p[0].x, l2->p[1].x) && FPeq(l1->p[0].y, l2->p[1].y)))
2291 {
2292 result->x = l1->p[0].x;
2293 result->y = l1->p[0].y;
2294 }
2295 else if ((FPeq(l1->p[1].x, l2->p[0].x) && FPeq(l1->p[1].y, l2->p[0].y)) ||
2296 (FPeq(l1->p[1].x, l2->p[1].x) && FPeq(l1->p[1].y, l2->p[1].y)))
2297 {
2298 result->x = l1->p[1].x;
2299 result->y = l1->p[1].y;
2300 }
2301
2302 return result;
2303 }
2304
2305 /* lseg_interpt -
2306 * Find the intersection point of two segments (if any).
2307 */
2308 Datum
lseg_interpt(PG_FUNCTION_ARGS)2309 lseg_interpt(PG_FUNCTION_ARGS)
2310 {
2311 LSEG *l1 = PG_GETARG_LSEG_P(0);
2312 LSEG *l2 = PG_GETARG_LSEG_P(1);
2313 Point *result;
2314
2315 result = lseg_interpt_internal(l1, l2);
2316 if (!PointerIsValid(result))
2317 PG_RETURN_NULL();
2318
2319 PG_RETURN_POINT_P(result);
2320 }
2321
2322 /***********************************************************************
2323 **
2324 ** Routines for position comparisons of differently-typed
2325 ** 2D objects.
2326 **
2327 ***********************************************************************/
2328
2329 /*---------------------------------------------------------------------
2330 * dist_
2331 * Minimum distance from one object to another.
2332 *-------------------------------------------------------------------*/
2333
2334 /*
2335 * Distance from a point to a line
2336 */
2337 Datum
dist_pl(PG_FUNCTION_ARGS)2338 dist_pl(PG_FUNCTION_ARGS)
2339 {
2340 Point *pt = PG_GETARG_POINT_P(0);
2341 LINE *line = PG_GETARG_LINE_P(1);
2342
2343 PG_RETURN_FLOAT8(dist_pl_internal(pt, line));
2344 }
2345
2346 static double
dist_pl_internal(Point * pt,LINE * line)2347 dist_pl_internal(Point *pt, LINE *line)
2348 {
2349 return fabs((line->A * pt->x + line->B * pt->y + line->C) /
2350 HYPOT(line->A, line->B));
2351 }
2352
2353 /*
2354 * Distance from a point to a lseg
2355 */
2356 Datum
dist_ps(PG_FUNCTION_ARGS)2357 dist_ps(PG_FUNCTION_ARGS)
2358 {
2359 Point *pt = PG_GETARG_POINT_P(0);
2360 LSEG *lseg = PG_GETARG_LSEG_P(1);
2361
2362 PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg));
2363 }
2364
2365 static double
dist_ps_internal(Point * pt,LSEG * lseg)2366 dist_ps_internal(Point *pt, LSEG *lseg)
2367 {
2368 double m; /* slope of perp. */
2369 LINE *ln;
2370 double result,
2371 tmpdist;
2372 Point *ip;
2373
2374 /*
2375 * Construct a line perpendicular to the input segment and through the
2376 * input point
2377 */
2378 if (lseg->p[1].x == lseg->p[0].x)
2379 m = 0;
2380 else if (lseg->p[1].y == lseg->p[0].y)
2381 m = (double) DBL_MAX; /* slope is infinite */
2382 else
2383 m = (lseg->p[0].x - lseg->p[1].x) / (lseg->p[1].y - lseg->p[0].y);
2384 ln = line_construct_pm(pt, m);
2385
2386 #ifdef GEODEBUG
2387 printf("dist_ps- line is A=%g B=%g C=%g from (point) slope (%f,%f) %g\n",
2388 ln->A, ln->B, ln->C, pt->x, pt->y, m);
2389 #endif
2390
2391 /*
2392 * Calculate distance to the line segment or to the nearest endpoint of
2393 * the segment.
2394 */
2395
2396 /* intersection is on the line segment? */
2397 if ((ip = interpt_sl(lseg, ln)) != NULL)
2398 {
2399 /* yes, so use distance to the intersection point */
2400 result = point_dt(pt, ip);
2401 #ifdef GEODEBUG
2402 printf("dist_ps- distance is %f to intersection point is (%f,%f)\n",
2403 result, ip->x, ip->y);
2404 #endif
2405 }
2406 else
2407 {
2408 /* no, so use distance to the nearer endpoint */
2409 result = point_dt(pt, &lseg->p[0]);
2410 tmpdist = point_dt(pt, &lseg->p[1]);
2411 if (tmpdist < result)
2412 result = tmpdist;
2413 }
2414
2415 return result;
2416 }
2417
2418 /*
2419 * Distance from a point to a path
2420 */
2421 Datum
dist_ppath(PG_FUNCTION_ARGS)2422 dist_ppath(PG_FUNCTION_ARGS)
2423 {
2424 Point *pt = PG_GETARG_POINT_P(0);
2425 PATH *path = PG_GETARG_PATH_P(1);
2426 float8 result = 0.0; /* keep compiler quiet */
2427 bool have_min = false;
2428 float8 tmp;
2429 int i;
2430 LSEG lseg;
2431
2432 switch (path->npts)
2433 {
2434 case 0:
2435 /* no points in path? then result is undefined... */
2436 PG_RETURN_NULL();
2437 case 1:
2438 /* one point in path? then get distance between two points... */
2439 result = point_dt(pt, &path->p[0]);
2440 break;
2441 default:
2442 /* make sure the path makes sense... */
2443 Assert(path->npts > 1);
2444
2445 /*
2446 * the distance from a point to a path is the smallest distance
2447 * from the point to any of its constituent segments.
2448 */
2449 for (i = 0; i < path->npts; i++)
2450 {
2451 int iprev;
2452
2453 if (i > 0)
2454 iprev = i - 1;
2455 else
2456 {
2457 if (!path->closed)
2458 continue;
2459 iprev = path->npts - 1; /* include the closure segment */
2460 }
2461
2462 statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
2463 tmp = dist_ps_internal(pt, &lseg);
2464 if (!have_min || tmp < result)
2465 {
2466 result = tmp;
2467 have_min = true;
2468 }
2469 }
2470 break;
2471 }
2472 PG_RETURN_FLOAT8(result);
2473 }
2474
2475 /*
2476 * Distance from a point to a box
2477 */
2478 Datum
dist_pb(PG_FUNCTION_ARGS)2479 dist_pb(PG_FUNCTION_ARGS)
2480 {
2481 Point *pt = PG_GETARG_POINT_P(0);
2482 BOX *box = PG_GETARG_BOX_P(1);
2483 float8 result;
2484 Point *near;
2485
2486 near = DatumGetPointP(DirectFunctionCall2(close_pb,
2487 PointPGetDatum(pt),
2488 BoxPGetDatum(box)));
2489 result = point_dt(near, pt);
2490
2491 PG_RETURN_FLOAT8(result);
2492 }
2493
2494 /*
2495 * Distance from a lseg to a line
2496 */
2497 Datum
dist_sl(PG_FUNCTION_ARGS)2498 dist_sl(PG_FUNCTION_ARGS)
2499 {
2500 LSEG *lseg = PG_GETARG_LSEG_P(0);
2501 LINE *line = PG_GETARG_LINE_P(1);
2502 float8 result,
2503 d2;
2504
2505 if (has_interpt_sl(lseg, line))
2506 result = 0.0;
2507 else
2508 {
2509 result = dist_pl_internal(&lseg->p[0], line);
2510 d2 = dist_pl_internal(&lseg->p[1], line);
2511 /* XXX shouldn't we take the min not max? */
2512 if (d2 > result)
2513 result = d2;
2514 }
2515
2516 PG_RETURN_FLOAT8(result);
2517 }
2518
2519 /*
2520 * Distance from a lseg to a box
2521 */
2522 Datum
dist_sb(PG_FUNCTION_ARGS)2523 dist_sb(PG_FUNCTION_ARGS)
2524 {
2525 LSEG *lseg = PG_GETARG_LSEG_P(0);
2526 BOX *box = PG_GETARG_BOX_P(1);
2527 Point *tmp;
2528 Datum result;
2529
2530 tmp = DatumGetPointP(DirectFunctionCall2(close_sb,
2531 LsegPGetDatum(lseg),
2532 BoxPGetDatum(box)));
2533 result = DirectFunctionCall2(dist_pb,
2534 PointPGetDatum(tmp),
2535 BoxPGetDatum(box));
2536
2537 PG_RETURN_DATUM(result);
2538 }
2539
2540 /*
2541 * Distance from a line to a box
2542 */
2543 Datum
dist_lb(PG_FUNCTION_ARGS)2544 dist_lb(PG_FUNCTION_ARGS)
2545 {
2546 #ifdef NOT_USED
2547 LINE *line = PG_GETARG_LINE_P(0);
2548 BOX *box = PG_GETARG_BOX_P(1);
2549 #endif
2550
2551 /* need to think about this one for a while */
2552 ereport(ERROR,
2553 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2554 errmsg("function \"dist_lb\" not implemented")));
2555
2556 PG_RETURN_NULL();
2557 }
2558
2559 /*
2560 * Distance from a circle to a polygon
2561 */
2562 Datum
dist_cpoly(PG_FUNCTION_ARGS)2563 dist_cpoly(PG_FUNCTION_ARGS)
2564 {
2565 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
2566 POLYGON *poly = PG_GETARG_POLYGON_P(1);
2567 float8 result;
2568
2569 /* calculate distance to center, and subtract radius */
2570 result = dist_ppoly_internal(&circle->center, poly);
2571
2572 result -= circle->radius;
2573 if (result < 0)
2574 result = 0;
2575
2576 PG_RETURN_FLOAT8(result);
2577 }
2578
2579 /*
2580 * Distance from a point to a polygon
2581 */
2582 Datum
dist_ppoly(PG_FUNCTION_ARGS)2583 dist_ppoly(PG_FUNCTION_ARGS)
2584 {
2585 Point *point = PG_GETARG_POINT_P(0);
2586 POLYGON *poly = PG_GETARG_POLYGON_P(1);
2587 float8 result;
2588
2589 result = dist_ppoly_internal(point, poly);
2590
2591 PG_RETURN_FLOAT8(result);
2592 }
2593
2594 Datum
dist_polyp(PG_FUNCTION_ARGS)2595 dist_polyp(PG_FUNCTION_ARGS)
2596 {
2597 POLYGON *poly = PG_GETARG_POLYGON_P(0);
2598 Point *point = PG_GETARG_POINT_P(1);
2599 float8 result;
2600
2601 result = dist_ppoly_internal(point, poly);
2602
2603 PG_RETURN_FLOAT8(result);
2604 }
2605
2606 static double
dist_ppoly_internal(Point * pt,POLYGON * poly)2607 dist_ppoly_internal(Point *pt, POLYGON *poly)
2608 {
2609 float8 result;
2610 float8 d;
2611 int i;
2612 LSEG seg;
2613
2614 if (point_inside(pt, poly->npts, poly->p) != 0)
2615 {
2616 #ifdef GEODEBUG
2617 printf("dist_ppoly_internal- point inside of polygon\n");
2618 #endif
2619 return 0.0;
2620 }
2621
2622 /* initialize distance with segment between first and last points */
2623 seg.p[0].x = poly->p[0].x;
2624 seg.p[0].y = poly->p[0].y;
2625 seg.p[1].x = poly->p[poly->npts - 1].x;
2626 seg.p[1].y = poly->p[poly->npts - 1].y;
2627 result = dist_ps_internal(pt, &seg);
2628 #ifdef GEODEBUG
2629 printf("dist_ppoly_internal- segment 0/n distance is %f\n", result);
2630 #endif
2631
2632 /* check distances for other segments */
2633 for (i = 0; (i < poly->npts - 1); i++)
2634 {
2635 seg.p[0].x = poly->p[i].x;
2636 seg.p[0].y = poly->p[i].y;
2637 seg.p[1].x = poly->p[i + 1].x;
2638 seg.p[1].y = poly->p[i + 1].y;
2639 d = dist_ps_internal(pt, &seg);
2640 #ifdef GEODEBUG
2641 printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d);
2642 #endif
2643 if (d < result)
2644 result = d;
2645 }
2646
2647 return result;
2648 }
2649
2650
2651 /*---------------------------------------------------------------------
2652 * interpt_
2653 * Intersection point of objects.
2654 * We choose to ignore the "point" of intersection between
2655 * lines and boxes, since there are typically two.
2656 *-------------------------------------------------------------------*/
2657
2658 /* Get intersection point of lseg and line; returns NULL if no intersection */
2659 static Point *
interpt_sl(LSEG * lseg,LINE * line)2660 interpt_sl(LSEG *lseg, LINE *line)
2661 {
2662 LINE tmp;
2663 Point *p;
2664
2665 line_construct_pts(&tmp, &lseg->p[0], &lseg->p[1]);
2666 p = line_interpt_internal(&tmp, line);
2667 #ifdef GEODEBUG
2668 printf("interpt_sl- segment is (%.*g %.*g) (%.*g %.*g)\n",
2669 DBL_DIG, lseg->p[0].x, DBL_DIG, lseg->p[0].y, DBL_DIG, lseg->p[1].x, DBL_DIG, lseg->p[1].y);
2670 printf("interpt_sl- segment becomes line A=%.*g B=%.*g C=%.*g\n",
2671 DBL_DIG, tmp.A, DBL_DIG, tmp.B, DBL_DIG, tmp.C);
2672 #endif
2673 if (PointerIsValid(p))
2674 {
2675 #ifdef GEODEBUG
2676 printf("interpt_sl- intersection point is (%.*g %.*g)\n", DBL_DIG, p->x, DBL_DIG, p->y);
2677 #endif
2678 if (on_ps_internal(p, lseg))
2679 {
2680 #ifdef GEODEBUG
2681 printf("interpt_sl- intersection point is on segment\n");
2682 #endif
2683 }
2684 else
2685 p = NULL;
2686 }
2687
2688 return p;
2689 }
2690
2691 /* variant: just indicate if intersection point exists */
2692 static bool
has_interpt_sl(LSEG * lseg,LINE * line)2693 has_interpt_sl(LSEG *lseg, LINE *line)
2694 {
2695 Point *tmp;
2696
2697 tmp = interpt_sl(lseg, line);
2698 if (tmp)
2699 return true;
2700 return false;
2701 }
2702
2703 /*---------------------------------------------------------------------
2704 * close_
2705 * Point of closest proximity between objects.
2706 *-------------------------------------------------------------------*/
2707
2708 /* close_pl -
2709 * The intersection point of a perpendicular of the line
2710 * through the point.
2711 */
2712 Datum
close_pl(PG_FUNCTION_ARGS)2713 close_pl(PG_FUNCTION_ARGS)
2714 {
2715 Point *pt = PG_GETARG_POINT_P(0);
2716 LINE *line = PG_GETARG_LINE_P(1);
2717 Point *result;
2718 LINE *tmp;
2719 double invm;
2720
2721 result = (Point *) palloc(sizeof(Point));
2722
2723 if (FPzero(line->B)) /* vertical? */
2724 {
2725 result->x = line->C;
2726 result->y = pt->y;
2727 PG_RETURN_POINT_P(result);
2728 }
2729 if (FPzero(line->A)) /* horizontal? */
2730 {
2731 result->x = pt->x;
2732 result->y = line->C;
2733 PG_RETURN_POINT_P(result);
2734 }
2735 /* drop a perpendicular and find the intersection point */
2736
2737 /* invert and flip the sign on the slope to get a perpendicular */
2738 invm = line->B / line->A;
2739 tmp = line_construct_pm(pt, invm);
2740 result = line_interpt_internal(tmp, line);
2741 Assert(result != NULL);
2742 PG_RETURN_POINT_P(result);
2743 }
2744
2745
2746 /* close_ps()
2747 * Closest point on line segment to specified point.
2748 * Take the closest endpoint if the point is left, right,
2749 * above, or below the segment, otherwise find the intersection
2750 * point of the segment and its perpendicular through the point.
2751 *
2752 * Some tricky code here, relying on boolean expressions
2753 * evaluating to only zero or one to use as an array index.
2754 * bug fixes by gthaker@atl.lmco.com; May 1, 1998
2755 */
2756 Datum
close_ps(PG_FUNCTION_ARGS)2757 close_ps(PG_FUNCTION_ARGS)
2758 {
2759 Point *pt = PG_GETARG_POINT_P(0);
2760 LSEG *lseg = PG_GETARG_LSEG_P(1);
2761 Point *result = NULL;
2762 LINE *tmp;
2763 double invm;
2764 int xh,
2765 yh;
2766
2767 #ifdef GEODEBUG
2768 printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f lseg(1).x %f lseg(1).y %f\n",
2769 pt->x, pt->y, lseg->p[0].x, lseg->p[0].y,
2770 lseg->p[1].x, lseg->p[1].y);
2771 #endif
2772
2773 /* xh (or yh) is the index of upper x( or y) end point of lseg */
2774 /* !xh (or !yh) is the index of lower x( or y) end point of lseg */
2775 xh = lseg->p[0].x < lseg->p[1].x;
2776 yh = lseg->p[0].y < lseg->p[1].y;
2777
2778 if (FPeq(lseg->p[0].x, lseg->p[1].x)) /* vertical? */
2779 {
2780 #ifdef GEODEBUG
2781 printf("close_ps- segment is vertical\n");
2782 #endif
2783 /* first check if point is below or above the entire lseg. */
2784 if (pt->y < lseg->p[!yh].y)
2785 result = point_copy(&lseg->p[!yh]); /* below the lseg */
2786 else if (pt->y > lseg->p[yh].y)
2787 result = point_copy(&lseg->p[yh]); /* above the lseg */
2788 if (result != NULL)
2789 PG_RETURN_POINT_P(result);
2790
2791 /* point lines along (to left or right) of the vertical lseg. */
2792
2793 result = (Point *) palloc(sizeof(Point));
2794 result->x = lseg->p[0].x;
2795 result->y = pt->y;
2796 PG_RETURN_POINT_P(result);
2797 }
2798 else if (FPeq(lseg->p[0].y, lseg->p[1].y)) /* horizontal? */
2799 {
2800 #ifdef GEODEBUG
2801 printf("close_ps- segment is horizontal\n");
2802 #endif
2803 /* first check if point is left or right of the entire lseg. */
2804 if (pt->x < lseg->p[!xh].x)
2805 result = point_copy(&lseg->p[!xh]); /* left of the lseg */
2806 else if (pt->x > lseg->p[xh].x)
2807 result = point_copy(&lseg->p[xh]); /* right of the lseg */
2808 if (result != NULL)
2809 PG_RETURN_POINT_P(result);
2810
2811 /* point lines along (at top or below) the horiz. lseg. */
2812 result = (Point *) palloc(sizeof(Point));
2813 result->x = pt->x;
2814 result->y = lseg->p[0].y;
2815 PG_RETURN_POINT_P(result);
2816 }
2817
2818 /*
2819 * vert. and horiz. cases are down, now check if the closest point is one
2820 * of the end points or someplace on the lseg.
2821 */
2822
2823 invm = -1.0 / point_sl(&(lseg->p[0]), &(lseg->p[1]));
2824 tmp = line_construct_pm(&lseg->p[!yh], invm); /* lower edge of the
2825 * "band" */
2826 if (pt->y < (tmp->A * pt->x + tmp->C))
2827 { /* we are below the lower edge */
2828 result = point_copy(&lseg->p[!yh]); /* below the lseg, take lower end
2829 * pt */
2830 #ifdef GEODEBUG
2831 printf("close_ps below: tmp A %f B %f C %f\n",
2832 tmp->A, tmp->B, tmp->C);
2833 #endif
2834 PG_RETURN_POINT_P(result);
2835 }
2836 tmp = line_construct_pm(&lseg->p[yh], invm); /* upper edge of the
2837 * "band" */
2838 if (pt->y > (tmp->A * pt->x + tmp->C))
2839 { /* we are below the lower edge */
2840 result = point_copy(&lseg->p[yh]); /* above the lseg, take higher end
2841 * pt */
2842 #ifdef GEODEBUG
2843 printf("close_ps above: tmp A %f B %f C %f\n",
2844 tmp->A, tmp->B, tmp->C);
2845 #endif
2846 PG_RETURN_POINT_P(result);
2847 }
2848
2849 /*
2850 * at this point the "normal" from point will hit lseg. The closest point
2851 * will be somewhere on the lseg
2852 */
2853 tmp = line_construct_pm(pt, invm);
2854 #ifdef GEODEBUG
2855 printf("close_ps- tmp A %f B %f C %f\n",
2856 tmp->A, tmp->B, tmp->C);
2857 #endif
2858 result = interpt_sl(lseg, tmp);
2859
2860 /*
2861 * ordinarily we should always find an intersection point, but that could
2862 * fail in the presence of NaN coordinates, and perhaps even from simple
2863 * roundoff issues. Return a SQL NULL if so.
2864 */
2865 if (result == NULL)
2866 PG_RETURN_NULL();
2867
2868 #ifdef GEODEBUG
2869 printf("close_ps- result.x %f result.y %f\n", result->x, result->y);
2870 #endif
2871 PG_RETURN_POINT_P(result);
2872 }
2873
2874
2875 /* close_lseg()
2876 * Closest point to l1 on l2.
2877 */
2878 Datum
close_lseg(PG_FUNCTION_ARGS)2879 close_lseg(PG_FUNCTION_ARGS)
2880 {
2881 LSEG *l1 = PG_GETARG_LSEG_P(0);
2882 LSEG *l2 = PG_GETARG_LSEG_P(1);
2883 Point *result = NULL;
2884 Point point;
2885 double dist;
2886 double d;
2887
2888 d = dist_ps_internal(&l1->p[0], l2);
2889 dist = d;
2890 memcpy(&point, &l1->p[0], sizeof(Point));
2891
2892 if ((d = dist_ps_internal(&l1->p[1], l2)) < dist)
2893 {
2894 dist = d;
2895 memcpy(&point, &l1->p[1], sizeof(Point));
2896 }
2897
2898 if (dist_ps_internal(&l2->p[0], l1) < dist)
2899 {
2900 result = DatumGetPointP(DirectFunctionCall2(close_ps,
2901 PointPGetDatum(&l2->p[0]),
2902 LsegPGetDatum(l1)));
2903 memcpy(&point, result, sizeof(Point));
2904 result = DatumGetPointP(DirectFunctionCall2(close_ps,
2905 PointPGetDatum(&point),
2906 LsegPGetDatum(l2)));
2907 }
2908
2909 if (dist_ps_internal(&l2->p[1], l1) < dist)
2910 {
2911 result = DatumGetPointP(DirectFunctionCall2(close_ps,
2912 PointPGetDatum(&l2->p[1]),
2913 LsegPGetDatum(l1)));
2914 memcpy(&point, result, sizeof(Point));
2915 result = DatumGetPointP(DirectFunctionCall2(close_ps,
2916 PointPGetDatum(&point),
2917 LsegPGetDatum(l2)));
2918 }
2919
2920 if (result == NULL)
2921 result = point_copy(&point);
2922
2923 PG_RETURN_POINT_P(result);
2924 }
2925
2926 /* close_pb()
2927 * Closest point on or in box to specified point.
2928 */
2929 Datum
close_pb(PG_FUNCTION_ARGS)2930 close_pb(PG_FUNCTION_ARGS)
2931 {
2932 Point *pt = PG_GETARG_POINT_P(0);
2933 BOX *box = PG_GETARG_BOX_P(1);
2934 LSEG lseg,
2935 seg;
2936 Point point;
2937 double dist,
2938 d;
2939
2940 if (DatumGetBool(DirectFunctionCall2(on_pb,
2941 PointPGetDatum(pt),
2942 BoxPGetDatum(box))))
2943 PG_RETURN_POINT_P(pt);
2944
2945 /* pairwise check lseg distances */
2946 point.x = box->low.x;
2947 point.y = box->high.y;
2948 statlseg_construct(&lseg, &box->low, &point);
2949 dist = dist_ps_internal(pt, &lseg);
2950
2951 statlseg_construct(&seg, &box->high, &point);
2952 if ((d = dist_ps_internal(pt, &seg)) < dist)
2953 {
2954 dist = d;
2955 memcpy(&lseg, &seg, sizeof(lseg));
2956 }
2957
2958 point.x = box->high.x;
2959 point.y = box->low.y;
2960 statlseg_construct(&seg, &box->low, &point);
2961 if ((d = dist_ps_internal(pt, &seg)) < dist)
2962 {
2963 dist = d;
2964 memcpy(&lseg, &seg, sizeof(lseg));
2965 }
2966
2967 statlseg_construct(&seg, &box->high, &point);
2968 if ((d = dist_ps_internal(pt, &seg)) < dist)
2969 {
2970 dist = d;
2971 memcpy(&lseg, &seg, sizeof(lseg));
2972 }
2973
2974 PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
2975 PointPGetDatum(pt),
2976 LsegPGetDatum(&lseg)));
2977 }
2978
2979 /* close_sl()
2980 * Closest point on line to line segment.
2981 *
2982 * XXX THIS CODE IS WRONG
2983 * The code is actually calculating the point on the line segment
2984 * which is backwards from the routine naming convention.
2985 * Copied code to new routine close_ls() but haven't fixed this one yet.
2986 * - thomas 1998-01-31
2987 */
2988 Datum
close_sl(PG_FUNCTION_ARGS)2989 close_sl(PG_FUNCTION_ARGS)
2990 {
2991 #ifdef NOT_USED
2992 LSEG *lseg = PG_GETARG_LSEG_P(0);
2993 LINE *line = PG_GETARG_LINE_P(1);
2994 Point *result;
2995 float8 d1,
2996 d2;
2997
2998 result = interpt_sl(lseg, line);
2999 if (result)
3000 PG_RETURN_POINT_P(result);
3001
3002 d1 = dist_pl_internal(&lseg->p[0], line);
3003 d2 = dist_pl_internal(&lseg->p[1], line);
3004 if (d1 < d2)
3005 result = point_copy(&lseg->p[0]);
3006 else
3007 result = point_copy(&lseg->p[1]);
3008
3009 PG_RETURN_POINT_P(result);
3010 #endif
3011
3012 ereport(ERROR,
3013 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3014 errmsg("function \"close_sl\" not implemented")));
3015
3016 PG_RETURN_NULL();
3017 }
3018
3019 /* close_ls()
3020 * Closest point on line segment to line.
3021 */
3022 Datum
close_ls(PG_FUNCTION_ARGS)3023 close_ls(PG_FUNCTION_ARGS)
3024 {
3025 LINE *line = PG_GETARG_LINE_P(0);
3026 LSEG *lseg = PG_GETARG_LSEG_P(1);
3027 Point *result;
3028 float8 d1,
3029 d2;
3030
3031 result = interpt_sl(lseg, line);
3032 if (result)
3033 PG_RETURN_POINT_P(result);
3034
3035 d1 = dist_pl_internal(&lseg->p[0], line);
3036 d2 = dist_pl_internal(&lseg->p[1], line);
3037 if (d1 < d2)
3038 result = point_copy(&lseg->p[0]);
3039 else
3040 result = point_copy(&lseg->p[1]);
3041
3042 PG_RETURN_POINT_P(result);
3043 }
3044
3045 /* close_sb()
3046 * Closest point on or in box to line segment.
3047 */
3048 Datum
close_sb(PG_FUNCTION_ARGS)3049 close_sb(PG_FUNCTION_ARGS)
3050 {
3051 LSEG *lseg = PG_GETARG_LSEG_P(0);
3052 BOX *box = PG_GETARG_BOX_P(1);
3053 Point point;
3054 LSEG bseg,
3055 seg;
3056 double dist,
3057 d;
3058
3059 /* segment intersects box? then just return closest point to center */
3060 if (DatumGetBool(DirectFunctionCall2(inter_sb,
3061 LsegPGetDatum(lseg),
3062 BoxPGetDatum(box))))
3063 {
3064 box_cn(&point, box);
3065 PG_RETURN_DATUM(DirectFunctionCall2(close_ps,
3066 PointPGetDatum(&point),
3067 LsegPGetDatum(lseg)));
3068 }
3069
3070 /* pairwise check lseg distances */
3071 point.x = box->low.x;
3072 point.y = box->high.y;
3073 statlseg_construct(&bseg, &box->low, &point);
3074 dist = lseg_dt(lseg, &bseg);
3075
3076 statlseg_construct(&seg, &box->high, &point);
3077 if ((d = lseg_dt(lseg, &seg)) < dist)
3078 {
3079 dist = d;
3080 memcpy(&bseg, &seg, sizeof(bseg));
3081 }
3082
3083 point.x = box->high.x;
3084 point.y = box->low.y;
3085 statlseg_construct(&seg, &box->low, &point);
3086 if ((d = lseg_dt(lseg, &seg)) < dist)
3087 {
3088 dist = d;
3089 memcpy(&bseg, &seg, sizeof(bseg));
3090 }
3091
3092 statlseg_construct(&seg, &box->high, &point);
3093 if ((d = lseg_dt(lseg, &seg)) < dist)
3094 {
3095 dist = d;
3096 memcpy(&bseg, &seg, sizeof(bseg));
3097 }
3098
3099 /* OK, we now have the closest line segment on the box boundary */
3100 PG_RETURN_DATUM(DirectFunctionCall2(close_lseg,
3101 LsegPGetDatum(lseg),
3102 LsegPGetDatum(&bseg)));
3103 }
3104
3105 Datum
close_lb(PG_FUNCTION_ARGS)3106 close_lb(PG_FUNCTION_ARGS)
3107 {
3108 #ifdef NOT_USED
3109 LINE *line = PG_GETARG_LINE_P(0);
3110 BOX *box = PG_GETARG_BOX_P(1);
3111 #endif
3112
3113 /* think about this one for a while */
3114 ereport(ERROR,
3115 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3116 errmsg("function \"close_lb\" not implemented")));
3117
3118 PG_RETURN_NULL();
3119 }
3120
3121 /*---------------------------------------------------------------------
3122 * on_
3123 * Whether one object lies completely within another.
3124 *-------------------------------------------------------------------*/
3125
3126 /* on_pl -
3127 * Does the point satisfy the equation?
3128 */
3129 Datum
on_pl(PG_FUNCTION_ARGS)3130 on_pl(PG_FUNCTION_ARGS)
3131 {
3132 Point *pt = PG_GETARG_POINT_P(0);
3133 LINE *line = PG_GETARG_LINE_P(1);
3134
3135 PG_RETURN_BOOL(FPzero(line->A * pt->x + line->B * pt->y + line->C));
3136 }
3137
3138
3139 /* on_ps -
3140 * Determine colinearity by detecting a triangle inequality.
3141 * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09
3142 */
3143 Datum
on_ps(PG_FUNCTION_ARGS)3144 on_ps(PG_FUNCTION_ARGS)
3145 {
3146 Point *pt = PG_GETARG_POINT_P(0);
3147 LSEG *lseg = PG_GETARG_LSEG_P(1);
3148
3149 PG_RETURN_BOOL(on_ps_internal(pt, lseg));
3150 }
3151
3152 static bool
on_ps_internal(Point * pt,LSEG * lseg)3153 on_ps_internal(Point *pt, LSEG *lseg)
3154 {
3155 return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]),
3156 point_dt(&lseg->p[0], &lseg->p[1]));
3157 }
3158
3159 Datum
on_pb(PG_FUNCTION_ARGS)3160 on_pb(PG_FUNCTION_ARGS)
3161 {
3162 Point *pt = PG_GETARG_POINT_P(0);
3163 BOX *box = PG_GETARG_BOX_P(1);
3164
3165 PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
3166 pt->y <= box->high.y && pt->y >= box->low.y);
3167 }
3168
3169 Datum
box_contain_pt(PG_FUNCTION_ARGS)3170 box_contain_pt(PG_FUNCTION_ARGS)
3171 {
3172 BOX *box = PG_GETARG_BOX_P(0);
3173 Point *pt = PG_GETARG_POINT_P(1);
3174
3175 PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x &&
3176 pt->y <= box->high.y && pt->y >= box->low.y);
3177 }
3178
3179 /* on_ppath -
3180 * Whether a point lies within (on) a polyline.
3181 * If open, we have to (groan) check each segment.
3182 * (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
3183 * If closed, we use the old O(n) ray method for point-in-polygon.
3184 * The ray is horizontal, from pt out to the right.
3185 * Each segment that crosses the ray counts as an
3186 * intersection; note that an endpoint or edge may touch
3187 * but not cross.
3188 * (we can do p-in-p in lg(n), but it takes preprocessing)
3189 */
3190 Datum
on_ppath(PG_FUNCTION_ARGS)3191 on_ppath(PG_FUNCTION_ARGS)
3192 {
3193 Point *pt = PG_GETARG_POINT_P(0);
3194 PATH *path = PG_GETARG_PATH_P(1);
3195 int i,
3196 n;
3197 double a,
3198 b;
3199
3200 /*-- OPEN --*/
3201 if (!path->closed)
3202 {
3203 n = path->npts - 1;
3204 a = point_dt(pt, &path->p[0]);
3205 for (i = 0; i < n; i++)
3206 {
3207 b = point_dt(pt, &path->p[i + 1]);
3208 if (FPeq(a + b,
3209 point_dt(&path->p[i], &path->p[i + 1])))
3210 PG_RETURN_BOOL(true);
3211 a = b;
3212 }
3213 PG_RETURN_BOOL(false);
3214 }
3215
3216 /*-- CLOSED --*/
3217 PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
3218 }
3219
3220 Datum
on_sl(PG_FUNCTION_ARGS)3221 on_sl(PG_FUNCTION_ARGS)
3222 {
3223 LSEG *lseg = PG_GETARG_LSEG_P(0);
3224 LINE *line = PG_GETARG_LINE_P(1);
3225
3226 PG_RETURN_BOOL(DatumGetBool(DirectFunctionCall2(on_pl,
3227 PointPGetDatum(&lseg->p[0]),
3228 LinePGetDatum(line))) &&
3229 DatumGetBool(DirectFunctionCall2(on_pl,
3230 PointPGetDatum(&lseg->p[1]),
3231 LinePGetDatum(line))));
3232 }
3233
3234 Datum
on_sb(PG_FUNCTION_ARGS)3235 on_sb(PG_FUNCTION_ARGS)
3236 {
3237 LSEG *lseg = PG_GETARG_LSEG_P(0);
3238 BOX *box = PG_GETARG_BOX_P(1);
3239
3240 PG_RETURN_BOOL(DatumGetBool(DirectFunctionCall2(on_pb,
3241 PointPGetDatum(&lseg->p[0]),
3242 BoxPGetDatum(box))) &&
3243 DatumGetBool(DirectFunctionCall2(on_pb,
3244 PointPGetDatum(&lseg->p[1]),
3245 BoxPGetDatum(box))));
3246 }
3247
3248 /*---------------------------------------------------------------------
3249 * inter_
3250 * Whether one object intersects another.
3251 *-------------------------------------------------------------------*/
3252
3253 Datum
inter_sl(PG_FUNCTION_ARGS)3254 inter_sl(PG_FUNCTION_ARGS)
3255 {
3256 LSEG *lseg = PG_GETARG_LSEG_P(0);
3257 LINE *line = PG_GETARG_LINE_P(1);
3258
3259 PG_RETURN_BOOL(has_interpt_sl(lseg, line));
3260 }
3261
3262 /* inter_sb()
3263 * Do line segment and box intersect?
3264 *
3265 * Segment completely inside box counts as intersection.
3266 * If you want only segments crossing box boundaries,
3267 * try converting box to path first.
3268 *
3269 * Optimize for non-intersection by checking for box intersection first.
3270 * - thomas 1998-01-30
3271 */
3272 Datum
inter_sb(PG_FUNCTION_ARGS)3273 inter_sb(PG_FUNCTION_ARGS)
3274 {
3275 LSEG *lseg = PG_GETARG_LSEG_P(0);
3276 BOX *box = PG_GETARG_BOX_P(1);
3277 BOX lbox;
3278 LSEG bseg;
3279 Point point;
3280
3281 lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x);
3282 lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y);
3283 lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x);
3284 lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y);
3285
3286 /* nothing close to overlap? then not going to intersect */
3287 if (!box_ov(&lbox, box))
3288 PG_RETURN_BOOL(false);
3289
3290 /* an endpoint of segment is inside box? then clearly intersects */
3291 if (DatumGetBool(DirectFunctionCall2(on_pb,
3292 PointPGetDatum(&lseg->p[0]),
3293 BoxPGetDatum(box))) ||
3294 DatumGetBool(DirectFunctionCall2(on_pb,
3295 PointPGetDatum(&lseg->p[1]),
3296 BoxPGetDatum(box))))
3297 PG_RETURN_BOOL(true);
3298
3299 /* pairwise check lseg intersections */
3300 point.x = box->low.x;
3301 point.y = box->high.y;
3302 statlseg_construct(&bseg, &box->low, &point);
3303 if (lseg_intersect_internal(&bseg, lseg))
3304 PG_RETURN_BOOL(true);
3305
3306 statlseg_construct(&bseg, &box->high, &point);
3307 if (lseg_intersect_internal(&bseg, lseg))
3308 PG_RETURN_BOOL(true);
3309
3310 point.x = box->high.x;
3311 point.y = box->low.y;
3312 statlseg_construct(&bseg, &box->low, &point);
3313 if (lseg_intersect_internal(&bseg, lseg))
3314 PG_RETURN_BOOL(true);
3315
3316 statlseg_construct(&bseg, &box->high, &point);
3317 if (lseg_intersect_internal(&bseg, lseg))
3318 PG_RETURN_BOOL(true);
3319
3320 /* if we dropped through, no two segs intersected */
3321 PG_RETURN_BOOL(false);
3322 }
3323
3324 /* inter_lb()
3325 * Do line and box intersect?
3326 */
3327 Datum
inter_lb(PG_FUNCTION_ARGS)3328 inter_lb(PG_FUNCTION_ARGS)
3329 {
3330 LINE *line = PG_GETARG_LINE_P(0);
3331 BOX *box = PG_GETARG_BOX_P(1);
3332 LSEG bseg;
3333 Point p1,
3334 p2;
3335
3336 /* pairwise check lseg intersections */
3337 p1.x = box->low.x;
3338 p1.y = box->low.y;
3339 p2.x = box->low.x;
3340 p2.y = box->high.y;
3341 statlseg_construct(&bseg, &p1, &p2);
3342 if (has_interpt_sl(&bseg, line))
3343 PG_RETURN_BOOL(true);
3344 p1.x = box->high.x;
3345 p1.y = box->high.y;
3346 statlseg_construct(&bseg, &p1, &p2);
3347 if (has_interpt_sl(&bseg, line))
3348 PG_RETURN_BOOL(true);
3349 p2.x = box->high.x;
3350 p2.y = box->low.y;
3351 statlseg_construct(&bseg, &p1, &p2);
3352 if (has_interpt_sl(&bseg, line))
3353 PG_RETURN_BOOL(true);
3354 p1.x = box->low.x;
3355 p1.y = box->low.y;
3356 statlseg_construct(&bseg, &p1, &p2);
3357 if (has_interpt_sl(&bseg, line))
3358 PG_RETURN_BOOL(true);
3359
3360 /* if we dropped through, no intersection */
3361 PG_RETURN_BOOL(false);
3362 }
3363
3364 /*------------------------------------------------------------------
3365 * The following routines define a data type and operator class for
3366 * POLYGONS .... Part of which (the polygon's bounding box) is built on
3367 * top of the BOX data type.
3368 *
3369 * make_bound_box - create the bounding box for the input polygon
3370 *------------------------------------------------------------------*/
3371
3372 /*---------------------------------------------------------------------
3373 * Make the smallest bounding box for the given polygon.
3374 *---------------------------------------------------------------------*/
3375 static void
make_bound_box(POLYGON * poly)3376 make_bound_box(POLYGON *poly)
3377 {
3378 int i;
3379 double x1,
3380 y1,
3381 x2,
3382 y2;
3383
3384 if (poly->npts > 0)
3385 {
3386 x2 = x1 = poly->p[0].x;
3387 y2 = y1 = poly->p[0].y;
3388 for (i = 1; i < poly->npts; i++)
3389 {
3390 if (poly->p[i].x < x1)
3391 x1 = poly->p[i].x;
3392 if (poly->p[i].x > x2)
3393 x2 = poly->p[i].x;
3394 if (poly->p[i].y < y1)
3395 y1 = poly->p[i].y;
3396 if (poly->p[i].y > y2)
3397 y2 = poly->p[i].y;
3398 }
3399
3400 box_fill(&(poly->boundbox), x1, x2, y1, y2);
3401 }
3402 else
3403 ereport(ERROR,
3404 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
3405 errmsg("cannot create bounding box for empty polygon")));
3406 }
3407
3408 /*------------------------------------------------------------------
3409 * poly_in - read in the polygon from a string specification
3410 *
3411 * External format:
3412 * "((x0,y0),...,(xn,yn))"
3413 * "x0,y0,...,xn,yn"
3414 * also supports the older style "(x1,...,xn,y1,...yn)"
3415 *------------------------------------------------------------------*/
3416 Datum
poly_in(PG_FUNCTION_ARGS)3417 poly_in(PG_FUNCTION_ARGS)
3418 {
3419 char *str = PG_GETARG_CSTRING(0);
3420 POLYGON *poly;
3421 int npts;
3422 int size;
3423 int base_size;
3424 bool isopen;
3425
3426 if ((npts = pair_count(str, ',')) <= 0)
3427 ereport(ERROR,
3428 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3429 errmsg("invalid input syntax for type %s: \"%s\"",
3430 "polygon", str)));
3431
3432 base_size = sizeof(poly->p[0]) * npts;
3433 size = offsetof(POLYGON, p) + base_size;
3434
3435 /* Check for integer overflow */
3436 if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
3437 ereport(ERROR,
3438 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
3439 errmsg("too many points requested")));
3440
3441 poly = (POLYGON *) palloc0(size); /* zero any holes */
3442
3443 SET_VARSIZE(poly, size);
3444 poly->npts = npts;
3445
3446 path_decode(str, false, npts, &(poly->p[0]), &isopen, NULL, "polygon", str);
3447
3448 make_bound_box(poly);
3449
3450 PG_RETURN_POLYGON_P(poly);
3451 }
3452
3453 /*---------------------------------------------------------------
3454 * poly_out - convert internal POLYGON representation to the
3455 * character string format "((f8,f8),...,(f8,f8))"
3456 *---------------------------------------------------------------*/
3457 Datum
poly_out(PG_FUNCTION_ARGS)3458 poly_out(PG_FUNCTION_ARGS)
3459 {
3460 POLYGON *poly = PG_GETARG_POLYGON_P(0);
3461
3462 PG_RETURN_CSTRING(path_encode(PATH_CLOSED, poly->npts, poly->p));
3463 }
3464
3465 /*
3466 * poly_recv - converts external binary format to polygon
3467 *
3468 * External representation is int32 number of points, and the points.
3469 * We recompute the bounding box on read, instead of trusting it to
3470 * be valid. (Checking it would take just as long, so may as well
3471 * omit it from external representation.)
3472 */
3473 Datum
poly_recv(PG_FUNCTION_ARGS)3474 poly_recv(PG_FUNCTION_ARGS)
3475 {
3476 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
3477 POLYGON *poly;
3478 int32 npts;
3479 int32 i;
3480 int size;
3481
3482 npts = pq_getmsgint(buf, sizeof(int32));
3483 if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p)) / sizeof(Point)))
3484 ereport(ERROR,
3485 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
3486 errmsg("invalid number of points in external \"polygon\" value")));
3487
3488 size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * npts;
3489 poly = (POLYGON *) palloc0(size); /* zero any holes */
3490
3491 SET_VARSIZE(poly, size);
3492 poly->npts = npts;
3493
3494 for (i = 0; i < npts; i++)
3495 {
3496 poly->p[i].x = pq_getmsgfloat8(buf);
3497 poly->p[i].y = pq_getmsgfloat8(buf);
3498 }
3499
3500 make_bound_box(poly);
3501
3502 PG_RETURN_POLYGON_P(poly);
3503 }
3504
3505 /*
3506 * poly_send - converts polygon to binary format
3507 */
3508 Datum
poly_send(PG_FUNCTION_ARGS)3509 poly_send(PG_FUNCTION_ARGS)
3510 {
3511 POLYGON *poly = PG_GETARG_POLYGON_P(0);
3512 StringInfoData buf;
3513 int32 i;
3514
3515 pq_begintypsend(&buf);
3516 pq_sendint32(&buf, poly->npts);
3517 for (i = 0; i < poly->npts; i++)
3518 {
3519 pq_sendfloat8(&buf, poly->p[i].x);
3520 pq_sendfloat8(&buf, poly->p[i].y);
3521 }
3522 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
3523 }
3524
3525
3526 /*-------------------------------------------------------
3527 * Is polygon A strictly left of polygon B? i.e. is
3528 * the right most point of A left of the left most point
3529 * of B?
3530 *-------------------------------------------------------*/
3531 Datum
poly_left(PG_FUNCTION_ARGS)3532 poly_left(PG_FUNCTION_ARGS)
3533 {
3534 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3535 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3536 bool result;
3537
3538 result = polya->boundbox.high.x < polyb->boundbox.low.x;
3539
3540 /*
3541 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3542 */
3543 PG_FREE_IF_COPY(polya, 0);
3544 PG_FREE_IF_COPY(polyb, 1);
3545
3546 PG_RETURN_BOOL(result);
3547 }
3548
3549 /*-------------------------------------------------------
3550 * Is polygon A overlapping or left of polygon B? i.e. is
3551 * the right most point of A at or left of the right most point
3552 * of B?
3553 *-------------------------------------------------------*/
3554 Datum
poly_overleft(PG_FUNCTION_ARGS)3555 poly_overleft(PG_FUNCTION_ARGS)
3556 {
3557 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3558 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3559 bool result;
3560
3561 result = polya->boundbox.high.x <= polyb->boundbox.high.x;
3562
3563 /*
3564 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3565 */
3566 PG_FREE_IF_COPY(polya, 0);
3567 PG_FREE_IF_COPY(polyb, 1);
3568
3569 PG_RETURN_BOOL(result);
3570 }
3571
3572 /*-------------------------------------------------------
3573 * Is polygon A strictly right of polygon B? i.e. is
3574 * the left most point of A right of the right most point
3575 * of B?
3576 *-------------------------------------------------------*/
3577 Datum
poly_right(PG_FUNCTION_ARGS)3578 poly_right(PG_FUNCTION_ARGS)
3579 {
3580 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3581 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3582 bool result;
3583
3584 result = polya->boundbox.low.x > polyb->boundbox.high.x;
3585
3586 /*
3587 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3588 */
3589 PG_FREE_IF_COPY(polya, 0);
3590 PG_FREE_IF_COPY(polyb, 1);
3591
3592 PG_RETURN_BOOL(result);
3593 }
3594
3595 /*-------------------------------------------------------
3596 * Is polygon A overlapping or right of polygon B? i.e. is
3597 * the left most point of A at or right of the left most point
3598 * of B?
3599 *-------------------------------------------------------*/
3600 Datum
poly_overright(PG_FUNCTION_ARGS)3601 poly_overright(PG_FUNCTION_ARGS)
3602 {
3603 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3604 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3605 bool result;
3606
3607 result = polya->boundbox.low.x >= polyb->boundbox.low.x;
3608
3609 /*
3610 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3611 */
3612 PG_FREE_IF_COPY(polya, 0);
3613 PG_FREE_IF_COPY(polyb, 1);
3614
3615 PG_RETURN_BOOL(result);
3616 }
3617
3618 /*-------------------------------------------------------
3619 * Is polygon A strictly below polygon B? i.e. is
3620 * the upper most point of A below the lower most point
3621 * of B?
3622 *-------------------------------------------------------*/
3623 Datum
poly_below(PG_FUNCTION_ARGS)3624 poly_below(PG_FUNCTION_ARGS)
3625 {
3626 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3627 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3628 bool result;
3629
3630 result = polya->boundbox.high.y < polyb->boundbox.low.y;
3631
3632 /*
3633 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3634 */
3635 PG_FREE_IF_COPY(polya, 0);
3636 PG_FREE_IF_COPY(polyb, 1);
3637
3638 PG_RETURN_BOOL(result);
3639 }
3640
3641 /*-------------------------------------------------------
3642 * Is polygon A overlapping or below polygon B? i.e. is
3643 * the upper most point of A at or below the upper most point
3644 * of B?
3645 *-------------------------------------------------------*/
3646 Datum
poly_overbelow(PG_FUNCTION_ARGS)3647 poly_overbelow(PG_FUNCTION_ARGS)
3648 {
3649 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3650 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3651 bool result;
3652
3653 result = polya->boundbox.high.y <= polyb->boundbox.high.y;
3654
3655 /*
3656 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3657 */
3658 PG_FREE_IF_COPY(polya, 0);
3659 PG_FREE_IF_COPY(polyb, 1);
3660
3661 PG_RETURN_BOOL(result);
3662 }
3663
3664 /*-------------------------------------------------------
3665 * Is polygon A strictly above polygon B? i.e. is
3666 * the lower most point of A above the upper most point
3667 * of B?
3668 *-------------------------------------------------------*/
3669 Datum
poly_above(PG_FUNCTION_ARGS)3670 poly_above(PG_FUNCTION_ARGS)
3671 {
3672 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3673 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3674 bool result;
3675
3676 result = polya->boundbox.low.y > polyb->boundbox.high.y;
3677
3678 /*
3679 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3680 */
3681 PG_FREE_IF_COPY(polya, 0);
3682 PG_FREE_IF_COPY(polyb, 1);
3683
3684 PG_RETURN_BOOL(result);
3685 }
3686
3687 /*-------------------------------------------------------
3688 * Is polygon A overlapping or above polygon B? i.e. is
3689 * the lower most point of A at or above the lower most point
3690 * of B?
3691 *-------------------------------------------------------*/
3692 Datum
poly_overabove(PG_FUNCTION_ARGS)3693 poly_overabove(PG_FUNCTION_ARGS)
3694 {
3695 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3696 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3697 bool result;
3698
3699 result = polya->boundbox.low.y >= polyb->boundbox.low.y;
3700
3701 /*
3702 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3703 */
3704 PG_FREE_IF_COPY(polya, 0);
3705 PG_FREE_IF_COPY(polyb, 1);
3706
3707 PG_RETURN_BOOL(result);
3708 }
3709
3710
3711 /*-------------------------------------------------------
3712 * Is polygon A the same as polygon B? i.e. are all the
3713 * points the same?
3714 * Check all points for matches in both forward and reverse
3715 * direction since polygons are non-directional and are
3716 * closed shapes.
3717 *-------------------------------------------------------*/
3718 Datum
poly_same(PG_FUNCTION_ARGS)3719 poly_same(PG_FUNCTION_ARGS)
3720 {
3721 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3722 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3723 bool result;
3724
3725 if (polya->npts != polyb->npts)
3726 result = false;
3727 else
3728 result = plist_same(polya->npts, polya->p, polyb->p);
3729
3730 /*
3731 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3732 */
3733 PG_FREE_IF_COPY(polya, 0);
3734 PG_FREE_IF_COPY(polyb, 1);
3735
3736 PG_RETURN_BOOL(result);
3737 }
3738
3739 /*-----------------------------------------------------------------
3740 * Determine if polygon A overlaps polygon B
3741 *-----------------------------------------------------------------*/
3742 Datum
poly_overlap(PG_FUNCTION_ARGS)3743 poly_overlap(PG_FUNCTION_ARGS)
3744 {
3745 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3746 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3747 bool result;
3748
3749 /* Quick check by bounding box */
3750 result = (polya->npts > 0 && polyb->npts > 0 &&
3751 box_ov(&polya->boundbox, &polyb->boundbox)) ? true : false;
3752
3753 /*
3754 * Brute-force algorithm - try to find intersected edges, if so then
3755 * polygons are overlapped else check is one polygon inside other or not
3756 * by testing single point of them.
3757 */
3758 if (result)
3759 {
3760 int ia,
3761 ib;
3762 LSEG sa,
3763 sb;
3764
3765 /* Init first of polya's edge with last point */
3766 sa.p[0] = polya->p[polya->npts - 1];
3767 result = false;
3768
3769 for (ia = 0; ia < polya->npts && result == false; ia++)
3770 {
3771 /* Second point of polya's edge is a current one */
3772 sa.p[1] = polya->p[ia];
3773
3774 /* Init first of polyb's edge with last point */
3775 sb.p[0] = polyb->p[polyb->npts - 1];
3776
3777 for (ib = 0; ib < polyb->npts && result == false; ib++)
3778 {
3779 sb.p[1] = polyb->p[ib];
3780 result = lseg_intersect_internal(&sa, &sb);
3781 sb.p[0] = sb.p[1];
3782 }
3783
3784 /*
3785 * move current endpoint to the first point of next edge
3786 */
3787 sa.p[0] = sa.p[1];
3788 }
3789
3790 if (result == false)
3791 {
3792 result = (point_inside(polya->p, polyb->npts, polyb->p)
3793 ||
3794 point_inside(polyb->p, polya->npts, polya->p));
3795 }
3796 }
3797
3798 /*
3799 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3800 */
3801 PG_FREE_IF_COPY(polya, 0);
3802 PG_FREE_IF_COPY(polyb, 1);
3803
3804 PG_RETURN_BOOL(result);
3805 }
3806
3807 /*
3808 * Tests special kind of segment for in/out of polygon.
3809 * Special kind means:
3810 * - point a should be on segment s
3811 * - segment (a,b) should not be contained by s
3812 * Returns true if:
3813 * - segment (a,b) is collinear to s and (a,b) is in polygon
3814 * - segment (a,b) s not collinear to s. Note: that doesn't
3815 * mean that segment is in polygon!
3816 */
3817
3818 static bool
touched_lseg_inside_poly(Point * a,Point * b,LSEG * s,POLYGON * poly,int start)3819 touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
3820 {
3821 /* point a is on s, b is not */
3822 LSEG t;
3823
3824 t.p[0] = *a;
3825 t.p[1] = *b;
3826
3827 #define POINTEQ(pt1, pt2) (FPeq((pt1)->x, (pt2)->x) && FPeq((pt1)->y, (pt2)->y))
3828 if (POINTEQ(a, s->p))
3829 {
3830 if (on_ps_internal(s->p + 1, &t))
3831 return lseg_inside_poly(b, s->p + 1, poly, start);
3832 }
3833 else if (POINTEQ(a, s->p + 1))
3834 {
3835 if (on_ps_internal(s->p, &t))
3836 return lseg_inside_poly(b, s->p, poly, start);
3837 }
3838 else if (on_ps_internal(s->p, &t))
3839 {
3840 return lseg_inside_poly(b, s->p, poly, start);
3841 }
3842 else if (on_ps_internal(s->p + 1, &t))
3843 {
3844 return lseg_inside_poly(b, s->p + 1, poly, start);
3845 }
3846
3847 return true; /* may be not true, but that will check later */
3848 }
3849
3850 /*
3851 * Returns true if segment (a,b) is in polygon, option
3852 * start is used for optimization - function checks
3853 * polygon's edges started from start
3854 */
3855 static bool
lseg_inside_poly(Point * a,Point * b,POLYGON * poly,int start)3856 lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
3857 {
3858 LSEG s,
3859 t;
3860 int i;
3861 bool res = true,
3862 intersection = false;
3863
3864 t.p[0] = *a;
3865 t.p[1] = *b;
3866 s.p[0] = poly->p[(start == 0) ? (poly->npts - 1) : (start - 1)];
3867
3868 for (i = start; i < poly->npts && res; i++)
3869 {
3870 Point *interpt;
3871
3872 CHECK_FOR_INTERRUPTS();
3873
3874 s.p[1] = poly->p[i];
3875
3876 if (on_ps_internal(t.p, &s))
3877 {
3878 if (on_ps_internal(t.p + 1, &s))
3879 return true; /* t is contained by s */
3880
3881 /* Y-cross */
3882 res = touched_lseg_inside_poly(t.p, t.p + 1, &s, poly, i + 1);
3883 }
3884 else if (on_ps_internal(t.p + 1, &s))
3885 {
3886 /* Y-cross */
3887 res = touched_lseg_inside_poly(t.p + 1, t.p, &s, poly, i + 1);
3888 }
3889 else if ((interpt = lseg_interpt_internal(&t, &s)) != NULL)
3890 {
3891 /*
3892 * segments are X-crossing, go to check each subsegment
3893 */
3894
3895 intersection = true;
3896 res = lseg_inside_poly(t.p, interpt, poly, i + 1);
3897 if (res)
3898 res = lseg_inside_poly(t.p + 1, interpt, poly, i + 1);
3899 pfree(interpt);
3900 }
3901
3902 s.p[0] = s.p[1];
3903 }
3904
3905 if (res && !intersection)
3906 {
3907 Point p;
3908
3909 /*
3910 * if X-intersection wasn't found then check central point of tested
3911 * segment. In opposite case we already check all subsegments
3912 */
3913 p.x = (t.p[0].x + t.p[1].x) / 2.0;
3914 p.y = (t.p[0].y + t.p[1].y) / 2.0;
3915
3916 res = point_inside(&p, poly->npts, poly->p);
3917 }
3918
3919 return res;
3920 }
3921
3922 /*-----------------------------------------------------------------
3923 * Determine if polygon A contains polygon B.
3924 *-----------------------------------------------------------------*/
3925 Datum
poly_contain(PG_FUNCTION_ARGS)3926 poly_contain(PG_FUNCTION_ARGS)
3927 {
3928 POLYGON *polya = PG_GETARG_POLYGON_P(0);
3929 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3930 bool result;
3931
3932 /*
3933 * Quick check to see if bounding box is contained.
3934 */
3935 if (polya->npts > 0 && polyb->npts > 0 &&
3936 DatumGetBool(DirectFunctionCall2(box_contain,
3937 BoxPGetDatum(&polya->boundbox),
3938 BoxPGetDatum(&polyb->boundbox))))
3939 {
3940 int i;
3941 LSEG s;
3942
3943 s.p[0] = polyb->p[polyb->npts - 1];
3944 result = true;
3945
3946 for (i = 0; i < polyb->npts && result; i++)
3947 {
3948 s.p[1] = polyb->p[i];
3949 result = lseg_inside_poly(s.p, s.p + 1, polya, 0);
3950 s.p[0] = s.p[1];
3951 }
3952 }
3953 else
3954 {
3955 result = false;
3956 }
3957
3958 /*
3959 * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3960 */
3961 PG_FREE_IF_COPY(polya, 0);
3962 PG_FREE_IF_COPY(polyb, 1);
3963
3964 PG_RETURN_BOOL(result);
3965 }
3966
3967
3968 /*-----------------------------------------------------------------
3969 * Determine if polygon A is contained by polygon B
3970 *-----------------------------------------------------------------*/
3971 Datum
poly_contained(PG_FUNCTION_ARGS)3972 poly_contained(PG_FUNCTION_ARGS)
3973 {
3974 Datum polya = PG_GETARG_DATUM(0);
3975 Datum polyb = PG_GETARG_DATUM(1);
3976
3977 /* Just switch the arguments and pass it off to poly_contain */
3978 PG_RETURN_DATUM(DirectFunctionCall2(poly_contain, polyb, polya));
3979 }
3980
3981
3982 Datum
poly_contain_pt(PG_FUNCTION_ARGS)3983 poly_contain_pt(PG_FUNCTION_ARGS)
3984 {
3985 POLYGON *poly = PG_GETARG_POLYGON_P(0);
3986 Point *p = PG_GETARG_POINT_P(1);
3987
3988 PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
3989 }
3990
3991 Datum
pt_contained_poly(PG_FUNCTION_ARGS)3992 pt_contained_poly(PG_FUNCTION_ARGS)
3993 {
3994 Point *p = PG_GETARG_POINT_P(0);
3995 POLYGON *poly = PG_GETARG_POLYGON_P(1);
3996
3997 PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
3998 }
3999
4000
4001 Datum
poly_distance(PG_FUNCTION_ARGS)4002 poly_distance(PG_FUNCTION_ARGS)
4003 {
4004 #ifdef NOT_USED
4005 POLYGON *polya = PG_GETARG_POLYGON_P(0);
4006 POLYGON *polyb = PG_GETARG_POLYGON_P(1);
4007 #endif
4008
4009 ereport(ERROR,
4010 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4011 errmsg("function \"poly_distance\" not implemented")));
4012
4013 PG_RETURN_NULL();
4014 }
4015
4016
4017 /***********************************************************************
4018 **
4019 ** Routines for 2D points.
4020 **
4021 ***********************************************************************/
4022
4023 Datum
construct_point(PG_FUNCTION_ARGS)4024 construct_point(PG_FUNCTION_ARGS)
4025 {
4026 float8 x = PG_GETARG_FLOAT8(0);
4027 float8 y = PG_GETARG_FLOAT8(1);
4028
4029 PG_RETURN_POINT_P(point_construct(x, y));
4030 }
4031
4032 Datum
point_add(PG_FUNCTION_ARGS)4033 point_add(PG_FUNCTION_ARGS)
4034 {
4035 Point *p1 = PG_GETARG_POINT_P(0);
4036 Point *p2 = PG_GETARG_POINT_P(1);
4037 Point *result;
4038
4039 result = (Point *) palloc(sizeof(Point));
4040
4041 result->x = (p1->x + p2->x);
4042 result->y = (p1->y + p2->y);
4043
4044 PG_RETURN_POINT_P(result);
4045 }
4046
4047 Datum
point_sub(PG_FUNCTION_ARGS)4048 point_sub(PG_FUNCTION_ARGS)
4049 {
4050 Point *p1 = PG_GETARG_POINT_P(0);
4051 Point *p2 = PG_GETARG_POINT_P(1);
4052 Point *result;
4053
4054 result = (Point *) palloc(sizeof(Point));
4055
4056 result->x = (p1->x - p2->x);
4057 result->y = (p1->y - p2->y);
4058
4059 PG_RETURN_POINT_P(result);
4060 }
4061
4062 Datum
point_mul(PG_FUNCTION_ARGS)4063 point_mul(PG_FUNCTION_ARGS)
4064 {
4065 Point *p1 = PG_GETARG_POINT_P(0);
4066 Point *p2 = PG_GETARG_POINT_P(1);
4067 Point *result;
4068
4069 result = (Point *) palloc(sizeof(Point));
4070
4071 result->x = (p1->x * p2->x) - (p1->y * p2->y);
4072 result->y = (p1->x * p2->y) + (p1->y * p2->x);
4073
4074 PG_RETURN_POINT_P(result);
4075 }
4076
4077 Datum
point_div(PG_FUNCTION_ARGS)4078 point_div(PG_FUNCTION_ARGS)
4079 {
4080 Point *p1 = PG_GETARG_POINT_P(0);
4081 Point *p2 = PG_GETARG_POINT_P(1);
4082 Point *result;
4083 double div;
4084
4085 result = (Point *) palloc(sizeof(Point));
4086
4087 div = (p2->x * p2->x) + (p2->y * p2->y);
4088
4089 if (div == 0.0)
4090 ereport(ERROR,
4091 (errcode(ERRCODE_DIVISION_BY_ZERO),
4092 errmsg("division by zero")));
4093
4094 result->x = ((p1->x * p2->x) + (p1->y * p2->y)) / div;
4095 result->y = ((p2->x * p1->y) - (p2->y * p1->x)) / div;
4096
4097 PG_RETURN_POINT_P(result);
4098 }
4099
4100
4101 /***********************************************************************
4102 **
4103 ** Routines for 2D boxes.
4104 **
4105 ***********************************************************************/
4106
4107 Datum
points_box(PG_FUNCTION_ARGS)4108 points_box(PG_FUNCTION_ARGS)
4109 {
4110 Point *p1 = PG_GETARG_POINT_P(0);
4111 Point *p2 = PG_GETARG_POINT_P(1);
4112
4113 PG_RETURN_BOX_P(box_construct(p1->x, p2->x, p1->y, p2->y));
4114 }
4115
4116 Datum
box_add(PG_FUNCTION_ARGS)4117 box_add(PG_FUNCTION_ARGS)
4118 {
4119 BOX *box = PG_GETARG_BOX_P(0);
4120 Point *p = PG_GETARG_POINT_P(1);
4121
4122 PG_RETURN_BOX_P(box_construct((box->high.x + p->x),
4123 (box->low.x + p->x),
4124 (box->high.y + p->y),
4125 (box->low.y + p->y)));
4126 }
4127
4128 Datum
box_sub(PG_FUNCTION_ARGS)4129 box_sub(PG_FUNCTION_ARGS)
4130 {
4131 BOX *box = PG_GETARG_BOX_P(0);
4132 Point *p = PG_GETARG_POINT_P(1);
4133
4134 PG_RETURN_BOX_P(box_construct((box->high.x - p->x),
4135 (box->low.x - p->x),
4136 (box->high.y - p->y),
4137 (box->low.y - p->y)));
4138 }
4139
4140 Datum
box_mul(PG_FUNCTION_ARGS)4141 box_mul(PG_FUNCTION_ARGS)
4142 {
4143 BOX *box = PG_GETARG_BOX_P(0);
4144 Point *p = PG_GETARG_POINT_P(1);
4145 BOX *result;
4146 Point *high,
4147 *low;
4148
4149 high = DatumGetPointP(DirectFunctionCall2(point_mul,
4150 PointPGetDatum(&box->high),
4151 PointPGetDatum(p)));
4152 low = DatumGetPointP(DirectFunctionCall2(point_mul,
4153 PointPGetDatum(&box->low),
4154 PointPGetDatum(p)));
4155
4156 result = box_construct(high->x, low->x, high->y, low->y);
4157
4158 PG_RETURN_BOX_P(result);
4159 }
4160
4161 Datum
box_div(PG_FUNCTION_ARGS)4162 box_div(PG_FUNCTION_ARGS)
4163 {
4164 BOX *box = PG_GETARG_BOX_P(0);
4165 Point *p = PG_GETARG_POINT_P(1);
4166 BOX *result;
4167 Point *high,
4168 *low;
4169
4170 high = DatumGetPointP(DirectFunctionCall2(point_div,
4171 PointPGetDatum(&box->high),
4172 PointPGetDatum(p)));
4173 low = DatumGetPointP(DirectFunctionCall2(point_div,
4174 PointPGetDatum(&box->low),
4175 PointPGetDatum(p)));
4176
4177 result = box_construct(high->x, low->x, high->y, low->y);
4178
4179 PG_RETURN_BOX_P(result);
4180 }
4181
4182 /*
4183 * Convert point to empty box
4184 */
4185 Datum
point_box(PG_FUNCTION_ARGS)4186 point_box(PG_FUNCTION_ARGS)
4187 {
4188 Point *pt = PG_GETARG_POINT_P(0);
4189 BOX *box;
4190
4191 box = (BOX *) palloc(sizeof(BOX));
4192
4193 box->high.x = pt->x;
4194 box->low.x = pt->x;
4195 box->high.y = pt->y;
4196 box->low.y = pt->y;
4197
4198 PG_RETURN_BOX_P(box);
4199 }
4200
4201 /*
4202 * Smallest bounding box that includes both of the given boxes
4203 */
4204 Datum
boxes_bound_box(PG_FUNCTION_ARGS)4205 boxes_bound_box(PG_FUNCTION_ARGS)
4206 {
4207 BOX *box1 = PG_GETARG_BOX_P(0),
4208 *box2 = PG_GETARG_BOX_P(1),
4209 *container;
4210
4211 container = (BOX *) palloc(sizeof(BOX));
4212
4213 container->high.x = Max(box1->high.x, box2->high.x);
4214 container->low.x = Min(box1->low.x, box2->low.x);
4215 container->high.y = Max(box1->high.y, box2->high.y);
4216 container->low.y = Min(box1->low.y, box2->low.y);
4217
4218 PG_RETURN_BOX_P(container);
4219 }
4220
4221
4222 /***********************************************************************
4223 **
4224 ** Routines for 2D paths.
4225 **
4226 ***********************************************************************/
4227
4228 /* path_add()
4229 * Concatenate two paths (only if they are both open).
4230 */
4231 Datum
path_add(PG_FUNCTION_ARGS)4232 path_add(PG_FUNCTION_ARGS)
4233 {
4234 PATH *p1 = PG_GETARG_PATH_P(0);
4235 PATH *p2 = PG_GETARG_PATH_P(1);
4236 PATH *result;
4237 int size,
4238 base_size;
4239 int i;
4240
4241 if (p1->closed || p2->closed)
4242 PG_RETURN_NULL();
4243
4244 base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts);
4245 size = offsetof(PATH, p) + base_size;
4246
4247 /* Check for integer overflow */
4248 if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) ||
4249 size <= base_size)
4250 ereport(ERROR,
4251 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
4252 errmsg("too many points requested")));
4253
4254 result = (PATH *) palloc(size);
4255
4256 SET_VARSIZE(result, size);
4257 result->npts = (p1->npts + p2->npts);
4258 result->closed = p1->closed;
4259 /* prevent instability in unused pad bytes */
4260 result->dummy = 0;
4261
4262 for (i = 0; i < p1->npts; i++)
4263 {
4264 result->p[i].x = p1->p[i].x;
4265 result->p[i].y = p1->p[i].y;
4266 }
4267 for (i = 0; i < p2->npts; i++)
4268 {
4269 result->p[i + p1->npts].x = p2->p[i].x;
4270 result->p[i + p1->npts].y = p2->p[i].y;
4271 }
4272
4273 PG_RETURN_PATH_P(result);
4274 }
4275
4276 /* path_add_pt()
4277 * Translation operators.
4278 */
4279 Datum
path_add_pt(PG_FUNCTION_ARGS)4280 path_add_pt(PG_FUNCTION_ARGS)
4281 {
4282 PATH *path = PG_GETARG_PATH_P_COPY(0);
4283 Point *point = PG_GETARG_POINT_P(1);
4284 int i;
4285
4286 for (i = 0; i < path->npts; i++)
4287 {
4288 path->p[i].x += point->x;
4289 path->p[i].y += point->y;
4290 }
4291
4292 PG_RETURN_PATH_P(path);
4293 }
4294
4295 Datum
path_sub_pt(PG_FUNCTION_ARGS)4296 path_sub_pt(PG_FUNCTION_ARGS)
4297 {
4298 PATH *path = PG_GETARG_PATH_P_COPY(0);
4299 Point *point = PG_GETARG_POINT_P(1);
4300 int i;
4301
4302 for (i = 0; i < path->npts; i++)
4303 {
4304 path->p[i].x -= point->x;
4305 path->p[i].y -= point->y;
4306 }
4307
4308 PG_RETURN_PATH_P(path);
4309 }
4310
4311 /* path_mul_pt()
4312 * Rotation and scaling operators.
4313 */
4314 Datum
path_mul_pt(PG_FUNCTION_ARGS)4315 path_mul_pt(PG_FUNCTION_ARGS)
4316 {
4317 PATH *path = PG_GETARG_PATH_P_COPY(0);
4318 Point *point = PG_GETARG_POINT_P(1);
4319 Point *p;
4320 int i;
4321
4322 for (i = 0; i < path->npts; i++)
4323 {
4324 p = DatumGetPointP(DirectFunctionCall2(point_mul,
4325 PointPGetDatum(&path->p[i]),
4326 PointPGetDatum(point)));
4327 path->p[i].x = p->x;
4328 path->p[i].y = p->y;
4329 }
4330
4331 PG_RETURN_PATH_P(path);
4332 }
4333
4334 Datum
path_div_pt(PG_FUNCTION_ARGS)4335 path_div_pt(PG_FUNCTION_ARGS)
4336 {
4337 PATH *path = PG_GETARG_PATH_P_COPY(0);
4338 Point *point = PG_GETARG_POINT_P(1);
4339 Point *p;
4340 int i;
4341
4342 for (i = 0; i < path->npts; i++)
4343 {
4344 p = DatumGetPointP(DirectFunctionCall2(point_div,
4345 PointPGetDatum(&path->p[i]),
4346 PointPGetDatum(point)));
4347 path->p[i].x = p->x;
4348 path->p[i].y = p->y;
4349 }
4350
4351 PG_RETURN_PATH_P(path);
4352 }
4353
4354
4355 Datum
path_center(PG_FUNCTION_ARGS)4356 path_center(PG_FUNCTION_ARGS)
4357 {
4358 #ifdef NOT_USED
4359 PATH *path = PG_GETARG_PATH_P(0);
4360 #endif
4361
4362 ereport(ERROR,
4363 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4364 errmsg("function \"path_center\" not implemented")));
4365
4366 PG_RETURN_NULL();
4367 }
4368
4369 Datum
path_poly(PG_FUNCTION_ARGS)4370 path_poly(PG_FUNCTION_ARGS)
4371 {
4372 PATH *path = PG_GETARG_PATH_P(0);
4373 POLYGON *poly;
4374 int size;
4375 int i;
4376
4377 /* This is not very consistent --- other similar cases return NULL ... */
4378 if (!path->closed)
4379 ereport(ERROR,
4380 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
4381 errmsg("open path cannot be converted to polygon")));
4382
4383 /*
4384 * Never overflows: the old size fit in MaxAllocSize, and the new size is
4385 * just a small constant larger.
4386 */
4387 size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * path->npts;
4388 poly = (POLYGON *) palloc(size);
4389
4390 SET_VARSIZE(poly, size);
4391 poly->npts = path->npts;
4392
4393 for (i = 0; i < path->npts; i++)
4394 {
4395 poly->p[i].x = path->p[i].x;
4396 poly->p[i].y = path->p[i].y;
4397 }
4398
4399 make_bound_box(poly);
4400
4401 PG_RETURN_POLYGON_P(poly);
4402 }
4403
4404
4405 /***********************************************************************
4406 **
4407 ** Routines for 2D polygons.
4408 **
4409 ***********************************************************************/
4410
4411 Datum
poly_npoints(PG_FUNCTION_ARGS)4412 poly_npoints(PG_FUNCTION_ARGS)
4413 {
4414 POLYGON *poly = PG_GETARG_POLYGON_P(0);
4415
4416 PG_RETURN_INT32(poly->npts);
4417 }
4418
4419
4420 Datum
poly_center(PG_FUNCTION_ARGS)4421 poly_center(PG_FUNCTION_ARGS)
4422 {
4423 POLYGON *poly = PG_GETARG_POLYGON_P(0);
4424 Datum result;
4425 CIRCLE *circle;
4426
4427 circle = DatumGetCircleP(DirectFunctionCall1(poly_circle,
4428 PolygonPGetDatum(poly)));
4429 result = DirectFunctionCall1(circle_center,
4430 CirclePGetDatum(circle));
4431
4432 PG_RETURN_DATUM(result);
4433 }
4434
4435
4436 Datum
poly_box(PG_FUNCTION_ARGS)4437 poly_box(PG_FUNCTION_ARGS)
4438 {
4439 POLYGON *poly = PG_GETARG_POLYGON_P(0);
4440 BOX *box;
4441
4442 if (poly->npts < 1)
4443 PG_RETURN_NULL();
4444
4445 box = box_copy(&poly->boundbox);
4446
4447 PG_RETURN_BOX_P(box);
4448 }
4449
4450
4451 /* box_poly()
4452 * Convert a box to a polygon.
4453 */
4454 Datum
box_poly(PG_FUNCTION_ARGS)4455 box_poly(PG_FUNCTION_ARGS)
4456 {
4457 BOX *box = PG_GETARG_BOX_P(0);
4458 POLYGON *poly;
4459 int size;
4460
4461 /* map four corners of the box to a polygon */
4462 size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * 4;
4463 poly = (POLYGON *) palloc(size);
4464
4465 SET_VARSIZE(poly, size);
4466 poly->npts = 4;
4467
4468 poly->p[0].x = box->low.x;
4469 poly->p[0].y = box->low.y;
4470 poly->p[1].x = box->low.x;
4471 poly->p[1].y = box->high.y;
4472 poly->p[2].x = box->high.x;
4473 poly->p[2].y = box->high.y;
4474 poly->p[3].x = box->high.x;
4475 poly->p[3].y = box->low.y;
4476
4477 box_fill(&poly->boundbox, box->high.x, box->low.x,
4478 box->high.y, box->low.y);
4479
4480 PG_RETURN_POLYGON_P(poly);
4481 }
4482
4483
4484 Datum
poly_path(PG_FUNCTION_ARGS)4485 poly_path(PG_FUNCTION_ARGS)
4486 {
4487 POLYGON *poly = PG_GETARG_POLYGON_P(0);
4488 PATH *path;
4489 int size;
4490 int i;
4491
4492 /*
4493 * Never overflows: the old size fit in MaxAllocSize, and the new size is
4494 * smaller by a small constant.
4495 */
4496 size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts;
4497 path = (PATH *) palloc(size);
4498
4499 SET_VARSIZE(path, size);
4500 path->npts = poly->npts;
4501 path->closed = true;
4502 /* prevent instability in unused pad bytes */
4503 path->dummy = 0;
4504
4505 for (i = 0; i < poly->npts; i++)
4506 {
4507 path->p[i].x = poly->p[i].x;
4508 path->p[i].y = poly->p[i].y;
4509 }
4510
4511 PG_RETURN_PATH_P(path);
4512 }
4513
4514
4515 /***********************************************************************
4516 **
4517 ** Routines for circles.
4518 **
4519 ***********************************************************************/
4520
4521 /*----------------------------------------------------------
4522 * Formatting and conversion routines.
4523 *---------------------------------------------------------*/
4524
4525 /* circle_in - convert a string to internal form.
4526 *
4527 * External format: (center and radius of circle)
4528 * "<(f8,f8),f8>"
4529 * also supports quick entry style "f8,f8,f8"
4530 */
4531 Datum
circle_in(PG_FUNCTION_ARGS)4532 circle_in(PG_FUNCTION_ARGS)
4533 {
4534 char *str = PG_GETARG_CSTRING(0);
4535 CIRCLE *circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4536 char *s,
4537 *cp;
4538 int depth = 0;
4539
4540 s = str;
4541 while (isspace((unsigned char) *s))
4542 s++;
4543 if (*s == LDELIM_C)
4544 depth++, s++;
4545 else if (*s == LDELIM)
4546 {
4547 /* If there are two left parens, consume the first one */
4548 cp = (s + 1);
4549 while (isspace((unsigned char) *cp))
4550 cp++;
4551 if (*cp == LDELIM)
4552 depth++, s = cp;
4553 }
4554
4555 /* pair_decode will consume parens around the pair, if any */
4556 pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
4557
4558 if (*s == DELIM)
4559 s++;
4560
4561 circle->radius = single_decode(s, &s, "circle", str);
4562 if (circle->radius < 0)
4563 ereport(ERROR,
4564 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4565 errmsg("invalid input syntax for type %s: \"%s\"",
4566 "circle", str)));
4567
4568 while (depth > 0)
4569 {
4570 if ((*s == RDELIM)
4571 || ((*s == RDELIM_C) && (depth == 1)))
4572 {
4573 depth--;
4574 s++;
4575 while (isspace((unsigned char) *s))
4576 s++;
4577 }
4578 else
4579 ereport(ERROR,
4580 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4581 errmsg("invalid input syntax for type %s: \"%s\"",
4582 "circle", str)));
4583 }
4584
4585 if (*s != '\0')
4586 ereport(ERROR,
4587 (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4588 errmsg("invalid input syntax for type %s: \"%s\"",
4589 "circle", str)));
4590
4591 PG_RETURN_CIRCLE_P(circle);
4592 }
4593
4594 /* circle_out - convert a circle to external form.
4595 */
4596 Datum
circle_out(PG_FUNCTION_ARGS)4597 circle_out(PG_FUNCTION_ARGS)
4598 {
4599 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4600 StringInfoData str;
4601
4602 initStringInfo(&str);
4603
4604 appendStringInfoChar(&str, LDELIM_C);
4605 appendStringInfoChar(&str, LDELIM);
4606 pair_encode(circle->center.x, circle->center.y, &str);
4607 appendStringInfoChar(&str, RDELIM);
4608 appendStringInfoChar(&str, DELIM);
4609 single_encode(circle->radius, &str);
4610 appendStringInfoChar(&str, RDELIM_C);
4611
4612 PG_RETURN_CSTRING(str.data);
4613 }
4614
4615 /*
4616 * circle_recv - converts external binary format to circle
4617 */
4618 Datum
circle_recv(PG_FUNCTION_ARGS)4619 circle_recv(PG_FUNCTION_ARGS)
4620 {
4621 StringInfo buf = (StringInfo) PG_GETARG_POINTER(0);
4622 CIRCLE *circle;
4623
4624 circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4625
4626 circle->center.x = pq_getmsgfloat8(buf);
4627 circle->center.y = pq_getmsgfloat8(buf);
4628 circle->radius = pq_getmsgfloat8(buf);
4629
4630 if (circle->radius < 0)
4631 ereport(ERROR,
4632 (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
4633 errmsg("invalid radius in external \"circle\" value")));
4634
4635 PG_RETURN_CIRCLE_P(circle);
4636 }
4637
4638 /*
4639 * circle_send - converts circle to binary format
4640 */
4641 Datum
circle_send(PG_FUNCTION_ARGS)4642 circle_send(PG_FUNCTION_ARGS)
4643 {
4644 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4645 StringInfoData buf;
4646
4647 pq_begintypsend(&buf);
4648 pq_sendfloat8(&buf, circle->center.x);
4649 pq_sendfloat8(&buf, circle->center.y);
4650 pq_sendfloat8(&buf, circle->radius);
4651 PG_RETURN_BYTEA_P(pq_endtypsend(&buf));
4652 }
4653
4654
4655 /*----------------------------------------------------------
4656 * Relational operators for CIRCLEs.
4657 * <, >, <=, >=, and == are based on circle area.
4658 *---------------------------------------------------------*/
4659
4660 /* circles identical?
4661 */
4662 Datum
circle_same(PG_FUNCTION_ARGS)4663 circle_same(PG_FUNCTION_ARGS)
4664 {
4665 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4666 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4667
4668 PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
4669 FPeq(circle1->center.x, circle2->center.x) &&
4670 FPeq(circle1->center.y, circle2->center.y));
4671 }
4672
4673 /* circle_overlap - does circle1 overlap circle2?
4674 */
4675 Datum
circle_overlap(PG_FUNCTION_ARGS)4676 circle_overlap(PG_FUNCTION_ARGS)
4677 {
4678 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4679 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4680
4681 PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4682 circle1->radius + circle2->radius));
4683 }
4684
4685 /* circle_overleft - is the right edge of circle1 at or left of
4686 * the right edge of circle2?
4687 */
4688 Datum
circle_overleft(PG_FUNCTION_ARGS)4689 circle_overleft(PG_FUNCTION_ARGS)
4690 {
4691 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4692 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4693
4694 PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
4695 (circle2->center.x + circle2->radius)));
4696 }
4697
4698 /* circle_left - is circle1 strictly left of circle2?
4699 */
4700 Datum
circle_left(PG_FUNCTION_ARGS)4701 circle_left(PG_FUNCTION_ARGS)
4702 {
4703 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4704 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4705
4706 PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
4707 (circle2->center.x - circle2->radius)));
4708 }
4709
4710 /* circle_right - is circle1 strictly right of circle2?
4711 */
4712 Datum
circle_right(PG_FUNCTION_ARGS)4713 circle_right(PG_FUNCTION_ARGS)
4714 {
4715 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4716 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4717
4718 PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
4719 (circle2->center.x + circle2->radius)));
4720 }
4721
4722 /* circle_overright - is the left edge of circle1 at or right of
4723 * the left edge of circle2?
4724 */
4725 Datum
circle_overright(PG_FUNCTION_ARGS)4726 circle_overright(PG_FUNCTION_ARGS)
4727 {
4728 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4729 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4730
4731 PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
4732 (circle2->center.x - circle2->radius)));
4733 }
4734
4735 /* circle_contained - is circle1 contained by circle2?
4736 */
4737 Datum
circle_contained(PG_FUNCTION_ARGS)4738 circle_contained(PG_FUNCTION_ARGS)
4739 {
4740 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4741 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4742
4743 PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle1->radius), circle2->radius));
4744 }
4745
4746 /* circle_contain - does circle1 contain circle2?
4747 */
4748 Datum
circle_contain(PG_FUNCTION_ARGS)4749 circle_contain(PG_FUNCTION_ARGS)
4750 {
4751 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4752 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4753
4754 PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle2->radius), circle1->radius));
4755 }
4756
4757
4758 /* circle_below - is circle1 strictly below circle2?
4759 */
4760 Datum
circle_below(PG_FUNCTION_ARGS)4761 circle_below(PG_FUNCTION_ARGS)
4762 {
4763 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4764 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4765
4766 PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
4767 (circle2->center.y - circle2->radius)));
4768 }
4769
4770 /* circle_above - is circle1 strictly above circle2?
4771 */
4772 Datum
circle_above(PG_FUNCTION_ARGS)4773 circle_above(PG_FUNCTION_ARGS)
4774 {
4775 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4776 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4777
4778 PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
4779 (circle2->center.y + circle2->radius)));
4780 }
4781
4782 /* circle_overbelow - is the upper edge of circle1 at or below
4783 * the upper edge of circle2?
4784 */
4785 Datum
circle_overbelow(PG_FUNCTION_ARGS)4786 circle_overbelow(PG_FUNCTION_ARGS)
4787 {
4788 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4789 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4790
4791 PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
4792 (circle2->center.y + circle2->radius)));
4793 }
4794
4795 /* circle_overabove - is the lower edge of circle1 at or above
4796 * the lower edge of circle2?
4797 */
4798 Datum
circle_overabove(PG_FUNCTION_ARGS)4799 circle_overabove(PG_FUNCTION_ARGS)
4800 {
4801 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4802 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4803
4804 PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
4805 (circle2->center.y - circle2->radius)));
4806 }
4807
4808
4809 /* circle_relop - is area(circle1) relop area(circle2), within
4810 * our accuracy constraint?
4811 */
4812 Datum
circle_eq(PG_FUNCTION_ARGS)4813 circle_eq(PG_FUNCTION_ARGS)
4814 {
4815 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4816 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4817
4818 PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
4819 }
4820
4821 Datum
circle_ne(PG_FUNCTION_ARGS)4822 circle_ne(PG_FUNCTION_ARGS)
4823 {
4824 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4825 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4826
4827 PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
4828 }
4829
4830 Datum
circle_lt(PG_FUNCTION_ARGS)4831 circle_lt(PG_FUNCTION_ARGS)
4832 {
4833 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4834 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4835
4836 PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
4837 }
4838
4839 Datum
circle_gt(PG_FUNCTION_ARGS)4840 circle_gt(PG_FUNCTION_ARGS)
4841 {
4842 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4843 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4844
4845 PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
4846 }
4847
4848 Datum
circle_le(PG_FUNCTION_ARGS)4849 circle_le(PG_FUNCTION_ARGS)
4850 {
4851 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4852 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4853
4854 PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
4855 }
4856
4857 Datum
circle_ge(PG_FUNCTION_ARGS)4858 circle_ge(PG_FUNCTION_ARGS)
4859 {
4860 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4861 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4862
4863 PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
4864 }
4865
4866
4867 /*----------------------------------------------------------
4868 * "Arithmetic" operators on circles.
4869 *---------------------------------------------------------*/
4870
4871 static CIRCLE *
circle_copy(CIRCLE * circle)4872 circle_copy(CIRCLE *circle)
4873 {
4874 CIRCLE *result;
4875
4876 if (!PointerIsValid(circle))
4877 return NULL;
4878
4879 result = (CIRCLE *) palloc(sizeof(CIRCLE));
4880 memcpy((char *) result, (char *) circle, sizeof(CIRCLE));
4881 return result;
4882 }
4883
4884
4885 /* circle_add_pt()
4886 * Translation operator.
4887 */
4888 Datum
circle_add_pt(PG_FUNCTION_ARGS)4889 circle_add_pt(PG_FUNCTION_ARGS)
4890 {
4891 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4892 Point *point = PG_GETARG_POINT_P(1);
4893 CIRCLE *result;
4894
4895 result = circle_copy(circle);
4896
4897 result->center.x += point->x;
4898 result->center.y += point->y;
4899
4900 PG_RETURN_CIRCLE_P(result);
4901 }
4902
4903 Datum
circle_sub_pt(PG_FUNCTION_ARGS)4904 circle_sub_pt(PG_FUNCTION_ARGS)
4905 {
4906 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4907 Point *point = PG_GETARG_POINT_P(1);
4908 CIRCLE *result;
4909
4910 result = circle_copy(circle);
4911
4912 result->center.x -= point->x;
4913 result->center.y -= point->y;
4914
4915 PG_RETURN_CIRCLE_P(result);
4916 }
4917
4918
4919 /* circle_mul_pt()
4920 * Rotation and scaling operators.
4921 */
4922 Datum
circle_mul_pt(PG_FUNCTION_ARGS)4923 circle_mul_pt(PG_FUNCTION_ARGS)
4924 {
4925 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4926 Point *point = PG_GETARG_POINT_P(1);
4927 CIRCLE *result;
4928 Point *p;
4929
4930 result = circle_copy(circle);
4931
4932 p = DatumGetPointP(DirectFunctionCall2(point_mul,
4933 PointPGetDatum(&circle->center),
4934 PointPGetDatum(point)));
4935 result->center.x = p->x;
4936 result->center.y = p->y;
4937 result->radius *= HYPOT(point->x, point->y);
4938
4939 PG_RETURN_CIRCLE_P(result);
4940 }
4941
4942 Datum
circle_div_pt(PG_FUNCTION_ARGS)4943 circle_div_pt(PG_FUNCTION_ARGS)
4944 {
4945 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4946 Point *point = PG_GETARG_POINT_P(1);
4947 CIRCLE *result;
4948 Point *p;
4949
4950 result = circle_copy(circle);
4951
4952 p = DatumGetPointP(DirectFunctionCall2(point_div,
4953 PointPGetDatum(&circle->center),
4954 PointPGetDatum(point)));
4955 result->center.x = p->x;
4956 result->center.y = p->y;
4957 result->radius /= HYPOT(point->x, point->y);
4958
4959 PG_RETURN_CIRCLE_P(result);
4960 }
4961
4962
4963 /* circle_area - returns the area of the circle.
4964 */
4965 Datum
circle_area(PG_FUNCTION_ARGS)4966 circle_area(PG_FUNCTION_ARGS)
4967 {
4968 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4969
4970 PG_RETURN_FLOAT8(circle_ar(circle));
4971 }
4972
4973
4974 /* circle_diameter - returns the diameter of the circle.
4975 */
4976 Datum
circle_diameter(PG_FUNCTION_ARGS)4977 circle_diameter(PG_FUNCTION_ARGS)
4978 {
4979 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4980
4981 PG_RETURN_FLOAT8(2 * circle->radius);
4982 }
4983
4984
4985 /* circle_radius - returns the radius of the circle.
4986 */
4987 Datum
circle_radius(PG_FUNCTION_ARGS)4988 circle_radius(PG_FUNCTION_ARGS)
4989 {
4990 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4991
4992 PG_RETURN_FLOAT8(circle->radius);
4993 }
4994
4995
4996 /* circle_distance - returns the distance between
4997 * two circles.
4998 */
4999 Datum
circle_distance(PG_FUNCTION_ARGS)5000 circle_distance(PG_FUNCTION_ARGS)
5001 {
5002 CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
5003 CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
5004 float8 result;
5005
5006 result = point_dt(&circle1->center, &circle2->center)
5007 - (circle1->radius + circle2->radius);
5008 if (result < 0)
5009 result = 0;
5010 PG_RETURN_FLOAT8(result);
5011 }
5012
5013
5014 Datum
circle_contain_pt(PG_FUNCTION_ARGS)5015 circle_contain_pt(PG_FUNCTION_ARGS)
5016 {
5017 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5018 Point *point = PG_GETARG_POINT_P(1);
5019 double d;
5020
5021 d = point_dt(&circle->center, point);
5022 PG_RETURN_BOOL(d <= circle->radius);
5023 }
5024
5025
5026 Datum
pt_contained_circle(PG_FUNCTION_ARGS)5027 pt_contained_circle(PG_FUNCTION_ARGS)
5028 {
5029 Point *point = PG_GETARG_POINT_P(0);
5030 CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5031 double d;
5032
5033 d = point_dt(&circle->center, point);
5034 PG_RETURN_BOOL(d <= circle->radius);
5035 }
5036
5037
5038 /* dist_pc - returns the distance between
5039 * a point and a circle.
5040 */
5041 Datum
dist_pc(PG_FUNCTION_ARGS)5042 dist_pc(PG_FUNCTION_ARGS)
5043 {
5044 Point *point = PG_GETARG_POINT_P(0);
5045 CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5046 float8 result;
5047
5048 result = point_dt(point, &circle->center) - circle->radius;
5049 if (result < 0)
5050 result = 0;
5051 PG_RETURN_FLOAT8(result);
5052 }
5053
5054 /*
5055 * Distance from a circle to a point
5056 */
5057 Datum
dist_cpoint(PG_FUNCTION_ARGS)5058 dist_cpoint(PG_FUNCTION_ARGS)
5059 {
5060 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5061 Point *point = PG_GETARG_POINT_P(1);
5062 float8 result;
5063
5064 result = point_dt(point, &circle->center) - circle->radius;
5065 if (result < 0)
5066 result = 0;
5067 PG_RETURN_FLOAT8(result);
5068 }
5069
5070 /* circle_center - returns the center point of the circle.
5071 */
5072 Datum
circle_center(PG_FUNCTION_ARGS)5073 circle_center(PG_FUNCTION_ARGS)
5074 {
5075 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5076 Point *result;
5077
5078 result = (Point *) palloc(sizeof(Point));
5079 result->x = circle->center.x;
5080 result->y = circle->center.y;
5081
5082 PG_RETURN_POINT_P(result);
5083 }
5084
5085
5086 /* circle_ar - returns the area of the circle.
5087 */
5088 static double
circle_ar(CIRCLE * circle)5089 circle_ar(CIRCLE *circle)
5090 {
5091 return M_PI * (circle->radius * circle->radius);
5092 }
5093
5094
5095 /*----------------------------------------------------------
5096 * Conversion operators.
5097 *---------------------------------------------------------*/
5098
5099 Datum
cr_circle(PG_FUNCTION_ARGS)5100 cr_circle(PG_FUNCTION_ARGS)
5101 {
5102 Point *center = PG_GETARG_POINT_P(0);
5103 float8 radius = PG_GETARG_FLOAT8(1);
5104 CIRCLE *result;
5105
5106 result = (CIRCLE *) palloc(sizeof(CIRCLE));
5107
5108 result->center.x = center->x;
5109 result->center.y = center->y;
5110 result->radius = radius;
5111
5112 PG_RETURN_CIRCLE_P(result);
5113 }
5114
5115 Datum
circle_box(PG_FUNCTION_ARGS)5116 circle_box(PG_FUNCTION_ARGS)
5117 {
5118 CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5119 BOX *box;
5120 double delta;
5121
5122 box = (BOX *) palloc(sizeof(BOX));
5123
5124 delta = circle->radius / sqrt(2.0);
5125
5126 box->high.x = circle->center.x + delta;
5127 box->low.x = circle->center.x - delta;
5128 box->high.y = circle->center.y + delta;
5129 box->low.y = circle->center.y - delta;
5130
5131 PG_RETURN_BOX_P(box);
5132 }
5133
5134 /* box_circle()
5135 * Convert a box to a circle.
5136 */
5137 Datum
box_circle(PG_FUNCTION_ARGS)5138 box_circle(PG_FUNCTION_ARGS)
5139 {
5140 BOX *box = PG_GETARG_BOX_P(0);
5141 CIRCLE *circle;
5142
5143 circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5144
5145 circle->center.x = (box->high.x + box->low.x) / 2;
5146 circle->center.y = (box->high.y + box->low.y) / 2;
5147
5148 circle->radius = point_dt(&circle->center, &box->high);
5149
5150 PG_RETURN_CIRCLE_P(circle);
5151 }
5152
5153
5154 Datum
circle_poly(PG_FUNCTION_ARGS)5155 circle_poly(PG_FUNCTION_ARGS)
5156 {
5157 int32 npts = PG_GETARG_INT32(0);
5158 CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5159 POLYGON *poly;
5160 int base_size,
5161 size;
5162 int i;
5163 double angle;
5164 double anglestep;
5165
5166 if (FPzero(circle->radius))
5167 ereport(ERROR,
5168 (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
5169 errmsg("cannot convert circle with radius zero to polygon")));
5170
5171 if (npts < 2)
5172 ereport(ERROR,
5173 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5174 errmsg("must request at least 2 points")));
5175
5176 base_size = sizeof(poly->p[0]) * npts;
5177 size = offsetof(POLYGON, p) + base_size;
5178
5179 /* Check for integer overflow */
5180 if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
5181 ereport(ERROR,
5182 (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
5183 errmsg("too many points requested")));
5184
5185 poly = (POLYGON *) palloc0(size); /* zero any holes */
5186 SET_VARSIZE(poly, size);
5187 poly->npts = npts;
5188
5189 anglestep = (2.0 * M_PI) / npts;
5190
5191 for (i = 0; i < npts; i++)
5192 {
5193 angle = i * anglestep;
5194 poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
5195 poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
5196 }
5197
5198 make_bound_box(poly);
5199
5200 PG_RETURN_POLYGON_P(poly);
5201 }
5202
5203 /* poly_circle - convert polygon to circle
5204 *
5205 * XXX This algorithm should use weighted means of line segments
5206 * rather than straight average values of points - tgl 97/01/21.
5207 */
5208 Datum
poly_circle(PG_FUNCTION_ARGS)5209 poly_circle(PG_FUNCTION_ARGS)
5210 {
5211 POLYGON *poly = PG_GETARG_POLYGON_P(0);
5212 CIRCLE *circle;
5213 int i;
5214
5215 if (poly->npts < 1)
5216 ereport(ERROR,
5217 (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5218 errmsg("cannot convert empty polygon to circle")));
5219
5220 circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5221
5222 circle->center.x = 0;
5223 circle->center.y = 0;
5224 circle->radius = 0;
5225
5226 for (i = 0; i < poly->npts; i++)
5227 {
5228 circle->center.x += poly->p[i].x;
5229 circle->center.y += poly->p[i].y;
5230 }
5231 circle->center.x /= poly->npts;
5232 circle->center.y /= poly->npts;
5233
5234 for (i = 0; i < poly->npts; i++)
5235 circle->radius += point_dt(&poly->p[i], &circle->center);
5236 circle->radius /= poly->npts;
5237
5238 PG_RETURN_CIRCLE_P(circle);
5239 }
5240
5241
5242 /***********************************************************************
5243 **
5244 ** Private routines for multiple types.
5245 **
5246 ***********************************************************************/
5247
5248 /*
5249 * Test to see if the point is inside the polygon, returns 1/0, or 2 if
5250 * the point is on the polygon.
5251 * Code adapted but not copied from integer-based routines in WN: A
5252 * Server for the HTTP
5253 * version 1.15.1, file wn/image.c
5254 * http://hopf.math.northwestern.edu/index.html
5255 * Description of algorithm: http://www.linuxjournal.com/article/2197
5256 * http://www.linuxjournal.com/article/2029
5257 */
5258
5259 #define POINT_ON_POLYGON INT_MAX
5260
5261 static int
point_inside(Point * p,int npts,Point * plist)5262 point_inside(Point *p, int npts, Point *plist)
5263 {
5264 double x0,
5265 y0;
5266 double prev_x,
5267 prev_y;
5268 int i = 0;
5269 double x,
5270 y;
5271 int cross,
5272 total_cross = 0;
5273
5274 if (npts <= 0)
5275 return 0;
5276
5277 /* compute first polygon point relative to single point */
5278 x0 = plist[0].x - p->x;
5279 y0 = plist[0].y - p->y;
5280
5281 prev_x = x0;
5282 prev_y = y0;
5283 /* loop over polygon points and aggregate total_cross */
5284 for (i = 1; i < npts; i++)
5285 {
5286 /* compute next polygon point relative to single point */
5287 x = plist[i].x - p->x;
5288 y = plist[i].y - p->y;
5289
5290 /* compute previous to current point crossing */
5291 if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
5292 return 2;
5293 total_cross += cross;
5294
5295 prev_x = x;
5296 prev_y = y;
5297 }
5298
5299 /* now do the first point */
5300 if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
5301 return 2;
5302 total_cross += cross;
5303
5304 if (total_cross != 0)
5305 return 1;
5306 return 0;
5307 }
5308
5309
5310 /* lseg_crossing()
5311 * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
5312 * Returns +/-1 if one point is on the positive X-axis.
5313 * Returns 0 if both points are on the positive X-axis, or there is no crossing.
5314 * Returns POINT_ON_POLYGON if the segment contains (0,0).
5315 * Wow, that is one confusing API, but it is used above, and when summed,
5316 * can tell is if a point is in a polygon.
5317 */
5318
5319 static int
lseg_crossing(double x,double y,double prev_x,double prev_y)5320 lseg_crossing(double x, double y, double prev_x, double prev_y)
5321 {
5322 double z;
5323 int y_sign;
5324
5325 if (FPzero(y))
5326 { /* y == 0, on X axis */
5327 if (FPzero(x)) /* (x,y) is (0,0)? */
5328 return POINT_ON_POLYGON;
5329 else if (FPgt(x, 0))
5330 { /* x > 0 */
5331 if (FPzero(prev_y)) /* y and prev_y are zero */
5332 /* prev_x > 0? */
5333 return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5334 return FPlt(prev_y, 0) ? 1 : -1;
5335 }
5336 else
5337 { /* x < 0, x not on positive X axis */
5338 if (FPzero(prev_y))
5339 /* prev_x < 0? */
5340 return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5341 return 0;
5342 }
5343 }
5344 else
5345 { /* y != 0 */
5346 /* compute y crossing direction from previous point */
5347 y_sign = FPgt(y, 0) ? 1 : -1;
5348
5349 if (FPzero(prev_y))
5350 /* previous point was on X axis, so new point is either off or on */
5351 return FPlt(prev_x, 0) ? 0 : y_sign;
5352 else if (FPgt(y_sign * prev_y, 0))
5353 /* both above or below X axis */
5354 return 0; /* same sign */
5355 else
5356 { /* y and prev_y cross X-axis */
5357 if (FPge(x, 0) && FPgt(prev_x, 0))
5358 /* both non-negative so cross positive X-axis */
5359 return 2 * y_sign;
5360 if (FPlt(x, 0) && FPle(prev_x, 0))
5361 /* both non-positive so do not cross positive X-axis */
5362 return 0;
5363
5364 /* x and y cross axises, see URL above point_inside() */
5365 z = (x - prev_x) * y - (y - prev_y) * x;
5366 if (FPzero(z))
5367 return POINT_ON_POLYGON;
5368 return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
5369 }
5370 }
5371 }
5372
5373
5374 static bool
plist_same(int npts,Point * p1,Point * p2)5375 plist_same(int npts, Point *p1, Point *p2)
5376 {
5377 int i,
5378 ii,
5379 j;
5380
5381 /* find match for first point */
5382 for (i = 0; i < npts; i++)
5383 {
5384 if ((FPeq(p2[i].x, p1[0].x))
5385 && (FPeq(p2[i].y, p1[0].y)))
5386 {
5387
5388 /* match found? then look forward through remaining points */
5389 for (ii = 1, j = i + 1; ii < npts; ii++, j++)
5390 {
5391 if (j >= npts)
5392 j = 0;
5393 if ((!FPeq(p2[j].x, p1[ii].x))
5394 || (!FPeq(p2[j].y, p1[ii].y)))
5395 {
5396 #ifdef GEODEBUG
5397 printf("plist_same- %d failed forward match with %d\n", j, ii);
5398 #endif
5399 break;
5400 }
5401 }
5402 #ifdef GEODEBUG
5403 printf("plist_same- ii = %d/%d after forward match\n", ii, npts);
5404 #endif
5405 if (ii == npts)
5406 return true;
5407
5408 /* match not found forwards? then look backwards */
5409 for (ii = 1, j = i - 1; ii < npts; ii++, j--)
5410 {
5411 if (j < 0)
5412 j = (npts - 1);
5413 if ((!FPeq(p2[j].x, p1[ii].x))
5414 || (!FPeq(p2[j].y, p1[ii].y)))
5415 {
5416 #ifdef GEODEBUG
5417 printf("plist_same- %d failed reverse match with %d\n", j, ii);
5418 #endif
5419 break;
5420 }
5421 }
5422 #ifdef GEODEBUG
5423 printf("plist_same- ii = %d/%d after reverse match\n", ii, npts);
5424 #endif
5425 if (ii == npts)
5426 return true;
5427 }
5428 }
5429
5430 return false;
5431 }
5432
5433
5434 /*-------------------------------------------------------------------------
5435 * Determine the hypotenuse.
5436 *
5437 * If required, x and y are swapped to make x the larger number. The
5438 * traditional formula of x^2+y^2 is rearranged to factor x outside the
5439 * sqrt. This allows computation of the hypotenuse for significantly
5440 * larger values, and with a higher precision than when using the naive
5441 * formula. In particular, this cannot overflow unless the final result
5442 * would be out-of-range.
5443 *
5444 * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
5445 * = x * sqrt( 1 + y^2/x^2 )
5446 * = x * sqrt( 1 + y/x * y/x )
5447 *
5448 * It is expected that this routine will eventually be replaced with the
5449 * C99 hypot() function.
5450 *
5451 * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
5452 * case of hypot(inf,nan) results in INF, and not NAN.
5453 *-----------------------------------------------------------------------
5454 */
5455 double
pg_hypot(double x,double y)5456 pg_hypot(double x, double y)
5457 {
5458 double yx;
5459
5460 /* Handle INF and NaN properly */
5461 if (isinf(x) || isinf(y))
5462 return get_float8_infinity();
5463
5464 if (isnan(x) || isnan(y))
5465 return get_float8_nan();
5466
5467 /* Else, drop any minus signs */
5468 x = fabs(x);
5469 y = fabs(y);
5470
5471 /* Swap x and y if needed to make x the larger one */
5472 if (x < y)
5473 {
5474 double temp = x;
5475
5476 x = y;
5477 y = temp;
5478 }
5479
5480 /*
5481 * If y is zero, the hypotenuse is x. This test saves a few cycles in
5482 * such cases, but more importantly it also protects against
5483 * divide-by-zero errors, since now x >= y.
5484 */
5485 if (y == 0.0)
5486 return x;
5487
5488 /* Determine the hypotenuse */
5489 yx = y / x;
5490 return x * sqrt(1.0 + (yx * yx));
5491 }
5492