PostgreSQL Source Code  git master
geo_ops.c
Go to the documentation of this file.
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 
37 {
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
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
132 {
133  char *xstr = float8out_internal(x);
134 
135  appendStringInfoString(str, xstr);
136  pfree(xstr);
137 } /* single_encode() */
138 
139 static void
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
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
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 *
270 {
272  int i;
273 
274  initStringInfo(&str);
275 
276  switch (path_delim)
277  {
278  case PATH_CLOSED:
280  break;
281  case PATH_OPEN:
283  break;
284  case PATH_NONE:
285  break;
286  }
287 
288  for (i = 0; i < npts; i++)
289  {
290  if (i > 0)
293  pair_encode(pt->x, pt->y, &str);
295  pt++;
296  }
297 
298  switch (path_delim)
299  {
300  case PATH_CLOSED:
302  break;
303  case PATH_OPEN:
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
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
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
382 {
383  BOX *box = PG_GETARG_BOX_P(0);
384 
386 }
387 
388 /*
389  * box_recv - converts external binary format to box
390  */
391 Datum
393 {
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
428 {
429  BOX *box = PG_GETARG_BOX_P(0);
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);
438 }
439 
440 
441 /* box_construct - fill in a new box.
442  */
443 static BOX *
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 *
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 *
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
971 {
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
989 {
990  LINE *line = PG_GETARG_LINE_P(0);
992 
993  pq_begintypsend(&buf);
994  pq_sendfloat8(&buf, line->A);
995  pq_sendfloat8(&buf, line->B);
996  pq_sendfloat8(&buf, line->C);
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 *
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
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
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
1095 {
1096  LINE *l1 = PG_GETARG_LINE_P(0);
1097  LINE *l2 = PG_GETARG_LINE_P(1);
1098 
1100  LinePGetDatum(l1),
1101  LinePGetDatum(l2))));
1102 }
1103 
1104 Datum
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
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
1132 {
1133  LINE *line = PG_GETARG_LINE_P(0);
1134 
1135  PG_RETURN_BOOL(FPzero(line->B));
1136 }
1137 
1138 Datum
1140 {
1141  LINE *line = PG_GETARG_LINE_P(0);
1142 
1143  PG_RETURN_BOOL(FPzero(line->A));
1144 }
1145 
1146 Datum
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
1177 {
1178  LINE *l1 = PG_GETARG_LINE_P(0);
1179  LINE *l2 = PG_GETARG_LINE_P(1);
1180  float8 result;
1181  Point *tmp;
1182 
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
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 *
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  */
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
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
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
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
1390 {
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
1428 {
1429  PATH *path = PG_GETARG_PATH_P(0);
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  }
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
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
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
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
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
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
1504 {
1505  PATH *path = PG_GETARG_PATH_P(0);
1506 
1507  PG_RETURN_BOOL(path->closed);
1508 }
1509 
1510 Datum
1512 {
1513  PATH *path = PG_GETARG_PATH_P(0);
1514 
1515  PG_RETURN_BOOL(!path->closed);
1516 }
1517 
1518 Datum
1520 {
1521  PATH *path = PG_GETARG_PATH_P(0);
1522 
1523  PG_RETURN_INT32(path->npts);
1524 }
1525 
1526 
1527 Datum
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
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
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
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 
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
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
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
1746 {
1747  Point *pt = PG_GETARG_POINT_P(0);
1748 
1750 }
1751 
1752 /*
1753  * point_recv - converts external binary format to point
1754  */
1755 Datum
1757 {
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
1772 {
1773  Point *pt = PG_GETARG_POINT_P(0);
1775 
1776  pq_begintypsend(&buf);
1777  pq_sendfloat8(&buf, pt->x);
1778  pq_sendfloat8(&buf, pt->y);
1780 }
1781 
1782 
1783 static Point *
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 *
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
1973 {
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
1992 {
1993  LSEG *ls = PG_GETARG_LSEG_P(0);
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);
2002 }
2003 
2004 
2005 /* lseg_construct -
2006  * form a LSEG from two Points.
2007  */
2008 Datum
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
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
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
2051 {
2052  LSEG *l1 = PG_GETARG_LSEG_P(0);
2053  LSEG *l2 = PG_GETARG_LSEG_P(1);
2054 
2056 }
2057 
2058 static bool
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
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
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
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
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
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
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
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
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
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
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
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
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
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 *
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
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
2339 {
2340  Point *pt = PG_GETARG_POINT_P(0);
2341  LINE *line = PG_GETARG_LINE_P(1);
2342 
2344 }
2345 
2346 static double
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
2358 {
2359  Point *pt = PG_GETARG_POINT_P(0);
2360  LSEG *lseg = PG_GETARG_LSEG_P(1);
2361 
2363 }
2364 
2365 static double
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
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
2480 {
2481  Point *pt = PG_GETARG_POINT_P(0);
2482  BOX *box = PG_GETARG_BOX_P(1);
2483  float8 result;
2484  Point *near;
2485 
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
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
2524 {
2525  LSEG *lseg = PG_GETARG_LSEG_P(0);
2526  BOX *box = PG_GETARG_BOX_P(1);
2527  Point *tmp;
2528  Datum result;
2529 
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
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
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
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
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
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 *
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
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
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
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
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  {
2901  PointPGetDatum(&l2->p[0]),
2902  LsegPGetDatum(l1)));
2903  memcpy(&point, result, sizeof(Point));
2905  PointPGetDatum(&point),
2906  LsegPGetDatum(l2)));
2907  }
2908 
2909  if (dist_ps_internal(&l2->p[1], l1) < dist)
2910  {
2912  PointPGetDatum(&l2->p[1]),
2913  LsegPGetDatum(l1)));
2914  memcpy(&point, result, sizeof(Point));
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
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 
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 
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
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
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
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 */
3061  LsegPGetDatum(lseg),
3062  BoxPGetDatum(box))))
3063  {
3064  box_cn(&point, box);
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 */
3101  LsegPGetDatum(lseg),
3102  LsegPGetDatum(&bseg)));
3103 }
3104 
3105 Datum
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
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
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
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
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
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
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
3222 {
3223  LSEG *lseg = PG_GETARG_LSEG_P(0);
3224  LINE *line = PG_GETARG_LINE_P(1);
3225 
3227  PointPGetDatum(&lseg->p[0]),
3228  LinePGetDatum(line))) &&
3230  PointPGetDatum(&lseg->p[1]),
3231  LinePGetDatum(line))));
3232 }
3233 
3234 Datum
3236 {
3237  LSEG *lseg = PG_GETARG_LSEG_P(0);
3238  BOX *box = PG_GETARG_BOX_P(1);
3239 
3241  PointPGetDatum(&lseg->p[0]),
3242  BoxPGetDatum(box))) &&
3244  PointPGetDatum(&lseg->p[1]),
3245  BoxPGetDatum(box))));
3246 }
3247 
3248 /*---------------------------------------------------------------------
3249  * inter_
3250  * Whether one object intersects another.
3251  *-------------------------------------------------------------------*/
3252 
3253 Datum
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
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 */
3292  PointPGetDatum(&lseg->p[0]),
3293  BoxPGetDatum(box))) ||
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
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
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
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
3459 {
3460  POLYGON *poly = PG_GETARG_POLYGON_P(0);
3461 
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
3475 {
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
3510 {
3511  POLYGON *poly = PG_GETARG_POLYGON_P(0);
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  }
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
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
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
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
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
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
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
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
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
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
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
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
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 
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
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 &&
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
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 */
3979 }
3980 
3981 
3982 Datum
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
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
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
4025 {
4026  float8 x = PG_GETARG_FLOAT8(0);
4027  float8 y = PG_GETARG_FLOAT8(1);
4028 
4030 }
4031 
4032 Datum
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
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
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
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
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
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
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
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 
4150  PointPGetDatum(&box->high),
4151  PointPGetDatum(p)));
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
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 
4171  PointPGetDatum(&box->high),
4172  PointPGetDatum(p)));
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
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
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
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
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
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
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  {
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
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  {
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
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
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
4413 {
4414  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4415 
4416  PG_RETURN_INT32(poly->npts);
4417 }
4418 
4419 
4420 Datum
4422 {
4423  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4424  Datum result;
4425  CIRCLE *circle;
4426 
4428  PolygonPGetDatum(poly)));
4430  CirclePGetDatum(circle));
4431 
4432  PG_RETURN_DATUM(result);
4433 }
4434 
4435 
4436 Datum
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
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
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
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) || (*s == LDELIM))
4544  {
4545  depth++;
4546  cp = (s + 1);
4547  while (isspace((unsigned char) *cp))
4548  cp++;
4549  if (*cp == LDELIM)
4550  s = cp;
4551  }
4552 
4553  pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
4554 
4555  if (*s == DELIM)
4556  s++;
4557 
4558  circle->radius = single_decode(s, &s, "circle", str);
4559  if (circle->radius < 0)
4560  ereport(ERROR,
4561  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4562  errmsg("invalid input syntax for type %s: \"%s\"",
4563  "circle", str)));
4564 
4565  while (depth > 0)
4566  {
4567  if ((*s == RDELIM)
4568  || ((*s == RDELIM_C) && (depth == 1)))
4569  {
4570  depth--;
4571  s++;
4572  while (isspace((unsigned char) *s))
4573  s++;
4574  }
4575  else
4576  ereport(ERROR,
4577  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4578  errmsg("invalid input syntax for type %s: \"%s\"",
4579  "circle", str)));
4580  }
4581 
4582  if (*s != '\0')
4583  ereport(ERROR,
4584  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4585  errmsg("invalid input syntax for type %s: \"%s\"",
4586  "circle", str)));
4587 
4588  PG_RETURN_CIRCLE_P(circle);
4589 }
4590 
4591 /* circle_out - convert a circle to external form.
4592  */
4593 Datum
4595 {
4596  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4598 
4599  initStringInfo(&str);
4600 
4603  pair_encode(circle->center.x, circle->center.y, &str);
4605  appendStringInfoChar(&str, DELIM);
4606  single_encode(circle->radius, &str);
4608 
4609  PG_RETURN_CSTRING(str.data);
4610 }
4611 
4612 /*
4613  * circle_recv - converts external binary format to circle
4614  */
4615 Datum
4617 {
4619  CIRCLE *circle;
4620 
4621  circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4622 
4623  circle->center.x = pq_getmsgfloat8(buf);
4624  circle->center.y = pq_getmsgfloat8(buf);
4625  circle->radius = pq_getmsgfloat8(buf);
4626 
4627  if (circle->radius < 0)
4628  ereport(ERROR,
4629  (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
4630  errmsg("invalid radius in external \"circle\" value")));
4631 
4632  PG_RETURN_CIRCLE_P(circle);
4633 }
4634 
4635 /*
4636  * circle_send - converts circle to binary format
4637  */
4638 Datum
4640 {
4641  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4643 
4644  pq_begintypsend(&buf);
4645  pq_sendfloat8(&buf, circle->center.x);
4646  pq_sendfloat8(&buf, circle->center.y);
4647  pq_sendfloat8(&buf, circle->radius);
4649 }
4650 
4651 
4652 /*----------------------------------------------------------
4653  * Relational operators for CIRCLEs.
4654  * <, >, <=, >=, and == are based on circle area.
4655  *---------------------------------------------------------*/
4656 
4657 /* circles identical?
4658  */
4659 Datum
4661 {
4662  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4663  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4664 
4665  PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) &&
4666  FPeq(circle1->center.x, circle2->center.x) &&
4667  FPeq(circle1->center.y, circle2->center.y));
4668 }
4669 
4670 /* circle_overlap - does circle1 overlap circle2?
4671  */
4672 Datum
4674 {
4675  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4676  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4677 
4678  PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4679  circle1->radius + circle2->radius));
4680 }
4681 
4682 /* circle_overleft - is the right edge of circle1 at or left of
4683  * the right edge of circle2?
4684  */
4685 Datum
4687 {
4688  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4689  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4690 
4691  PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius),
4692  (circle2->center.x + circle2->radius)));
4693 }
4694 
4695 /* circle_left - is circle1 strictly left of circle2?
4696  */
4697 Datum
4699 {
4700  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4701  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4702 
4703  PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius),
4704  (circle2->center.x - circle2->radius)));
4705 }
4706 
4707 /* circle_right - is circle1 strictly right of circle2?
4708  */
4709 Datum
4711 {
4712  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4713  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4714 
4715  PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius),
4716  (circle2->center.x + circle2->radius)));
4717 }
4718 
4719 /* circle_overright - is the left edge of circle1 at or right of
4720  * the left edge of circle2?
4721  */
4722 Datum
4724 {
4725  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4726  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4727 
4728  PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius),
4729  (circle2->center.x - circle2->radius)));
4730 }
4731 
4732 /* circle_contained - is circle1 contained by circle2?
4733  */
4734 Datum
4736 {
4737  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4738  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4739 
4740  PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle1->radius), circle2->radius));
4741 }
4742 
4743 /* circle_contain - does circle1 contain circle2?
4744  */
4745 Datum
4747 {
4748  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4749  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4750 
4751  PG_RETURN_BOOL(FPle((point_dt(&circle1->center, &circle2->center) + circle2->radius), circle1->radius));
4752 }
4753 
4754 
4755 /* circle_below - is circle1 strictly below circle2?
4756  */
4757 Datum
4759 {
4760  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4761  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4762 
4763  PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius),
4764  (circle2->center.y - circle2->radius)));
4765 }
4766 
4767 /* circle_above - is circle1 strictly above circle2?
4768  */
4769 Datum
4771 {
4772  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4773  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4774 
4775  PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius),
4776  (circle2->center.y + circle2->radius)));
4777 }
4778 
4779 /* circle_overbelow - is the upper edge of circle1 at or below
4780  * the upper edge of circle2?
4781  */
4782 Datum
4784 {
4785  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4786  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4787 
4788  PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius),
4789  (circle2->center.y + circle2->radius)));
4790 }
4791 
4792 /* circle_overabove - is the lower edge of circle1 at or above
4793  * the lower edge of circle2?
4794  */
4795 Datum
4797 {
4798  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4799  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4800 
4801  PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius),
4802  (circle2->center.y - circle2->radius)));
4803 }
4804 
4805 
4806 /* circle_relop - is area(circle1) relop area(circle2), within
4807  * our accuracy constraint?
4808  */
4809 Datum
4811 {
4812  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4813  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4814 
4815  PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
4816 }
4817 
4818 Datum
4820 {
4821  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4822  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4823 
4824  PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
4825 }
4826 
4827 Datum
4829 {
4830  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4831  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4832 
4833  PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
4834 }
4835 
4836 Datum
4838 {
4839  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4840  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4841 
4842  PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
4843 }
4844 
4845 Datum
4847 {
4848  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4849  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4850 
4851  PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
4852 }
4853 
4854 Datum
4856 {
4857  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4858  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4859 
4860  PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
4861 }
4862 
4863 
4864 /*----------------------------------------------------------
4865  * "Arithmetic" operators on circles.
4866  *---------------------------------------------------------*/
4867 
4868 static CIRCLE *
4870 {
4871  CIRCLE *result;
4872 
4873  if (!PointerIsValid(circle))
4874  return NULL;
4875 
4876  result = (CIRCLE *) palloc(sizeof(CIRCLE));
4877  memcpy((char *) result, (char *) circle, sizeof(CIRCLE));
4878  return result;
4879 }
4880 
4881 
4882 /* circle_add_pt()
4883  * Translation operator.
4884  */
4885 Datum
4887 {
4888  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4889  Point *point = PG_GETARG_POINT_P(1);
4890  CIRCLE *result;
4891 
4892  result = circle_copy(circle);
4893 
4894  result->center.x += point->x;
4895  result->center.y += point->y;
4896 
4897  PG_RETURN_CIRCLE_P(result);
4898 }
4899 
4900 Datum
4902 {
4903  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4904  Point *point = PG_GETARG_POINT_P(1);
4905  CIRCLE *result;
4906 
4907  result = circle_copy(circle);
4908 
4909  result->center.x -= point->x;
4910  result->center.y -= point->y;
4911 
4912  PG_RETURN_CIRCLE_P(result);
4913 }
4914 
4915 
4916 /* circle_mul_pt()
4917  * Rotation and scaling operators.
4918  */
4919 Datum
4921 {
4922  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4923  Point *point = PG_GETARG_POINT_P(1);
4924  CIRCLE *result;
4925  Point *p;
4926 
4927  result = circle_copy(circle);
4928 
4930  PointPGetDatum(&circle->center),
4931  PointPGetDatum(point)));
4932  result->center.x = p->x;
4933  result->center.y = p->y;
4934  result->radius *= HYPOT(point->x, point->y);
4935 
4936  PG_RETURN_CIRCLE_P(result);
4937 }
4938 
4939 Datum
4941 {
4942  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4943  Point *point = PG_GETARG_POINT_P(1);
4944  CIRCLE *result;
4945  Point *p;
4946 
4947  result = circle_copy(circle);
4948 
4950  PointPGetDatum(&circle->center),
4951  PointPGetDatum(point)));
4952  result->center.x = p->x;
4953  result->center.y = p->y;
4954  result->radius /= HYPOT(point->x, point->y);
4955 
4956  PG_RETURN_CIRCLE_P(result);
4957 }
4958 
4959 
4960 /* circle_area - returns the area of the circle.
4961  */
4962 Datum
4964 {
4965  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4966 
4967  PG_RETURN_FLOAT8(circle_ar(circle));
4968 }
4969 
4970 
4971 /* circle_diameter - returns the diameter of the circle.
4972  */
4973 Datum
4975 {
4976  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4977 
4978  PG_RETURN_FLOAT8(2 * circle->radius);
4979 }
4980 
4981 
4982 /* circle_radius - returns the radius of the circle.
4983  */
4984 Datum
4986 {
4987  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4988 
4989  PG_RETURN_FLOAT8(circle->radius);
4990 }
4991 
4992 
4993 /* circle_distance - returns the distance between
4994  * two circles.
4995  */
4996 Datum
4998 {
4999  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
5000  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
5001  float8 result;
5002 
5003  result = point_dt(&circle1->center, &circle2->center)
5004  - (circle1->radius + circle2->radius);
5005  if (result < 0)
5006  result = 0;
5007  PG_RETURN_FLOAT8(result);
5008 }
5009 
5010 
5011 Datum
5013 {
5014  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5015  Point *point = PG_GETARG_POINT_P(1);
5016  double d;
5017 
5018  d = point_dt(&circle->center, point);
5019  PG_RETURN_BOOL(d <= circle->radius);
5020 }
5021 
5022 
5023 Datum
5025 {
5026  Point *point = PG_GETARG_POINT_P(0);
5027  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5028  double d;
5029 
5030  d = point_dt(&circle->center, point);
5031  PG_RETURN_BOOL(d <= circle->radius);
5032 }
5033 
5034 
5035 /* dist_pc - returns the distance between
5036  * a point and a circle.
5037  */
5038 Datum
5040 {
5041  Point *point = PG_GETARG_POINT_P(0);
5042  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5043  float8 result;
5044 
5045  result = point_dt(point, &circle->center) - circle->radius;
5046  if (result < 0)
5047  result = 0;
5048  PG_RETURN_FLOAT8(result);
5049 }
5050 
5051 /*
5052  * Distance from a circle to a point
5053  */
5054 Datum
5056 {
5057  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5058  Point *point = PG_GETARG_POINT_P(1);
5059  float8 result;
5060 
5061  result = point_dt(point, &circle->center) - circle->radius;
5062  if (result < 0)
5063  result = 0;
5064  PG_RETURN_FLOAT8(result);
5065 }
5066 
5067 /* circle_center - returns the center point of the circle.
5068  */
5069 Datum
5071 {
5072  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5073  Point *result;
5074 
5075  result = (Point *) palloc(sizeof(Point));
5076  result->x = circle->center.x;
5077  result->y = circle->center.y;
5078 
5079  PG_RETURN_POINT_P(result);
5080 }
5081 
5082 
5083 /* circle_ar - returns the area of the circle.
5084  */
5085 static double
5087 {
5088  return M_PI * (circle->radius * circle->radius);
5089 }
5090 
5091 
5092 /*----------------------------------------------------------
5093  * Conversion operators.
5094  *---------------------------------------------------------*/
5095 
5096 Datum
5098 {
5099  Point *center = PG_GETARG_POINT_P(0);
5100  float8 radius = PG_GETARG_FLOAT8(1);
5101  CIRCLE *result;
5102 
5103  result = (CIRCLE *) palloc(sizeof(CIRCLE));
5104 
5105  result->center.x = center->x;
5106  result->center.y = center->y;
5107  result->radius = radius;
5108 
5109  PG_RETURN_CIRCLE_P(result);
5110 }
5111 
5112 Datum
5114 {
5115  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5116  BOX *box;
5117  double delta;
5118 
5119  box = (BOX *) palloc(sizeof(BOX));
5120 
5121  delta = circle->radius / sqrt(2.0);
5122 
5123  box->high.x = circle->center.x + delta;
5124  box->low.x = circle->center.x - delta;
5125  box->high.y = circle->center.y + delta;
5126  box->low.y = circle->center.y - delta;
5127 
5128  PG_RETURN_BOX_P(box);
5129 }
5130 
5131 /* box_circle()
5132  * Convert a box to a circle.
5133  */
5134 Datum
5136 {
5137  BOX *box = PG_GETARG_BOX_P(0);
5138  CIRCLE *circle;
5139 
5140  circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5141 
5142  circle->center.x = (box->high.x + box->low.x) / 2;
5143  circle->center.y = (box->high.y + box->low.y) / 2;
5144 
5145  circle->radius = point_dt(&circle->center, &box->high);
5146 
5147  PG_RETURN_CIRCLE_P(circle);
5148 }
5149 
5150 
5151 Datum
5153 {
5154  int32 npts = PG_GETARG_INT32(0);
5155  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5156  POLYGON *poly;
5157  int base_size,
5158  size;
5159  int i;
5160  double angle;
5161  double anglestep;
5162 
5163  if (FPzero(circle->radius))
5164  ereport(ERROR,
5165  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
5166  errmsg("cannot convert circle with radius zero to polygon")));
5167 
5168  if (npts < 2)
5169  ereport(ERROR,
5170  (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5171  errmsg("must request at least 2 points")));
5172 
5173  base_size = sizeof(poly->p[0]) * npts;
5174  size = offsetof(POLYGON, p) + base_size;
5175 
5176  /* Check for integer overflow */
5177  if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
5178  ereport(ERROR,
5179  (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
5180  errmsg("too many points requested")));
5181 
5182  poly = (POLYGON *) palloc0(size); /* zero any holes */
5183  SET_VARSIZE(poly, size);
5184  poly->npts = npts;
5185 
5186  anglestep = (2.0 * M_PI) / npts;
5187 
5188  for (i = 0; i < npts; i++)
5189  {
5190  angle = i * anglestep;
5191  poly->p[i].x = circle->center.x - (circle->radius * cos(angle));
5192  poly->p[i].y = circle->center.y + (circle->radius * sin(angle));
5193  }
5194 
5195  make_bound_box(poly);
5196 
5197  PG_RETURN_POLYGON_P(poly);
5198 }
5199 
5200 /* poly_circle - convert polygon to circle
5201  *
5202  * XXX This algorithm should use weighted means of line segments
5203  * rather than straight average values of points - tgl 97/01/21.
5204  */
5205 Datum
5207 {
5208  POLYGON *poly = PG_GETARG_POLYGON_P(0);
5209  CIRCLE *circle;
5210  int i;
5211 
5212  if (poly->npts < 1)
5213  ereport(ERROR,
5214  (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5215  errmsg("cannot convert empty polygon to circle")));
5216 
5217  circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5218 
5219  circle->center.x = 0;
5220  circle->center.y = 0;
5221  circle->radius = 0;
5222 
5223  for (i = 0; i < poly->npts; i++)
5224  {
5225  circle->center.x += poly->p[i].x;
5226  circle->center.y += poly->p[i].y;
5227  }
5228  circle->center.x /= poly->npts;
5229  circle->center.y /= poly->npts;
5230 
5231  for (i = 0; i < poly->npts; i++)
5232  circle->radius += point_dt(&poly->p[i], &circle->center);
5233  circle->radius /= poly->npts;
5234 
5235  PG_RETURN_CIRCLE_P(circle);
5236 }
5237 
5238 
5239 /***********************************************************************
5240  **
5241  ** Private routines for multiple types.
5242  **
5243  ***********************************************************************/
5244 
5245 /*
5246  * Test to see if the point is inside the polygon, returns 1/0, or 2 if
5247  * the point is on the polygon.
5248  * Code adapted but not copied from integer-based routines in WN: A
5249  * Server for the HTTP
5250  * version 1.15.1, file wn/image.c
5251  * http://hopf.math.northwestern.edu/index.html
5252  * Description of algorithm: http://www.linuxjournal.com/article/2197
5253  * http://www.linuxjournal.com/article/2029
5254  */
5255 
5256 #define POINT_ON_POLYGON INT_MAX
5257 
5258 static int
5259 point_inside(Point *p, int npts, Point *plist)
5260 {
5261  double x0,
5262  y0;
5263  double prev_x,
5264  prev_y;
5265  int i = 0;
5266  double x,
5267  y;
5268  int cross,
5269  total_cross = 0;
5270 
5271  if (npts <= 0)
5272  return 0;
5273 
5274  /* compute first polygon point relative to single point */
5275  x0 = plist[0].x - p->x;
5276  y0 = plist[0].y - p->y;
5277 
5278  prev_x = x0;
5279  prev_y = y0;
5280  /* loop over polygon points and aggregate total_cross */
5281  for (i = 1; i < npts; i++)
5282  {
5283  /* compute next polygon point relative to single point */
5284  x = plist[i].x - p->x;
5285  y = plist[i].y - p->y;
5286 
5287  /* compute previous to current point crossing */
5288  if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
5289  return 2;
5290  total_cross += cross;
5291 
5292  prev_x = x;
5293  prev_y = y;
5294  }
5295 
5296  /* now do the first point */
5297  if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
5298  return 2;
5299  total_cross += cross;
5300 
5301  if (total_cross != 0)
5302  return 1;
5303  return 0;
5304 }
5305 
5306 
5307 /* lseg_crossing()
5308  * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
5309  * Returns +/-1 if one point is on the positive X-axis.
5310  * Returns 0 if both points are on the positive X-axis, or there is no crossing.
5311  * Returns POINT_ON_POLYGON if the segment contains (0,0).
5312  * Wow, that is one confusing API, but it is used above, and when summed,
5313  * can tell is if a point is in a polygon.
5314  */
5315 
5316 static int
5317 lseg_crossing(double x, double y, double prev_x, double prev_y)
5318 {
5319  double z;
5320  int y_sign;
5321 
5322  if (FPzero(y))
5323  { /* y == 0, on X axis */
5324  if (FPzero(x)) /* (x,y) is (0,0)? */
5325  return POINT_ON_POLYGON;
5326  else if (FPgt(x, 0))
5327  { /* x > 0 */
5328  if (FPzero(prev_y)) /* y and prev_y are zero */
5329  /* prev_x > 0? */
5330  return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5331  return FPlt(prev_y, 0) ? 1 : -1;
5332  }
5333  else
5334  { /* x < 0, x not on positive X axis */
5335  if (FPzero(prev_y))
5336  /* prev_x < 0? */
5337  return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON;
5338  return 0;
5339  }
5340  }
5341  else
5342  { /* y != 0 */
5343  /* compute y crossing direction from previous point */
5344  y_sign = FPgt(y, 0) ? 1 : -1;
5345 
5346  if (FPzero(prev_y))
5347  /* previous point was on X axis, so new point is either off or on */
5348  return FPlt(prev_x, 0) ? 0 : y_sign;
5349  else if (FPgt(y_sign * prev_y, 0))
5350  /* both above or below X axis */
5351  return 0; /* same sign */
5352  else
5353  { /* y and prev_y cross X-axis */
5354  if (FPge(x, 0) && FPgt(prev_x, 0))
5355  /* both non-negative so cross positive X-axis */
5356  return 2 * y_sign;
5357  if (FPlt(x, 0) && FPle(prev_x, 0))
5358  /* both non-positive so do not cross positive X-axis */
5359  return 0;
5360 
5361  /* x and y cross axises, see URL above point_inside() */
5362  z = (x - prev_x) * y - (y - prev_y) * x;
5363  if (FPzero(z))
5364  return POINT_ON_POLYGON;
5365  return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign;
5366  }
5367  }
5368 }
5369 
5370 
5371 static bool
5372 plist_same(int npts, Point *p1, Point *p2)
5373 {
5374  int i,
5375  ii,
5376  j;
5377 
5378  /* find match for first point */
5379  for (i = 0; i < npts; i++)
5380  {
5381  if ((FPeq(p2[i].x, p1[0].x))
5382  && (FPeq(p2[i].y, p1[0].y)))
5383  {
5384 
5385  /* match found? then look forward through remaining points */
5386  for (ii = 1, j = i + 1; ii < npts; ii++, j++)
5387  {
5388  if (j >= npts)
5389  j = 0;
5390  if ((!FPeq(p2[j].x, p1[ii].x))
5391  || (!FPeq(p2[j].y, p1[ii].y)))
5392  {
5393 #ifdef GEODEBUG
5394  printf("plist_same- %d failed forward match with %d\n", j, ii);
5395 #endif
5396  break;
5397  }
5398  }
5399 #ifdef GEODEBUG
5400  printf("plist_same- ii = %d/%d after forward match\n", ii, npts);
5401 #endif
5402  if (ii == npts)
5403  return true;
5404 
5405  /* match not found forwards? then look backwards */
5406  for (ii = 1, j = i - 1; ii < npts; ii++, j--)
5407  {
5408  if (j < 0)
5409  j = (npts - 1);
5410  if ((!FPeq(p2[j].x, p1[ii].x))
5411  || (!FPeq(p2[j].y, p1[ii].y)))
5412  {
5413 #ifdef GEODEBUG
5414  printf("plist_same- %d failed reverse match with %d\n", j, ii);
5415 #endif
5416  break;
5417  }
5418  }
5419 #ifdef GEODEBUG
5420  printf("plist_same- ii = %d/%d after reverse match\n", ii, npts);
5421 #endif
5422  if (ii == npts)
5423  return true;
5424  }
5425  }
5426 
5427  return false;
5428 }
5429 
5430 
5431 /*-------------------------------------------------------------------------
5432  * Determine the hypotenuse.
5433  *
5434  * If required, x and y are swapped to make x the larger number. The
5435  * traditional formula of x^2+y^2 is rearranged to factor x outside the
5436  * sqrt. This allows computation of the hypotenuse for significantly
5437  * larger values, and with a higher precision than when using the naive
5438  * formula. In particular, this cannot overflow unless the final result
5439  * would be out-of-range.
5440  *
5441  * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
5442  * = x * sqrt( 1 + y^2/x^2 )
5443  * = x * sqrt( 1 + y/x * y/x )
5444  *
5445  * It is expected that this routine will eventually be replaced with the
5446  * C99 hypot() function.
5447  *
5448  * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
5449  * case of hypot(inf,nan) results in INF, and not NAN.
5450  *-----------------------------------------------------------------------
5451  */
5452 double
5453 pg_hypot(double x, double y)
5454 {
5455  double yx;
5456 
5457  /* Handle INF and NaN properly */
5458  if (isinf(x) || isinf(y))
5459  return get_float8_infinity();
5460 
5461  if (isnan(x) || isnan(y))
5462  return get_float8_nan();
5463 
5464  /* Else, drop any minus signs */
5465  x = fabs(x);
5466  y = fabs(y);
5467 
5468  /* Swap x and y if needed to make x the larger one */
5469  if (x < y)
5470  {
5471  double temp = x;
5472 
5473  x = y;
5474  y = temp;
5475  }
5476 
5477  /*
5478  * If y is zero, the hypotenuse is x. This test saves a few cycles in
5479  * such cases, but more importantly it also protects against
5480  * divide-by-zero errors, since now x >= y.
5481  */
5482  if (y == 0.0)
5483  return x;
5484 
5485  /* Determine the hypotenuse */
5486  yx = y / x;
5487  return x * sqrt(1.0 + (yx * yx));
5488 }
Datum boxes_bound_box(PG_FUNCTION_ARGS)
Definition: geo_ops.c:4205
Datum box_right(PG_FUNCTION_ARGS)
Definition: geo_ops.c:563