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