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