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-2021, 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 (isinf(m))
1059  {
1060  /* vertical - use "x = C" */
1061  result->A = -1.0;
1062  result->B = 0.0;
1063  result->C = pt->x;
1064  }
1065  else if (m == 0)
1066  {
1067  /* horizontal - use "y = C" */
1068  result->A = 0.0;
1069  result->B = -1.0;
1070  result->C = pt->y;
1071  }
1072  else
1073  {
1074  /* use "mx - y + yinter = 0" */
1075  result->A = m;
1076  result->B = -1.0;
1077  result->C = float8_mi(pt->y, float8_mul(m, pt->x));
1078  /* on some platforms, the preceding expression tends to produce -0 */
1079  if (result->C == 0.0)
1080  result->C = 0.0;
1081  }
1082 }
1083 
1084 /* line_construct_pp()
1085  * two points
1086  */
1087 Datum
1089 {
1090  Point *pt1 = PG_GETARG_POINT_P(0);
1091  Point *pt2 = PG_GETARG_POINT_P(1);
1092  LINE *result = (LINE *) palloc(sizeof(LINE));
1093 
1094  if (point_eq_point(pt1, pt2))
1095  ereport(ERROR,
1096  (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
1097  errmsg("invalid line specification: must be two distinct points")));
1098 
1099  line_construct(result, pt1, point_sl(pt1, pt2));
1100 
1101  PG_RETURN_LINE_P(result);
1102 }
1103 
1104 
1105 /*----------------------------------------------------------
1106  * Relative position routines.
1107  *---------------------------------------------------------*/
1108 
1109 Datum
1111 {
1112  LINE *l1 = PG_GETARG_LINE_P(0);
1113  LINE *l2 = PG_GETARG_LINE_P(1);
1114 
1115  PG_RETURN_BOOL(line_interpt_line(NULL, l1, l2));
1116 }
1117 
1118 Datum
1120 {
1121  LINE *l1 = PG_GETARG_LINE_P(0);
1122  LINE *l2 = PG_GETARG_LINE_P(1);
1123 
1124  PG_RETURN_BOOL(!line_interpt_line(NULL, l1, l2));
1125 }
1126 
1127 Datum
1129 {
1130  LINE *l1 = PG_GETARG_LINE_P(0);
1131  LINE *l2 = PG_GETARG_LINE_P(1);
1132 
1133  if (FPzero(l1->A))
1134  PG_RETURN_BOOL(FPzero(l2->B));
1135  if (FPzero(l2->A))
1136  PG_RETURN_BOOL(FPzero(l1->B));
1137  if (FPzero(l1->B))
1138  PG_RETURN_BOOL(FPzero(l2->A));
1139  if (FPzero(l2->B))
1140  PG_RETURN_BOOL(FPzero(l1->A));
1141 
1143  float8_mul(l1->B, l2->B)), -1.0));
1144 }
1145 
1146 Datum
1148 {
1149  LINE *line = PG_GETARG_LINE_P(0);
1150 
1151  PG_RETURN_BOOL(FPzero(line->B));
1152 }
1153 
1154 Datum
1156 {
1157  LINE *line = PG_GETARG_LINE_P(0);
1158 
1159  PG_RETURN_BOOL(FPzero(line->A));
1160 }
1161 
1162 
1163 /*
1164  * Check whether the two lines are the same
1165  */
1166 Datum
1168 {
1169  LINE *l1 = PG_GETARG_LINE_P(0);
1170  LINE *l2 = PG_GETARG_LINE_P(1);
1171  float8 ratio;
1172 
1173  /* If any NaNs are involved, insist on exact equality */
1174  if (unlikely(isnan(l1->A) || isnan(l1->B) || isnan(l1->C) ||
1175  isnan(l2->A) || isnan(l2->B) || isnan(l2->C)))
1176  {
1177  PG_RETURN_BOOL(float8_eq(l1->A, l2->A) &&
1178  float8_eq(l1->B, l2->B) &&
1179  float8_eq(l1->C, l2->C));
1180  }
1181 
1182  /* Otherwise, lines whose parameters are proportional are the same */
1183  if (!FPzero(l2->A))
1184  ratio = float8_div(l1->A, l2->A);
1185  else if (!FPzero(l2->B))
1186  ratio = float8_div(l1->B, l2->B);
1187  else if (!FPzero(l2->C))
1188  ratio = float8_div(l1->C, l2->C);
1189  else
1190  ratio = 1.0;
1191 
1192  PG_RETURN_BOOL(FPeq(l1->A, float8_mul(ratio, l2->A)) &&
1193  FPeq(l1->B, float8_mul(ratio, l2->B)) &&
1194  FPeq(l1->C, float8_mul(ratio, l2->C)));
1195 }
1196 
1197 
1198 /*----------------------------------------------------------
1199  * Line arithmetic routines.
1200  *---------------------------------------------------------*/
1201 
1202 /*
1203  * Return slope of the line
1204  */
1205 static inline float8
1207 {
1208  if (FPzero(line->A))
1209  return 0.0;
1210  if (FPzero(line->B))
1211  return get_float8_infinity();
1212  return float8_div(line->A, -line->B);
1213 }
1214 
1215 
1216 /*
1217  * Return inverse slope of the line
1218  */
1219 static inline float8
1221 {
1222  if (FPzero(line->A))
1223  return get_float8_infinity();
1224  if (FPzero(line->B))
1225  return 0.0;
1226  return float8_div(line->B, line->A);
1227 }
1228 
1229 
1230 /* line_distance()
1231  * Distance between two lines.
1232  */
1233 Datum
1235 {
1236  LINE *l1 = PG_GETARG_LINE_P(0);
1237  LINE *l2 = PG_GETARG_LINE_P(1);
1238  float8 ratio;
1239 
1240  if (line_interpt_line(NULL, l1, l2)) /* intersecting? */
1241  PG_RETURN_FLOAT8(0.0);
1242 
1243  if (!FPzero(l1->A) && !isnan(l1->A) && !FPzero(l2->A) && !isnan(l2->A))
1244  ratio = float8_div(l1->A, l2->A);
1245  else if (!FPzero(l1->B) && !isnan(l1->B) && !FPzero(l2->B) && !isnan(l2->B))
1246  ratio = float8_div(l1->B, l2->B);
1247  else
1248  ratio = 1.0;
1249 
1251  float8_mul(ratio, l2->C))),
1252  HYPOT(l1->A, l1->B)));
1253 }
1254 
1255 /* line_interpt()
1256  * Point where two lines l1, l2 intersect (if any)
1257  */
1258 Datum
1260 {
1261  LINE *l1 = PG_GETARG_LINE_P(0);
1262  LINE *l2 = PG_GETARG_LINE_P(1);
1263  Point *result;
1264 
1265  result = (Point *) palloc(sizeof(Point));
1266 
1267  if (!line_interpt_line(result, l1, l2))
1268  PG_RETURN_NULL();
1269  PG_RETURN_POINT_P(result);
1270 }
1271 
1272 /*
1273  * Internal version of line_interpt
1274  *
1275  * Return whether two lines intersect. If *result is not NULL, it is set to
1276  * the intersection point.
1277  *
1278  * NOTE: If the lines are identical then we will find they are parallel
1279  * and report "no intersection". This is a little weird, but since
1280  * there's no *unique* intersection, maybe it's appropriate behavior.
1281  *
1282  * If the lines have NaN constants, we will return true, and the intersection
1283  * point would have NaN coordinates. We shouldn't return false in this case
1284  * because that would mean the lines are parallel.
1285  */
1286 static bool
1287 line_interpt_line(Point *result, LINE *l1, LINE *l2)
1288 {
1289  float8 x,
1290  y;
1291 
1292  if (!FPzero(l1->B))
1293  {
1294  if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
1295  return false;
1296 
1297  x = float8_div(float8_mi(float8_mul(l1->B, l2->C),
1298  float8_mul(l2->B, l1->C)),
1299  float8_mi(float8_mul(l1->A, l2->B),
1300  float8_mul(l2->A, l1->B)));
1301  y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B);
1302  }
1303  else if (!FPzero(l2->B))
1304  {
1305  if (FPeq(l1->A, float8_mul(l2->A, float8_div(l1->B, l2->B))))
1306  return false;
1307 
1308  x = float8_div(float8_mi(float8_mul(l2->B, l1->C),
1309  float8_mul(l1->B, l2->C)),
1310  float8_mi(float8_mul(l2->A, l1->B),
1311  float8_mul(l1->A, l2->B)));
1312  y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B);
1313  }
1314  else
1315  return false;
1316 
1317  /* On some platforms, the preceding expressions tend to produce -0. */
1318  if (x == 0.0)
1319  x = 0.0;
1320  if (y == 0.0)
1321  y = 0.0;
1322 
1323  if (result != NULL)
1324  point_construct(result, x, y);
1325 
1326  return true;
1327 }
1328 
1329 
1330 /***********************************************************************
1331  **
1332  ** Routines for 2D paths (sequences of line segments, also
1333  ** called `polylines').
1334  **
1335  ** This is not a general package for geometric paths,
1336  ** which of course include polygons; the emphasis here
1337  ** is on (for example) usefulness in wire layout.
1338  **
1339  ***********************************************************************/
1340 
1341 /*----------------------------------------------------------
1342  * String to path / path to string conversion.
1343  * External format:
1344  * "((xcoord, ycoord),... )"
1345  * "[(xcoord, ycoord),... ]"
1346  * "(xcoord, ycoord),... "
1347  * "[xcoord, ycoord,... ]"
1348  * Also support older format:
1349  * "(closed, npts, xcoord, ycoord,... )"
1350  *---------------------------------------------------------*/
1351 
1352 Datum
1354 {
1355  PATH *path = PG_GETARG_PATH_P(0);
1356  float8 area = 0.0;
1357  int i,
1358  j;
1359 
1360  if (!path->closed)
1361  PG_RETURN_NULL();
1362 
1363  for (i = 0; i < path->npts; i++)
1364  {
1365  j = (i + 1) % path->npts;
1366  area = float8_pl(area, float8_mul(path->p[i].x, path->p[j].y));
1367  area = float8_mi(area, float8_mul(path->p[i].y, path->p[j].x));
1368  }
1369 
1370  PG_RETURN_FLOAT8(float8_div(fabs(area), 2.0));
1371 }
1372 
1373 
1374 Datum
1376 {
1377  char *str = PG_GETARG_CSTRING(0);
1378  PATH *path;
1379  bool isopen;
1380  char *s;
1381  int npts;
1382  int size;
1383  int base_size;
1384  int depth = 0;
1385 
1386  if ((npts = pair_count(str, ',')) <= 0)
1387  ereport(ERROR,
1388  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1389  errmsg("invalid input syntax for type %s: \"%s\"",
1390  "path", str)));
1391 
1392  s = str;
1393  while (isspace((unsigned char) *s))
1394  s++;
1395 
1396  /* skip single leading paren */
1397  if ((*s == LDELIM) && (strrchr(s, LDELIM) == s))
1398  {
1399  s++;
1400  depth++;
1401  }
1402 
1403  base_size = sizeof(path->p[0]) * npts;
1404  size = offsetof(PATH, p) + base_size;
1405 
1406  /* Check for integer overflow */
1407  if (base_size / npts != sizeof(path->p[0]) || size <= base_size)
1408  ereport(ERROR,
1409  (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
1410  errmsg("too many points requested")));
1411 
1412  path = (PATH *) palloc(size);
1413 
1414  SET_VARSIZE(path, size);
1415  path->npts = npts;
1416 
1417  path_decode(s, true, npts, &(path->p[0]), &isopen, &s, "path", str);
1418 
1419  if (depth >= 1)
1420  {
1421  if (*s++ != RDELIM)
1422  ereport(ERROR,
1423  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1424  errmsg("invalid input syntax for type %s: \"%s\"",
1425  "path", str)));
1426  while (isspace((unsigned char) *s))
1427  s++;
1428  }
1429  if (*s != '\0')
1430  ereport(ERROR,
1431  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
1432  errmsg("invalid input syntax for type %s: \"%s\"",
1433  "path", str)));
1434 
1435  path->closed = (!isopen);
1436  /* prevent instability in unused pad bytes */
1437  path->dummy = 0;
1438 
1439  PG_RETURN_PATH_P(path);
1440 }
1441 
1442 
1443 Datum
1445 {
1446  PATH *path = PG_GETARG_PATH_P(0);
1447 
1448  PG_RETURN_CSTRING(path_encode(path->closed ? PATH_CLOSED : PATH_OPEN, path->npts, path->p));
1449 }
1450 
1451 /*
1452  * path_recv - converts external binary format to path
1453  *
1454  * External representation is closed flag (a boolean byte), int32 number
1455  * of points, and the points.
1456  */
1457 Datum
1459 {
1461  PATH *path;
1462  int closed;
1463  int32 npts;
1464  int32 i;
1465  int size;
1466 
1467  closed = pq_getmsgbyte(buf);
1468  npts = pq_getmsgint(buf, sizeof(int32));
1469  if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(PATH, p)) / sizeof(Point)))
1470  ereport(ERROR,
1471  (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
1472  errmsg("invalid number of points in external \"path\" value")));
1473 
1474  size = offsetof(PATH, p) + sizeof(path->p[0]) * npts;
1475  path = (PATH *) palloc(size);
1476 
1477  SET_VARSIZE(path, size);
1478  path->npts = npts;
1479  path->closed = (closed ? 1 : 0);
1480  /* prevent instability in unused pad bytes */
1481  path->dummy = 0;
1482 
1483  for (i = 0; i < npts; i++)
1484  {
1485  path->p[i].x = pq_getmsgfloat8(buf);
1486  path->p[i].y = pq_getmsgfloat8(buf);
1487  }
1488 
1489  PG_RETURN_PATH_P(path);
1490 }
1491 
1492 /*
1493  * path_send - converts path to binary format
1494  */
1495 Datum
1497 {
1498  PATH *path = PG_GETARG_PATH_P(0);
1500  int32 i;
1501 
1502  pq_begintypsend(&buf);
1503  pq_sendbyte(&buf, path->closed ? 1 : 0);
1504  pq_sendint32(&buf, path->npts);
1505  for (i = 0; i < path->npts; i++)
1506  {
1507  pq_sendfloat8(&buf, path->p[i].x);
1508  pq_sendfloat8(&buf, path->p[i].y);
1509  }
1511 }
1512 
1513 
1514 /*----------------------------------------------------------
1515  * Relational operators.
1516  * These are based on the path cardinality,
1517  * as stupid as that sounds.
1518  *
1519  * Better relops and access methods coming soon.
1520  *---------------------------------------------------------*/
1521 
1522 Datum
1524 {
1525  PATH *p1 = PG_GETARG_PATH_P(0);
1526  PATH *p2 = PG_GETARG_PATH_P(1);
1527 
1528  PG_RETURN_BOOL(p1->npts < p2->npts);
1529 }
1530 
1531 Datum
1533 {
1534  PATH *p1 = PG_GETARG_PATH_P(0);
1535  PATH *p2 = PG_GETARG_PATH_P(1);
1536 
1537  PG_RETURN_BOOL(p1->npts > p2->npts);
1538 }
1539 
1540 Datum
1542 {
1543  PATH *p1 = PG_GETARG_PATH_P(0);
1544  PATH *p2 = PG_GETARG_PATH_P(1);
1545 
1546  PG_RETURN_BOOL(p1->npts == p2->npts);
1547 }
1548 
1549 Datum
1551 {
1552  PATH *p1 = PG_GETARG_PATH_P(0);
1553  PATH *p2 = PG_GETARG_PATH_P(1);
1554 
1555  PG_RETURN_BOOL(p1->npts <= p2->npts);
1556 }
1557 
1558 Datum
1560 {
1561  PATH *p1 = PG_GETARG_PATH_P(0);
1562  PATH *p2 = PG_GETARG_PATH_P(1);
1563 
1564  PG_RETURN_BOOL(p1->npts >= p2->npts);
1565 }
1566 
1567 /*----------------------------------------------------------
1568  * Conversion operators.
1569  *---------------------------------------------------------*/
1570 
1571 Datum
1573 {
1574  PATH *path = PG_GETARG_PATH_P(0);
1575 
1576  PG_RETURN_BOOL(path->closed);
1577 }
1578 
1579 Datum
1581 {
1582  PATH *path = PG_GETARG_PATH_P(0);
1583 
1584  PG_RETURN_BOOL(!path->closed);
1585 }
1586 
1587 Datum
1589 {
1590  PATH *path = PG_GETARG_PATH_P(0);
1591 
1592  PG_RETURN_INT32(path->npts);
1593 }
1594 
1595 
1596 Datum
1598 {
1599  PATH *path = PG_GETARG_PATH_P_COPY(0);
1600 
1601  path->closed = true;
1602 
1603  PG_RETURN_PATH_P(path);
1604 }
1605 
1606 Datum
1608 {
1609  PATH *path = PG_GETARG_PATH_P_COPY(0);
1610 
1611  path->closed = false;
1612 
1613  PG_RETURN_PATH_P(path);
1614 }
1615 
1616 
1617 /* path_inter -
1618  * Does p1 intersect p2 at any point?
1619  * Use bounding boxes for a quick (O(n)) check, then do a
1620  * O(n^2) iterative edge check.
1621  */
1622 Datum
1624 {
1625  PATH *p1 = PG_GETARG_PATH_P(0);
1626  PATH *p2 = PG_GETARG_PATH_P(1);
1627  BOX b1,
1628  b2;
1629  int i,
1630  j;
1631  LSEG seg1,
1632  seg2;
1633 
1634  Assert(p1->npts > 0 && p2->npts > 0);
1635 
1636  b1.high.x = b1.low.x = p1->p[0].x;
1637  b1.high.y = b1.low.y = p1->p[0].y;
1638  for (i = 1; i < p1->npts; i++)
1639  {
1640  b1.high.x = float8_max(p1->p[i].x, b1.high.x);
1641  b1.high.y = float8_max(p1->p[i].y, b1.high.y);
1642  b1.low.x = float8_min(p1->p[i].x, b1.low.x);
1643  b1.low.y = float8_min(p1->p[i].y, b1.low.y);
1644  }
1645  b2.high.x = b2.low.x = p2->p[0].x;
1646  b2.high.y = b2.low.y = p2->p[0].y;
1647  for (i = 1; i < p2->npts; i++)
1648  {
1649  b2.high.x = float8_max(p2->p[i].x, b2.high.x);
1650  b2.high.y = float8_max(p2->p[i].y, b2.high.y);
1651  b2.low.x = float8_min(p2->p[i].x, b2.low.x);
1652  b2.low.y = float8_min(p2->p[i].y, b2.low.y);
1653  }
1654  if (!box_ov(&b1, &b2))
1655  PG_RETURN_BOOL(false);
1656 
1657  /* pairwise check lseg intersections */
1658  for (i = 0; i < p1->npts; i++)
1659  {
1660  int iprev;
1661 
1662  if (i > 0)
1663  iprev = i - 1;
1664  else
1665  {
1666  if (!p1->closed)
1667  continue;
1668  iprev = p1->npts - 1; /* include the closure segment */
1669  }
1670 
1671  for (j = 0; j < p2->npts; j++)
1672  {
1673  int jprev;
1674 
1675  if (j > 0)
1676  jprev = j - 1;
1677  else
1678  {
1679  if (!p2->closed)
1680  continue;
1681  jprev = p2->npts - 1; /* include the closure segment */
1682  }
1683 
1684  statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1685  statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1686  if (lseg_interpt_lseg(NULL, &seg1, &seg2))
1687  PG_RETURN_BOOL(true);
1688  }
1689  }
1690 
1691  /* if we dropped through, no two segs intersected */
1692  PG_RETURN_BOOL(false);
1693 }
1694 
1695 /* path_distance()
1696  * This essentially does a cartesian product of the lsegs in the
1697  * two paths, and finds the min distance between any two lsegs
1698  */
1699 Datum
1701 {
1702  PATH *p1 = PG_GETARG_PATH_P(0);
1703  PATH *p2 = PG_GETARG_PATH_P(1);
1704  float8 min = 0.0; /* initialize to keep compiler quiet */
1705  bool have_min = false;
1706  float8 tmp;
1707  int i,
1708  j;
1709  LSEG seg1,
1710  seg2;
1711 
1712  for (i = 0; i < p1->npts; i++)
1713  {
1714  int iprev;
1715 
1716  if (i > 0)
1717  iprev = i - 1;
1718  else
1719  {
1720  if (!p1->closed)
1721  continue;
1722  iprev = p1->npts - 1; /* include the closure segment */
1723  }
1724 
1725  for (j = 0; j < p2->npts; j++)
1726  {
1727  int jprev;
1728 
1729  if (j > 0)
1730  jprev = j - 1;
1731  else
1732  {
1733  if (!p2->closed)
1734  continue;
1735  jprev = p2->npts - 1; /* include the closure segment */
1736  }
1737 
1738  statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]);
1739  statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]);
1740 
1741  tmp = lseg_closept_lseg(NULL, &seg1, &seg2);
1742  if (!have_min || float8_lt(tmp, min))
1743  {
1744  min = tmp;
1745  have_min = true;
1746  }
1747  }
1748  }
1749 
1750  if (!have_min)
1751  PG_RETURN_NULL();
1752 
1753  PG_RETURN_FLOAT8(min);
1754 }
1755 
1756 
1757 /*----------------------------------------------------------
1758  * "Arithmetic" operations.
1759  *---------------------------------------------------------*/
1760 
1761 Datum
1763 {
1764  PATH *path = PG_GETARG_PATH_P(0);
1765  float8 result = 0.0;
1766  int i;
1767 
1768  for (i = 0; i < path->npts; i++)
1769  {
1770  int iprev;
1771 
1772  if (i > 0)
1773  iprev = i - 1;
1774  else
1775  {
1776  if (!path->closed)
1777  continue;
1778  iprev = path->npts - 1; /* include the closure segment */
1779  }
1780 
1781  result = float8_pl(result, point_dt(&path->p[iprev], &path->p[i]));
1782  }
1783 
1784  PG_RETURN_FLOAT8(result);
1785 }
1786 
1787 /***********************************************************************
1788  **
1789  ** Routines for 2D points.
1790  **
1791  ***********************************************************************/
1792 
1793 /*----------------------------------------------------------
1794  * String to point, point to string conversion.
1795  * External format:
1796  * "(x,y)"
1797  * "x,y"
1798  *---------------------------------------------------------*/
1799 
1800 Datum
1802 {
1803  char *str = PG_GETARG_CSTRING(0);
1804  Point *point = (Point *) palloc(sizeof(Point));
1805 
1806  pair_decode(str, &point->x, &point->y, NULL, "point", str);
1807  PG_RETURN_POINT_P(point);
1808 }
1809 
1810 Datum
1812 {
1813  Point *pt = PG_GETARG_POINT_P(0);
1814 
1816 }
1817 
1818 /*
1819  * point_recv - converts external binary format to point
1820  */
1821 Datum
1823 {
1825  Point *point;
1826 
1827  point = (Point *) palloc(sizeof(Point));
1828  point->x = pq_getmsgfloat8(buf);
1829  point->y = pq_getmsgfloat8(buf);
1830  PG_RETURN_POINT_P(point);
1831 }
1832 
1833 /*
1834  * point_send - converts point to binary format
1835  */
1836 Datum
1838 {
1839  Point *pt = PG_GETARG_POINT_P(0);
1841 
1842  pq_begintypsend(&buf);
1843  pq_sendfloat8(&buf, pt->x);
1844  pq_sendfloat8(&buf, pt->y);
1846 }
1847 
1848 
1849 /*
1850  * Initialize a point
1851  */
1852 static inline void
1854 {
1855  result->x = x;
1856  result->y = y;
1857 }
1858 
1859 
1860 /*----------------------------------------------------------
1861  * Relational operators for Points.
1862  * Since we do have a sense of coordinates being
1863  * "equal" to a given accuracy (point_vert, point_horiz),
1864  * the other ops must preserve that sense. This means
1865  * that results may, strictly speaking, be a lie (unless
1866  * EPSILON = 0.0).
1867  *---------------------------------------------------------*/
1868 
1869 Datum
1871 {
1872  Point *pt1 = PG_GETARG_POINT_P(0);
1873  Point *pt2 = PG_GETARG_POINT_P(1);
1874 
1875  PG_RETURN_BOOL(FPlt(pt1->x, pt2->x));
1876 }
1877 
1878 Datum
1880 {
1881  Point *pt1 = PG_GETARG_POINT_P(0);
1882  Point *pt2 = PG_GETARG_POINT_P(1);
1883 
1884  PG_RETURN_BOOL(FPgt(pt1->x, pt2->x));
1885 }
1886 
1887 Datum
1889 {
1890  Point *pt1 = PG_GETARG_POINT_P(0);
1891  Point *pt2 = PG_GETARG_POINT_P(1);
1892 
1893  PG_RETURN_BOOL(FPgt(pt1->y, pt2->y));
1894 }
1895 
1896 Datum
1898 {
1899  Point *pt1 = PG_GETARG_POINT_P(0);
1900  Point *pt2 = PG_GETARG_POINT_P(1);
1901 
1902  PG_RETURN_BOOL(FPlt(pt1->y, pt2->y));
1903 }
1904 
1905 Datum
1907 {
1908  Point *pt1 = PG_GETARG_POINT_P(0);
1909  Point *pt2 = PG_GETARG_POINT_P(1);
1910 
1911  PG_RETURN_BOOL(FPeq(pt1->x, pt2->x));
1912 }
1913 
1914 Datum
1916 {
1917  Point *pt1 = PG_GETARG_POINT_P(0);
1918  Point *pt2 = PG_GETARG_POINT_P(1);
1919 
1920  PG_RETURN_BOOL(FPeq(pt1->y, pt2->y));
1921 }
1922 
1923 Datum
1925 {
1926  Point *pt1 = PG_GETARG_POINT_P(0);
1927  Point *pt2 = PG_GETARG_POINT_P(1);
1928 
1929  PG_RETURN_BOOL(point_eq_point(pt1, pt2));
1930 }
1931 
1932 Datum
1934 {
1935  Point *pt1 = PG_GETARG_POINT_P(0);
1936  Point *pt2 = PG_GETARG_POINT_P(1);
1937 
1938  PG_RETURN_BOOL(!point_eq_point(pt1, pt2));
1939 }
1940 
1941 
1942 /*
1943  * Check whether the two points are the same
1944  */
1945 static inline bool
1947 {
1948  /* If any NaNs are involved, insist on exact equality */
1949  if (unlikely(isnan(pt1->x) || isnan(pt1->y) ||
1950  isnan(pt2->x) || isnan(pt2->y)))
1951  return (float8_eq(pt1->x, pt2->x) && float8_eq(pt1->y, pt2->y));
1952 
1953  return (FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y));
1954 }
1955 
1956 
1957 /*----------------------------------------------------------
1958  * "Arithmetic" operators on points.
1959  *---------------------------------------------------------*/
1960 
1961 Datum
1963 {
1964  Point *pt1 = PG_GETARG_POINT_P(0);
1965  Point *pt2 = PG_GETARG_POINT_P(1);
1966 
1967  PG_RETURN_FLOAT8(point_dt(pt1, pt2));
1968 }
1969 
1970 static inline float8
1971 point_dt(Point *pt1, Point *pt2)
1972 {
1973  return HYPOT(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y));
1974 }
1975 
1976 Datum
1978 {
1979  Point *pt1 = PG_GETARG_POINT_P(0);
1980  Point *pt2 = PG_GETARG_POINT_P(1);
1981 
1982  PG_RETURN_FLOAT8(point_sl(pt1, pt2));
1983 }
1984 
1985 
1986 /*
1987  * Return slope of two points
1988  *
1989  * Note that this function returns Inf when the points are the same.
1990  */
1991 static inline float8
1992 point_sl(Point *pt1, Point *pt2)
1993 {
1994  if (FPeq(pt1->x, pt2->x))
1995  return get_float8_infinity();
1996  if (FPeq(pt1->y, pt2->y))
1997  return 0.0;
1998  return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x));
1999 }
2000 
2001 
2002 /*
2003  * Return inverse slope of two points
2004  *
2005  * Note that this function returns 0.0 when the points are the same.
2006  */
2007 static inline float8
2009 {
2010  if (FPeq(pt1->x, pt2->x))
2011  return 0.0;
2012  if (FPeq(pt1->y, pt2->y))
2013  return get_float8_infinity();
2014  return float8_div(float8_mi(pt1->x, pt2->x), float8_mi(pt2->y, pt1->y));
2015 }
2016 
2017 
2018 /***********************************************************************
2019  **
2020  ** Routines for 2D line segments.
2021  **
2022  ***********************************************************************/
2023 
2024 /*----------------------------------------------------------
2025  * String to lseg, lseg to string conversion.
2026  * External forms: "[(x1, y1), (x2, y2)]"
2027  * "(x1, y1), (x2, y2)"
2028  * "x1, y1, x2, y2"
2029  * closed form ok "((x1, y1), (x2, y2))"
2030  * (old form) "(x1, y1, x2, y2)"
2031  *---------------------------------------------------------*/
2032 
2033 Datum
2035 {
2036  char *str = PG_GETARG_CSTRING(0);
2037  LSEG *lseg = (LSEG *) palloc(sizeof(LSEG));
2038  bool isopen;
2039 
2040  path_decode(str, true, 2, &lseg->p[0], &isopen, NULL, "lseg", str);
2041  PG_RETURN_LSEG_P(lseg);
2042 }
2043 
2044 
2045 Datum
2047 {
2048  LSEG *ls = PG_GETARG_LSEG_P(0);
2049 
2050  PG_RETURN_CSTRING(path_encode(PATH_OPEN, 2, &ls->p[0]));
2051 }
2052 
2053 /*
2054  * lseg_recv - converts external binary format to lseg
2055  */
2056 Datum
2058 {
2060  LSEG *lseg;
2061 
2062  lseg = (LSEG *) palloc(sizeof(LSEG));
2063 
2064  lseg->p[0].x = pq_getmsgfloat8(buf);
2065  lseg->p[0].y = pq_getmsgfloat8(buf);
2066  lseg->p[1].x = pq_getmsgfloat8(buf);
2067  lseg->p[1].y = pq_getmsgfloat8(buf);
2068 
2069  PG_RETURN_LSEG_P(lseg);
2070 }
2071 
2072 /*
2073  * lseg_send - converts lseg to binary format
2074  */
2075 Datum
2077 {
2078  LSEG *ls = PG_GETARG_LSEG_P(0);
2080 
2081  pq_begintypsend(&buf);
2082  pq_sendfloat8(&buf, ls->p[0].x);
2083  pq_sendfloat8(&buf, ls->p[0].y);
2084  pq_sendfloat8(&buf, ls->p[1].x);
2085  pq_sendfloat8(&buf, ls->p[1].y);
2087 }
2088 
2089 
2090 /* lseg_construct -
2091  * form a LSEG from two Points.
2092  */
2093 Datum
2095 {
2096  Point *pt1 = PG_GETARG_POINT_P(0);
2097  Point *pt2 = PG_GETARG_POINT_P(1);
2098  LSEG *result = (LSEG *) palloc(sizeof(LSEG));
2099 
2100  statlseg_construct(result, pt1, pt2);
2101 
2102  PG_RETURN_LSEG_P(result);
2103 }
2104 
2105 /* like lseg_construct, but assume space already allocated */
2106 static inline void
2108 {
2109  lseg->p[0].x = pt1->x;
2110  lseg->p[0].y = pt1->y;
2111  lseg->p[1].x = pt2->x;
2112  lseg->p[1].y = pt2->y;
2113 }
2114 
2115 
2116 /*
2117  * Return slope of the line segment
2118  */
2119 static inline float8
2121 {
2122  return point_sl(&lseg->p[0], &lseg->p[1]);
2123 }
2124 
2125 
2126 /*
2127  * Return inverse slope of the line segment
2128  */
2129 static inline float8
2131 {
2132  return point_invsl(&lseg->p[0], &lseg->p[1]);
2133 }
2134 
2135 
2136 Datum
2138 {
2139  LSEG *lseg = PG_GETARG_LSEG_P(0);
2140 
2141  PG_RETURN_FLOAT8(point_dt(&lseg->p[0], &lseg->p[1]));
2142 }
2143 
2144 /*----------------------------------------------------------
2145  * Relative position routines.
2146  *---------------------------------------------------------*/
2147 
2148 /*
2149  ** find intersection of the two lines, and see if it falls on
2150  ** both segments.
2151  */
2152 Datum
2154 {
2155  LSEG *l1 = PG_GETARG_LSEG_P(0);
2156  LSEG *l2 = PG_GETARG_LSEG_P(1);
2157 
2158  PG_RETURN_BOOL(lseg_interpt_lseg(NULL, l1, l2));
2159 }
2160 
2161 
2162 Datum
2164 {
2165  LSEG *l1 = PG_GETARG_LSEG_P(0);
2166  LSEG *l2 = PG_GETARG_LSEG_P(1);
2167 
2168  PG_RETURN_BOOL(FPeq(lseg_sl(l1), lseg_sl(l2)));
2169 }
2170 
2171 /*
2172  * Determine if two line segments are perpendicular.
2173  */
2174 Datum
2176 {
2177  LSEG *l1 = PG_GETARG_LSEG_P(0);
2178  LSEG *l2 = PG_GETARG_LSEG_P(1);
2179 
2181 }
2182 
2183 Datum
2185 {
2186  LSEG *lseg = PG_GETARG_LSEG_P(0);
2187 
2188  PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
2189 }
2190 
2191 Datum
2193 {
2194  LSEG *lseg = PG_GETARG_LSEG_P(0);
2195 
2196  PG_RETURN_BOOL(FPeq(lseg->p[0].y, lseg->p[1].y));
2197 }
2198 
2199 
2200 Datum
2202 {
2203  LSEG *l1 = PG_GETARG_LSEG_P(0);
2204  LSEG *l2 = PG_GETARG_LSEG_P(1);
2205 
2206  PG_RETURN_BOOL(point_eq_point(&l1->p[0], &l2->p[0]) &&
2207  point_eq_point(&l1->p[1], &l2->p[1]));
2208 }
2209 
2210 Datum
2212 {
2213  LSEG *l1 = PG_GETARG_LSEG_P(0);
2214  LSEG *l2 = PG_GETARG_LSEG_P(1);
2215 
2216  PG_RETURN_BOOL(!point_eq_point(&l1->p[0], &l2->p[0]) ||
2217  !point_eq_point(&l1->p[1], &l2->p[1]));
2218 }
2219 
2220 Datum
2222 {
2223  LSEG *l1 = PG_GETARG_LSEG_P(0);
2224  LSEG *l2 = PG_GETARG_LSEG_P(1);
2225 
2226  PG_RETURN_BOOL(FPlt(point_dt(&l1->p[0], &l1->p[1]),
2227  point_dt(&l2->p[0], &l2->p[1])));
2228 }
2229 
2230 Datum
2232 {
2233  LSEG *l1 = PG_GETARG_LSEG_P(0);
2234  LSEG *l2 = PG_GETARG_LSEG_P(1);
2235 
2236  PG_RETURN_BOOL(FPle(point_dt(&l1->p[0], &l1->p[1]),
2237  point_dt(&l2->p[0], &l2->p[1])));
2238 }
2239 
2240 Datum
2242 {
2243  LSEG *l1 = PG_GETARG_LSEG_P(0);
2244  LSEG *l2 = PG_GETARG_LSEG_P(1);
2245 
2246  PG_RETURN_BOOL(FPgt(point_dt(&l1->p[0], &l1->p[1]),
2247  point_dt(&l2->p[0], &l2->p[1])));
2248 }
2249 
2250 Datum
2252 {
2253  LSEG *l1 = PG_GETARG_LSEG_P(0);
2254  LSEG *l2 = PG_GETARG_LSEG_P(1);
2255 
2256  PG_RETURN_BOOL(FPge(point_dt(&l1->p[0], &l1->p[1]),
2257  point_dt(&l2->p[0], &l2->p[1])));
2258 }
2259 
2260 
2261 /*----------------------------------------------------------
2262  * Line arithmetic routines.
2263  *---------------------------------------------------------*/
2264 
2265 /* lseg_distance -
2266  * If two segments don't intersect, then the closest
2267  * point will be from one of the endpoints to the other
2268  * segment.
2269  */
2270 Datum
2272 {
2273  LSEG *l1 = PG_GETARG_LSEG_P(0);
2274  LSEG *l2 = PG_GETARG_LSEG_P(1);
2275 
2276  PG_RETURN_FLOAT8(lseg_closept_lseg(NULL, l1, l2));
2277 }
2278 
2279 
2280 Datum
2282 {
2283  LSEG *lseg = PG_GETARG_LSEG_P(0);
2284  Point *result;
2285 
2286  result = (Point *) palloc(sizeof(Point));
2287 
2288  result->x = float8_div(float8_pl(lseg->p[0].x, lseg->p[1].x), 2.0);
2289  result->y = float8_div(float8_pl(lseg->p[0].y, lseg->p[1].y), 2.0);
2290 
2291  PG_RETURN_POINT_P(result);
2292 }
2293 
2294 
2295 /*
2296  * Return whether the two segments intersect. If *result is not NULL,
2297  * it is set to the intersection point.
2298  *
2299  * This function is almost perfectly symmetric, even though it doesn't look
2300  * like it. See lseg_interpt_line() for the other half of it.
2301  */
2302 static bool
2303 lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2)
2304 {
2305  Point interpt;
2306  LINE tmp;
2307 
2308  line_construct(&tmp, &l2->p[0], lseg_sl(l2));
2309  if (!lseg_interpt_line(&interpt, l1, &tmp))
2310  return false;
2311 
2312  /*
2313  * If the line intersection point isn't within l2, there is no valid
2314  * segment intersection point at all.
2315  */
2316  if (!lseg_contain_point(l2, &interpt))
2317  return false;
2318 
2319  if (result != NULL)
2320  *result = interpt;
2321 
2322  return true;
2323 }
2324 
2325 Datum
2327 {
2328  LSEG *l1 = PG_GETARG_LSEG_P(0);
2329  LSEG *l2 = PG_GETARG_LSEG_P(1);
2330  Point *result;
2331 
2332  result = (Point *) palloc(sizeof(Point));
2333 
2334  if (!lseg_interpt_lseg(result, l1, l2))
2335  PG_RETURN_NULL();
2336  PG_RETURN_POINT_P(result);
2337 }
2338 
2339 /***********************************************************************
2340  **
2341  ** Routines for position comparisons of differently-typed
2342  ** 2D objects.
2343  **
2344  ***********************************************************************/
2345 
2346 /*---------------------------------------------------------------------
2347  * dist_
2348  * Minimum distance from one object to another.
2349  *-------------------------------------------------------------------*/
2350 
2351 /*
2352  * Distance from a point to a line
2353  */
2354 Datum
2356 {
2357  Point *pt = PG_GETARG_POINT_P(0);
2358  LINE *line = PG_GETARG_LINE_P(1);
2359 
2360  PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt));
2361 }
2362 
2363 /*
2364  * Distance from a line to a point
2365  */
2366 Datum
2368 {
2369  LINE *line = PG_GETARG_LINE_P(0);
2370  Point *pt = PG_GETARG_POINT_P(1);
2371 
2372  PG_RETURN_FLOAT8(line_closept_point(NULL, line, pt));
2373 }
2374 
2375 /*
2376  * Distance from a point to a lseg
2377  */
2378 Datum
2380 {
2381  Point *pt = PG_GETARG_POINT_P(0);
2382  LSEG *lseg = PG_GETARG_LSEG_P(1);
2383 
2384  PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt));
2385 }
2386 
2387 /*
2388  * Distance from a lseg to a point
2389  */
2390 Datum
2392 {
2393  LSEG *lseg = PG_GETARG_LSEG_P(0);
2394  Point *pt = PG_GETARG_POINT_P(1);
2395 
2396  PG_RETURN_FLOAT8(lseg_closept_point(NULL, lseg, pt));
2397 }
2398 
2399 static float8
2401 {
2402  float8 result = 0.0; /* keep compiler quiet */
2403  bool have_min = false;
2404  float8 tmp;
2405  int i;
2406  LSEG lseg;
2407 
2408  Assert(path->npts > 0);
2409 
2410  /*
2411  * The distance from a point to a path is the smallest distance from the
2412  * point to any of its constituent segments.
2413  */
2414  for (i = 0; i < path->npts; i++)
2415  {
2416  int iprev;
2417 
2418  if (i > 0)
2419  iprev = i - 1;
2420  else
2421  {
2422  if (!path->closed)
2423  continue;
2424  iprev = path->npts - 1; /* Include the closure segment */
2425  }
2426 
2427  statlseg_construct(&lseg, &path->p[iprev], &path->p[i]);
2428  tmp = lseg_closept_point(NULL, &lseg, pt);
2429  if (!have_min || float8_lt(tmp, result))
2430  {
2431  result = tmp;
2432  have_min = true;
2433  }
2434  }
2435 
2436  return result;
2437 }
2438 
2439 /*
2440  * Distance from a point to a path
2441  */
2442 Datum
2444 {
2445  Point *pt = PG_GETARG_POINT_P(0);
2446  PATH *path = PG_GETARG_PATH_P(1);
2447 
2449 }
2450 
2451 /*
2452  * Distance from a path to a point
2453  */
2454 Datum
2456 {
2457  PATH *path = PG_GETARG_PATH_P(0);
2458  Point *pt = PG_GETARG_POINT_P(1);
2459 
2461 }
2462 
2463 /*
2464  * Distance from a point to a box
2465  */
2466 Datum
2468 {
2469  Point *pt = PG_GETARG_POINT_P(0);
2470  BOX *box = PG_GETARG_BOX_P(1);
2471 
2472  PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt));
2473 }
2474 
2475 /*
2476  * Distance from a box to a point
2477  */
2478 Datum
2480 {
2481  BOX *box = PG_GETARG_BOX_P(0);
2482  Point *pt = PG_GETARG_POINT_P(1);
2483 
2484  PG_RETURN_FLOAT8(box_closept_point(NULL, box, pt));
2485 }
2486 
2487 /*
2488  * Distance from a lseg to a line
2489  */
2490 Datum
2492 {
2493  LSEG *lseg = PG_GETARG_LSEG_P(0);
2494  LINE *line = PG_GETARG_LINE_P(1);
2495 
2496  PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line));
2497 }
2498 
2499 /*
2500  * Distance from a line to a lseg
2501  */
2502 Datum
2504 {
2505  LINE *line = PG_GETARG_LINE_P(0);
2506  LSEG *lseg = PG_GETARG_LSEG_P(1);
2507 
2508  PG_RETURN_FLOAT8(lseg_closept_line(NULL, lseg, line));
2509 }
2510 
2511 /*
2512  * Distance from a lseg to a box
2513  */
2514 Datum
2516 {
2517  LSEG *lseg = PG_GETARG_LSEG_P(0);
2518  BOX *box = PG_GETARG_BOX_P(1);
2519 
2520  PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg));
2521 }
2522 
2523 /*
2524  * Distance from a box to a lseg
2525  */
2526 Datum
2528 {
2529  BOX *box = PG_GETARG_BOX_P(0);
2530  LSEG *lseg = PG_GETARG_LSEG_P(1);
2531 
2532  PG_RETURN_FLOAT8(box_closept_lseg(NULL, box, lseg));
2533 }
2534 
2535 /*
2536  * Distance from a line to a box
2537  */
2538 Datum
2540 {
2541 #ifdef NOT_USED
2542  LINE *line = PG_GETARG_LINE_P(0);
2543  BOX *box = PG_GETARG_BOX_P(1);
2544 #endif
2545 
2546  /* need to think about this one for a while */
2547  ereport(ERROR,
2548  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2549  errmsg("function \"dist_lb\" not implemented")));
2550 
2551  PG_RETURN_NULL();
2552 }
2553 
2554 /*
2555  * Distance from a box to a line
2556  */
2557 Datum
2559 {
2560 #ifdef NOT_USED
2561  BOX *box = PG_GETARG_BOX_P(0);
2562  LINE *line = PG_GETARG_LINE_P(1);
2563 #endif
2564 
2565  /* need to think about this one for a while */
2566  ereport(ERROR,
2567  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2568  errmsg("function \"dist_bl\" not implemented")));
2569 
2570  PG_RETURN_NULL();
2571 }
2572 
2573 static float8
2575 {
2576  float8 result;
2577 
2578  /* calculate distance to center, and subtract radius */
2579  result = float8_mi(dist_ppoly_internal(&circle->center, poly),
2580  circle->radius);
2581  if (result < 0.0)
2582  result = 0.0;
2583 
2584  return result;
2585 }
2586 
2587 /*
2588  * Distance from a circle to a polygon
2589  */
2590 Datum
2592 {
2593  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
2594  POLYGON *poly = PG_GETARG_POLYGON_P(1);
2595 
2596  PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly));
2597 }
2598 
2599 /*
2600  * Distance from a polygon to a circle
2601  */
2602 Datum
2604 {
2605  POLYGON *poly = PG_GETARG_POLYGON_P(0);
2606  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
2607 
2608  PG_RETURN_FLOAT8(dist_cpoly_internal(circle, poly));
2609 }
2610 
2611 /*
2612  * Distance from a point to a polygon
2613  */
2614 Datum
2616 {
2617  Point *point = PG_GETARG_POINT_P(0);
2618  POLYGON *poly = PG_GETARG_POLYGON_P(1);
2619 
2620  PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly));
2621 }
2622 
2623 Datum
2625 {
2626  POLYGON *poly = PG_GETARG_POLYGON_P(0);
2627  Point *point = PG_GETARG_POINT_P(1);
2628 
2629  PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly));
2630 }
2631 
2632 static float8
2634 {
2635  float8 result;
2636  float8 d;
2637  int i;
2638  LSEG seg;
2639 
2640  if (point_inside(pt, poly->npts, poly->p) != 0)
2641  return 0.0;
2642 
2643  /* initialize distance with segment between first and last points */
2644  seg.p[0].x = poly->p[0].x;
2645  seg.p[0].y = poly->p[0].y;
2646  seg.p[1].x = poly->p[poly->npts - 1].x;
2647  seg.p[1].y = poly->p[poly->npts - 1].y;
2648  result = lseg_closept_point(NULL, &seg, pt);
2649 
2650  /* check distances for other segments */
2651  for (i = 0; i < poly->npts - 1; i++)
2652  {
2653  seg.p[0].x = poly->p[i].x;
2654  seg.p[0].y = poly->p[i].y;
2655  seg.p[1].x = poly->p[i + 1].x;
2656  seg.p[1].y = poly->p[i + 1].y;
2657  d = lseg_closept_point(NULL, &seg, pt);
2658  if (float8_lt(d, result))
2659  result = d;
2660  }
2661 
2662  return result;
2663 }
2664 
2665 
2666 /*---------------------------------------------------------------------
2667  * interpt_
2668  * Intersection point of objects.
2669  * We choose to ignore the "point" of intersection between
2670  * lines and boxes, since there are typically two.
2671  *-------------------------------------------------------------------*/
2672 
2673 /*
2674  * Return whether the line segment intersect with the line. If *result is not
2675  * NULL, it is set to the intersection point.
2676  */
2677 static bool
2678 lseg_interpt_line(Point *result, LSEG *lseg, LINE *line)
2679 {
2680  Point interpt;
2681  LINE tmp;
2682 
2683  /*
2684  * First, we promote the line segment to a line, because we know how to
2685  * find the intersection point of two lines. If they don't have an
2686  * intersection point, we are done.
2687  */
2688  line_construct(&tmp, &lseg->p[0], lseg_sl(lseg));
2689  if (!line_interpt_line(&interpt, &tmp, line))
2690  return false;
2691 
2692  /*
2693  * Then, we check whether the intersection point is actually on the line
2694  * segment.
2695  */
2696  if (!lseg_contain_point(lseg, &interpt))
2697  return false;
2698  if (result != NULL)
2699  {
2700  /*
2701  * If there is an intersection, then check explicitly for matching
2702  * endpoints since there may be rounding effects with annoying LSB
2703  * residue.
2704  */
2705  if (point_eq_point(&lseg->p[0], &interpt))
2706  *result = lseg->p[0];
2707  else if (point_eq_point(&lseg->p[1], &interpt))
2708  *result = lseg->p[1];
2709  else
2710  *result = interpt;
2711  }
2712 
2713  return true;
2714 }
2715 
2716 /*---------------------------------------------------------------------
2717  * close_
2718  * Point of closest proximity between objects.
2719  *-------------------------------------------------------------------*/
2720 
2721 /*
2722  * If *result is not NULL, it is set to the intersection point of a
2723  * perpendicular of the line through the point. Returns the distance
2724  * of those two points.
2725  */
2726 static float8
2727 line_closept_point(Point *result, LINE *line, Point *point)
2728 {
2729  Point closept;
2730  LINE tmp;
2731 
2732  /*
2733  * We drop a perpendicular to find the intersection point. Ordinarily we
2734  * should always find it, but that can fail in the presence of NaN
2735  * coordinates, and perhaps even from simple roundoff issues.
2736  */
2737  line_construct(&tmp, point, line_invsl(line));
2738  if (!line_interpt_line(&closept, &tmp, line))
2739  {
2740  if (result != NULL)
2741  *result = *point;
2742 
2743  return get_float8_nan();
2744  }
2745 
2746  if (result != NULL)
2747  *result = closept;
2748 
2749  return point_dt(&closept, point);
2750 }
2751 
2752 Datum
2754 {
2755  Point *pt = PG_GETARG_POINT_P(0);
2756  LINE *line = PG_GETARG_LINE_P(1);
2757  Point *result;
2758 
2759  result = (Point *) palloc(sizeof(Point));
2760 
2761  if (isnan(line_closept_point(result, line, pt)))
2762  PG_RETURN_NULL();
2763 
2764  PG_RETURN_POINT_P(result);
2765 }
2766 
2767 
2768 /*
2769  * Closest point on line segment to specified point.
2770  *
2771  * If *result is not NULL, set it to the closest point on the line segment
2772  * to the point. Returns the distance of the two points.
2773  */
2774 static float8
2775 lseg_closept_point(Point *result, LSEG *lseg, Point *pt)
2776 {
2777  Point closept;
2778  LINE tmp;
2779 
2780  /*
2781  * To find the closest point, we draw a perpendicular line from the point
2782  * to the line segment.
2783  */
2784  line_construct(&tmp, pt, point_invsl(&lseg->p[0], &lseg->p[1]));
2785  lseg_closept_line(&closept, lseg, &tmp);
2786 
2787  if (result != NULL)
2788  *result = closept;
2789 
2790  return point_dt(&closept, pt);
2791 }
2792 
2793 Datum
2795 {
2796  Point *pt = PG_GETARG_POINT_P(0);
2797  LSEG *lseg = PG_GETARG_LSEG_P(1);
2798  Point *result;
2799 
2800  result = (Point *) palloc(sizeof(Point));
2801 
2802  if (isnan(lseg_closept_point(result, lseg, pt)))
2803  PG_RETURN_NULL();
2804 
2805  PG_RETURN_POINT_P(result);
2806 }
2807 
2808 
2809 /*
2810  * Closest point on line segment to line segment
2811  */
2812 static float8
2813 lseg_closept_lseg(Point *result, LSEG *on_lseg, LSEG *to_lseg)
2814 {
2815  Point point;
2816  float8 dist,
2817  d;
2818 
2819  /* First, we handle the case when the line segments are intersecting. */
2820  if (lseg_interpt_lseg(result, on_lseg, to_lseg))
2821  return 0.0;
2822 
2823  /*
2824  * Then, we find the closest points from the endpoints of the second line
2825  * segment, and keep the closest one.
2826  */
2827  dist = lseg_closept_point(result, on_lseg, &to_lseg->p[0]);
2828  d = lseg_closept_point(&point, on_lseg, &to_lseg->p[1]);
2829  if (float8_lt(d, dist))
2830  {
2831  dist = d;
2832  if (result != NULL)
2833  *result = point;
2834  }
2835 
2836  /* The closest point can still be one of the endpoints, so we test them. */
2837  d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[0]);
2838  if (float8_lt(d, dist))
2839  {
2840  dist = d;
2841  if (result != NULL)
2842  *result = on_lseg->p[0];
2843  }
2844  d = lseg_closept_point(NULL, to_lseg, &on_lseg->p[1]);
2845  if (float8_lt(d, dist))
2846  {
2847  dist = d;
2848  if (result != NULL)
2849  *result = on_lseg->p[1];
2850  }
2851 
2852  return dist;
2853 }
2854 
2855 Datum
2857 {
2858  LSEG *l1 = PG_GETARG_LSEG_P(0);
2859  LSEG *l2 = PG_GETARG_LSEG_P(1);
2860  Point *result;
2861 
2862  if (lseg_sl(l1) == lseg_sl(l2))
2863  PG_RETURN_NULL();
2864 
2865  result = (Point *) palloc(sizeof(Point));
2866 
2867  if (isnan(lseg_closept_lseg(result, l2, l1)))
2868  PG_RETURN_NULL();
2869 
2870  PG_RETURN_POINT_P(result);
2871 }
2872 
2873 
2874 /*
2875  * Closest point on or in box to specified point.
2876  *
2877  * If *result is not NULL, set it to the closest point on the box to the
2878  * given point, and return the distance of the two points.
2879  */
2880 static float8
2881 box_closept_point(Point *result, BOX *box, Point *pt)
2882 {
2883  float8 dist,
2884  d;
2885  Point point,
2886  closept;
2887  LSEG lseg;
2888 
2889  if (box_contain_point(box, pt))
2890  {
2891  if (result != NULL)
2892  *result = *pt;
2893 
2894  return 0.0;
2895  }
2896 
2897  /* pairwise check lseg distances */
2898  point.x = box->low.x;
2899  point.y = box->high.y;
2900  statlseg_construct(&lseg, &box->low, &point);
2901  dist = lseg_closept_point(result, &lseg, pt);
2902 
2903  statlseg_construct(&lseg, &box->high, &point);
2904  d = lseg_closept_point(&closept, &lseg, pt);
2905  if (float8_lt(d, dist))
2906  {
2907  dist = d;
2908  if (result != NULL)
2909  *result = closept;
2910  }
2911 
2912  point.x = box->high.x;
2913  point.y = box->low.y;
2914  statlseg_construct(&lseg, &box->low, &point);
2915  d = lseg_closept_point(&closept, &lseg, pt);
2916  if (float8_lt(d, dist))
2917  {
2918  dist = d;
2919  if (result != NULL)
2920  *result = closept;
2921  }
2922 
2923  statlseg_construct(&lseg, &box->high, &point);
2924  d = lseg_closept_point(&closept, &lseg, pt);
2925  if (float8_lt(d, dist))
2926  {
2927  dist = d;
2928  if (result != NULL)
2929  *result = closept;
2930  }
2931 
2932  return dist;
2933 }
2934 
2935 Datum
2937 {
2938  Point *pt = PG_GETARG_POINT_P(0);
2939  BOX *box = PG_GETARG_BOX_P(1);
2940  Point *result;
2941 
2942  result = (Point *) palloc(sizeof(Point));
2943 
2944  if (isnan(box_closept_point(result, box, pt)))
2945  PG_RETURN_NULL();
2946 
2947  PG_RETURN_POINT_P(result);
2948 }
2949 
2950 
2951 /* close_sl()
2952  * Closest point on line to line segment.
2953  *
2954  * XXX THIS CODE IS WRONG
2955  * The code is actually calculating the point on the line segment
2956  * which is backwards from the routine naming convention.
2957  * Copied code to new routine close_ls() but haven't fixed this one yet.
2958  * - thomas 1998-01-31
2959  */
2960 Datum
2962 {
2963 #ifdef NOT_USED
2964  LSEG *lseg = PG_GETARG_LSEG_P(0);
2965  LINE *line = PG_GETARG_LINE_P(1);
2966  Point *result;
2967  float8 d1,
2968  d2;
2969 
2970  result = (Point *) palloc(sizeof(Point));
2971 
2972  if (lseg_interpt_line(result, lseg, line))
2973  PG_RETURN_POINT_P(result);
2974 
2975  d1 = line_closept_point(NULL, line, &lseg->p[0]);
2976  d2 = line_closept_point(NULL, line, &lseg->p[1]);
2977  if (float8_lt(d1, d2))
2978  *result = lseg->p[0];
2979  else
2980  *result = lseg->p[1];
2981 
2982  PG_RETURN_POINT_P(result);
2983 #endif
2984 
2985  ereport(ERROR,
2986  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
2987  errmsg("function \"close_sl\" not implemented")));
2988 
2989  PG_RETURN_NULL();
2990 }
2991 
2992 /*
2993  * Closest point on line segment to line.
2994  *
2995  * Return the distance between the line and the closest point of the line
2996  * segment to the line. If *result is not NULL, set it to that point.
2997  *
2998  * NOTE: When the lines are parallel, endpoints of one of the line segment
2999  * are FPeq(), in presence of NaN or Infinite coordinates, or perhaps =
3000  * even because of simple roundoff issues, there may not be a single closest
3001  * point. We are likely to set the result to the second endpoint in these
3002  * cases.
3003  */
3004 static float8
3005 lseg_closept_line(Point *result, LSEG *lseg, LINE *line)
3006 {
3007  float8 dist1,
3008  dist2;
3009 
3010  if (lseg_interpt_line(result, lseg, line))
3011  return 0.0;
3012 
3013  dist1 = line_closept_point(NULL, line, &lseg->p[0]);
3014  dist2 = line_closept_point(NULL, line, &lseg->p[1]);
3015 
3016  if (dist1 < dist2)
3017  {
3018  if (result != NULL)
3019  *result = lseg->p[0];
3020 
3021  return dist1;
3022  }
3023  else
3024  {
3025  if (result != NULL)
3026  *result = lseg->p[1];
3027 
3028  return dist2;
3029  }
3030 }
3031 
3032 Datum
3034 {
3035  LINE *line = PG_GETARG_LINE_P(0);
3036  LSEG *lseg = PG_GETARG_LSEG_P(1);
3037  Point *result;
3038 
3039  if (lseg_sl(lseg) == line_sl(line))
3040  PG_RETURN_NULL();
3041 
3042  result = (Point *) palloc(sizeof(Point));
3043 
3044  if (isnan(lseg_closept_line(result, lseg, line)))
3045  PG_RETURN_NULL();
3046 
3047  PG_RETURN_POINT_P(result);
3048 }
3049 
3050 
3051 /*
3052  * Closest point on or in box to line segment.
3053  *
3054  * Returns the distance between the closest point on or in the box to
3055  * the line segment. If *result is not NULL, it is set to that point.
3056  */
3057 static float8
3058 box_closept_lseg(Point *result, BOX *box, LSEG *lseg)
3059 {
3060  float8 dist,
3061  d;
3062  Point point,
3063  closept;
3064  LSEG bseg;
3065 
3066  if (box_interpt_lseg(result, box, lseg))
3067  return 0.0;
3068 
3069  /* pairwise check lseg distances */
3070  point.x = box->low.x;
3071  point.y = box->high.y;
3072  statlseg_construct(&bseg, &box->low, &point);
3073  dist = lseg_closept_lseg(result, &bseg, lseg);
3074 
3075  statlseg_construct(&bseg, &box->high, &point);
3076  d = lseg_closept_lseg(&closept, &bseg, lseg);
3077  if (float8_lt(d, dist))
3078  {
3079  dist = d;
3080  if (result != NULL)
3081  *result = closept;
3082  }
3083 
3084  point.x = box->high.x;
3085  point.y = box->low.y;
3086  statlseg_construct(&bseg, &box->low, &point);
3087  d = lseg_closept_lseg(&closept, &bseg, lseg);
3088  if (float8_lt(d, dist))
3089  {
3090  dist = d;
3091  if (result != NULL)
3092  *result = closept;
3093  }
3094 
3095  statlseg_construct(&bseg, &box->high, &point);
3096  d = lseg_closept_lseg(&closept, &bseg, lseg);
3097  if (float8_lt(d, dist))
3098  {
3099  dist = d;
3100  if (result != NULL)
3101  *result = closept;
3102  }
3103 
3104  return dist;
3105 }
3106 
3107 Datum
3109 {
3110  LSEG *lseg = PG_GETARG_LSEG_P(0);
3111  BOX *box = PG_GETARG_BOX_P(1);
3112  Point *result;
3113 
3114  result = (Point *) palloc(sizeof(Point));
3115 
3116  if (isnan(box_closept_lseg(result, box, lseg)))
3117  PG_RETURN_NULL();
3118 
3119  PG_RETURN_POINT_P(result);
3120 }
3121 
3122 
3123 Datum
3125 {
3126 #ifdef NOT_USED
3127  LINE *line = PG_GETARG_LINE_P(0);
3128  BOX *box = PG_GETARG_BOX_P(1);
3129 #endif
3130 
3131  /* think about this one for a while */
3132  ereport(ERROR,
3133  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
3134  errmsg("function \"close_lb\" not implemented")));
3135 
3136  PG_RETURN_NULL();
3137 }
3138 
3139 /*---------------------------------------------------------------------
3140  * on_
3141  * Whether one object lies completely within another.
3142  *-------------------------------------------------------------------*/
3143 
3144 /*
3145  * Does the point satisfy the equation?
3146  */
3147 static bool
3149 {
3150  return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x),
3151  float8_mul(line->B, point->y)),
3152  line->C));
3153 }
3154 
3155 Datum
3157 {
3158  Point *pt = PG_GETARG_POINT_P(0);
3159  LINE *line = PG_GETARG_LINE_P(1);
3160 
3162 }
3163 
3164 
3165 /*
3166  * Determine colinearity by detecting a triangle inequality.
3167  * This algorithm seems to behave nicely even with lsb residues - tgl 1997-07-09
3168  */
3169 static bool
3171 {
3172  return FPeq(point_dt(pt, &lseg->p[0]) +
3173  point_dt(pt, &lseg->p[1]),
3174  point_dt(&lseg->p[0], &lseg->p[1]));
3175 }
3176 
3177 Datum
3179 {
3180  Point *pt = PG_GETARG_POINT_P(0);
3181  LSEG *lseg = PG_GETARG_LSEG_P(1);
3182 
3184 }
3185 
3186 
3187 /*
3188  * Check whether the point is in the box or on its border
3189  */
3190 static bool
3192 {
3193  return box->high.x >= point->x && box->low.x <= point->x &&
3194  box->high.y >= point->y && box->low.y <= point->y;
3195 }
3196 
3197 Datum
3199 {
3200  Point *pt = PG_GETARG_POINT_P(0);
3201  BOX *box = PG_GETARG_BOX_P(1);
3202 
3204 }
3205 
3206 Datum
3208 {
3209  BOX *box = PG_GETARG_BOX_P(0);
3210  Point *pt = PG_GETARG_POINT_P(1);
3211 
3213 }
3214 
3215 /* on_ppath -
3216  * Whether a point lies within (on) a polyline.
3217  * If open, we have to (groan) check each segment.
3218  * (uses same algorithm as for point intersecting segment - tgl 1997-07-09)
3219  * If closed, we use the old O(n) ray method for point-in-polygon.
3220  * The ray is horizontal, from pt out to the right.
3221  * Each segment that crosses the ray counts as an
3222  * intersection; note that an endpoint or edge may touch
3223  * but not cross.
3224  * (we can do p-in-p in lg(n), but it takes preprocessing)
3225  */
3226 Datum
3228 {
3229  Point *pt = PG_GETARG_POINT_P(0);
3230  PATH *path = PG_GETARG_PATH_P(1);
3231  int i,
3232  n;
3233  float8 a,
3234  b;
3235 
3236  /*-- OPEN --*/
3237  if (!path->closed)
3238  {
3239  n = path->npts - 1;
3240  a = point_dt(pt, &path->p[0]);
3241  for (i = 0; i < n; i++)
3242  {
3243  b = point_dt(pt, &path->p[i + 1]);
3244  if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1])))
3245  PG_RETURN_BOOL(true);
3246  a = b;
3247  }
3248  PG_RETURN_BOOL(false);
3249  }
3250 
3251  /*-- CLOSED --*/
3252  PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0);
3253 }
3254 
3255 
3256 /*
3257  * Check whether the line segment is on the line or close enough
3258  *
3259  * It is, if both of its points are on the line or close enough.
3260  */
3261 Datum
3263 {
3264  LSEG *lseg = PG_GETARG_LSEG_P(0);
3265  LINE *line = PG_GETARG_LINE_P(1);
3266 
3267  PG_RETURN_BOOL(line_contain_point(line, &lseg->p[0]) &&
3268  line_contain_point(line, &lseg->p[1]));
3269 }
3270 
3271 
3272 /*
3273  * Check whether the line segment is in the box or on its border
3274  *
3275  * It is, if both of its points are in the box or on its border.
3276  */
3277 static bool
3279 {
3280  return box_contain_point(box, &lseg->p[0]) &&
3281  box_contain_point(box, &lseg->p[1]);
3282 }
3283 
3284 Datum
3286 {
3287  LSEG *lseg = PG_GETARG_LSEG_P(0);
3288  BOX *box = PG_GETARG_BOX_P(1);
3289 
3290  PG_RETURN_BOOL(box_contain_lseg(box, lseg));
3291 }
3292 
3293 /*---------------------------------------------------------------------
3294  * inter_
3295  * Whether one object intersects another.
3296  *-------------------------------------------------------------------*/
3297 
3298 Datum
3300 {
3301  LSEG *lseg = PG_GETARG_LSEG_P(0);
3302  LINE *line = PG_GETARG_LINE_P(1);
3303 
3304  PG_RETURN_BOOL(lseg_interpt_line(NULL, lseg, line));
3305 }
3306 
3307 
3308 /*
3309  * Do line segment and box intersect?
3310  *
3311  * Segment completely inside box counts as intersection.
3312  * If you want only segments crossing box boundaries,
3313  * try converting box to path first.
3314  *
3315  * This function also sets the *result to the closest point on the line
3316  * segment to the center of the box when they overlap and the result is
3317  * not NULL. It is somewhat arbitrary, but maybe the best we can do as
3318  * there are typically two points they intersect.
3319  *
3320  * Optimize for non-intersection by checking for box intersection first.
3321  * - thomas 1998-01-30
3322  */
3323 static bool
3324 box_interpt_lseg(Point *result, BOX *box, LSEG *lseg)
3325 {
3326  BOX lbox;
3327  LSEG bseg;
3328  Point point;
3329 
3330  lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x);
3331  lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y);
3332  lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x);
3333  lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y);
3334 
3335  /* nothing close to overlap? then not going to intersect */
3336  if (!box_ov(&lbox, box))
3337  return false;
3338 
3339  if (result != NULL)
3340  {
3341  box_cn(&point, box);
3342  lseg_closept_point(result, lseg, &point);
3343  }
3344 
3345  /* an endpoint of segment is inside box? then clearly intersects */
3346  if (box_contain_point(box, &lseg->p[0]) ||
3347  box_contain_point(box, &lseg->p[1]))
3348  return true;
3349 
3350  /* pairwise check lseg intersections */
3351  point.x = box->low.x;
3352  point.y = box->high.y;
3353  statlseg_construct(&bseg, &box->low, &point);
3354  if (lseg_interpt_lseg(NULL, &bseg, lseg))
3355  return true;
3356 
3357  statlseg_construct(&bseg, &box->high, &point);
3358  if (lseg_interpt_lseg(NULL, &bseg, lseg))
3359  return true;
3360 
3361  point.x = box->high.x;
3362  point.y = box->low.y;
3363  statlseg_construct(&bseg, &box->low, &point);
3364  if (lseg_interpt_lseg(NULL, &bseg, lseg))
3365  return true;
3366 
3367  statlseg_construct(&bseg, &box->high, &point);
3368  if (lseg_interpt_lseg(NULL, &bseg, lseg))
3369  return true;
3370 
3371  /* if we dropped through, no two segs intersected */
3372  return false;
3373 }
3374 
3375 Datum
3377 {
3378  LSEG *lseg = PG_GETARG_LSEG_P(0);
3379  BOX *box = PG_GETARG_BOX_P(1);
3380 
3381  PG_RETURN_BOOL(box_interpt_lseg(NULL, box, lseg));
3382 }
3383 
3384 
3385 /* inter_lb()
3386  * Do line and box intersect?
3387  */
3388 Datum
3390 {
3391  LINE *line = PG_GETARG_LINE_P(0);
3392  BOX *box = PG_GETARG_BOX_P(1);
3393  LSEG bseg;
3394  Point p1,
3395  p2;
3396 
3397  /* pairwise check lseg intersections */
3398  p1.x = box->low.x;
3399  p1.y = box->low.y;
3400  p2.x = box->low.x;
3401  p2.y = box->high.y;
3402  statlseg_construct(&bseg, &p1, &p2);
3403  if (lseg_interpt_line(NULL, &bseg, line))
3404  PG_RETURN_BOOL(true);
3405  p1.x = box->high.x;
3406  p1.y = box->high.y;
3407  statlseg_construct(&bseg, &p1, &p2);
3408  if (lseg_interpt_line(NULL, &bseg, line))
3409  PG_RETURN_BOOL(true);
3410  p2.x = box->high.x;
3411  p2.y = box->low.y;
3412  statlseg_construct(&bseg, &p1, &p2);
3413  if (lseg_interpt_line(NULL, &bseg, line))
3414  PG_RETURN_BOOL(true);
3415  p1.x = box->low.x;
3416  p1.y = box->low.y;
3417  statlseg_construct(&bseg, &p1, &p2);
3418  if (lseg_interpt_line(NULL, &bseg, line))
3419  PG_RETURN_BOOL(true);
3420 
3421  /* if we dropped through, no intersection */
3422  PG_RETURN_BOOL(false);
3423 }
3424 
3425 /*------------------------------------------------------------------
3426  * The following routines define a data type and operator class for
3427  * POLYGONS .... Part of which (the polygon's bounding box) is built on
3428  * top of the BOX data type.
3429  *
3430  * make_bound_box - create the bounding box for the input polygon
3431  *------------------------------------------------------------------*/
3432 
3433 /*---------------------------------------------------------------------
3434  * Make the smallest bounding box for the given polygon.
3435  *---------------------------------------------------------------------*/
3436 static void
3438 {
3439  int i;
3440  float8 x1,
3441  y1,
3442  x2,
3443  y2;
3444 
3445  Assert(poly->npts > 0);
3446 
3447  x1 = x2 = poly->p[0].x;
3448  y2 = y1 = poly->p[0].y;
3449  for (i = 1; i < poly->npts; i++)
3450  {
3451  if (float8_lt(poly->p[i].x, x1))
3452  x1 = poly->p[i].x;
3453  if (float8_gt(poly->p[i].x, x2))
3454  x2 = poly->p[i].x;
3455  if (float8_lt(poly->p[i].y, y1))
3456  y1 = poly->p[i].y;
3457  if (float8_gt(poly->p[i].y, y2))
3458  y2 = poly->p[i].y;
3459  }
3460 
3461  poly->boundbox.low.x = x1;
3462  poly->boundbox.high.x = x2;
3463  poly->boundbox.low.y = y1;
3464  poly->boundbox.high.y = y2;
3465 }
3466 
3467 /*------------------------------------------------------------------
3468  * poly_in - read in the polygon from a string specification
3469  *
3470  * External format:
3471  * "((x0,y0),...,(xn,yn))"
3472  * "x0,y0,...,xn,yn"
3473  * also supports the older style "(x1,...,xn,y1,...yn)"
3474  *------------------------------------------------------------------*/
3475 Datum
3477 {
3478  char *str = PG_GETARG_CSTRING(0);
3479  POLYGON *poly;
3480  int npts;
3481  int size;
3482  int base_size;
3483  bool isopen;
3484 
3485  if ((npts = pair_count(str, ',')) <= 0)
3486  ereport(ERROR,
3487  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
3488  errmsg("invalid input syntax for type %s: \"%s\"",
3489  "polygon", str)));
3490 
3491  base_size = sizeof(poly->p[0]) * npts;
3492  size = offsetof(POLYGON, p) + base_size;
3493 
3494  /* Check for integer overflow */
3495  if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
3496  ereport(ERROR,
3497  (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
3498  errmsg("too many points requested")));
3499 
3500  poly = (POLYGON *) palloc0(size); /* zero any holes */
3501 
3502  SET_VARSIZE(poly, size);
3503  poly->npts = npts;
3504 
3505  path_decode(str, false, npts, &(poly->p[0]), &isopen, NULL, "polygon", str);
3506 
3507  make_bound_box(poly);
3508 
3509  PG_RETURN_POLYGON_P(poly);
3510 }
3511 
3512 /*---------------------------------------------------------------
3513  * poly_out - convert internal POLYGON representation to the
3514  * character string format "((f8,f8),...,(f8,f8))"
3515  *---------------------------------------------------------------*/
3516 Datum
3518 {
3519  POLYGON *poly = PG_GETARG_POLYGON_P(0);
3520 
3522 }
3523 
3524 /*
3525  * poly_recv - converts external binary format to polygon
3526  *
3527  * External representation is int32 number of points, and the points.
3528  * We recompute the bounding box on read, instead of trusting it to
3529  * be valid. (Checking it would take just as long, so may as well
3530  * omit it from external representation.)
3531  */
3532 Datum
3534 {
3536  POLYGON *poly;
3537  int32 npts;
3538  int32 i;
3539  int size;
3540 
3541  npts = pq_getmsgint(buf, sizeof(int32));
3542  if (npts <= 0 || npts >= (int32) ((INT_MAX - offsetof(POLYGON, p)) / sizeof(Point)))
3543  ereport(ERROR,
3544  (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
3545  errmsg("invalid number of points in external \"polygon\" value")));
3546 
3547  size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * npts;
3548  poly = (POLYGON *) palloc0(size); /* zero any holes */
3549 
3550  SET_VARSIZE(poly, size);
3551  poly->npts = npts;
3552 
3553  for (i = 0; i < npts; i++)
3554  {
3555  poly->p[i].x = pq_getmsgfloat8(buf);
3556  poly->p[i].y = pq_getmsgfloat8(buf);
3557  }
3558 
3559  make_bound_box(poly);
3560 
3561  PG_RETURN_POLYGON_P(poly);
3562 }
3563 
3564 /*
3565  * poly_send - converts polygon to binary format
3566  */
3567 Datum
3569 {
3570  POLYGON *poly = PG_GETARG_POLYGON_P(0);
3572  int32 i;
3573 
3574  pq_begintypsend(&buf);
3575  pq_sendint32(&buf, poly->npts);
3576  for (i = 0; i < poly->npts; i++)
3577  {
3578  pq_sendfloat8(&buf, poly->p[i].x);
3579  pq_sendfloat8(&buf, poly->p[i].y);
3580  }
3582 }
3583 
3584 
3585 /*-------------------------------------------------------
3586  * Is polygon A strictly left of polygon B? i.e. is
3587  * the right most point of A left of the left most point
3588  * of B?
3589  *-------------------------------------------------------*/
3590 Datum
3592 {
3593  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3594  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3595  bool result;
3596 
3597  result = polya->boundbox.high.x < polyb->boundbox.low.x;
3598 
3599  /*
3600  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3601  */
3602  PG_FREE_IF_COPY(polya, 0);
3603  PG_FREE_IF_COPY(polyb, 1);
3604 
3605  PG_RETURN_BOOL(result);
3606 }
3607 
3608 /*-------------------------------------------------------
3609  * Is polygon A overlapping or left of polygon B? i.e. is
3610  * the right most point of A at or left of the right most point
3611  * of B?
3612  *-------------------------------------------------------*/
3613 Datum
3615 {
3616  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3617  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3618  bool result;
3619 
3620  result = polya->boundbox.high.x <= polyb->boundbox.high.x;
3621 
3622  /*
3623  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3624  */
3625  PG_FREE_IF_COPY(polya, 0);
3626  PG_FREE_IF_COPY(polyb, 1);
3627 
3628  PG_RETURN_BOOL(result);
3629 }
3630 
3631 /*-------------------------------------------------------
3632  * Is polygon A strictly right of polygon B? i.e. is
3633  * the left most point of A right of the right most point
3634  * of B?
3635  *-------------------------------------------------------*/
3636 Datum
3638 {
3639  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3640  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3641  bool result;
3642 
3643  result = polya->boundbox.low.x > polyb->boundbox.high.x;
3644 
3645  /*
3646  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3647  */
3648  PG_FREE_IF_COPY(polya, 0);
3649  PG_FREE_IF_COPY(polyb, 1);
3650 
3651  PG_RETURN_BOOL(result);
3652 }
3653 
3654 /*-------------------------------------------------------
3655  * Is polygon A overlapping or right of polygon B? i.e. is
3656  * the left most point of A at or right of the left most point
3657  * of B?
3658  *-------------------------------------------------------*/
3659 Datum
3661 {
3662  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3663  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3664  bool result;
3665 
3666  result = polya->boundbox.low.x >= polyb->boundbox.low.x;
3667 
3668  /*
3669  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3670  */
3671  PG_FREE_IF_COPY(polya, 0);
3672  PG_FREE_IF_COPY(polyb, 1);
3673 
3674  PG_RETURN_BOOL(result);
3675 }
3676 
3677 /*-------------------------------------------------------
3678  * Is polygon A strictly below polygon B? i.e. is
3679  * the upper most point of A below the lower most point
3680  * of B?
3681  *-------------------------------------------------------*/
3682 Datum
3684 {
3685  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3686  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3687  bool result;
3688 
3689  result = polya->boundbox.high.y < polyb->boundbox.low.y;
3690 
3691  /*
3692  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3693  */
3694  PG_FREE_IF_COPY(polya, 0);
3695  PG_FREE_IF_COPY(polyb, 1);
3696 
3697  PG_RETURN_BOOL(result);
3698 }
3699 
3700 /*-------------------------------------------------------
3701  * Is polygon A overlapping or below polygon B? i.e. is
3702  * the upper most point of A at or below the upper most point
3703  * of B?
3704  *-------------------------------------------------------*/
3705 Datum
3707 {
3708  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3709  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3710  bool result;
3711 
3712  result = polya->boundbox.high.y <= polyb->boundbox.high.y;
3713 
3714  /*
3715  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3716  */
3717  PG_FREE_IF_COPY(polya, 0);
3718  PG_FREE_IF_COPY(polyb, 1);
3719 
3720  PG_RETURN_BOOL(result);
3721 }
3722 
3723 /*-------------------------------------------------------
3724  * Is polygon A strictly above polygon B? i.e. is
3725  * the lower most point of A above the upper most point
3726  * of B?
3727  *-------------------------------------------------------*/
3728 Datum
3730 {
3731  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3732  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3733  bool result;
3734 
3735  result = polya->boundbox.low.y > polyb->boundbox.high.y;
3736 
3737  /*
3738  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3739  */
3740  PG_FREE_IF_COPY(polya, 0);
3741  PG_FREE_IF_COPY(polyb, 1);
3742 
3743  PG_RETURN_BOOL(result);
3744 }
3745 
3746 /*-------------------------------------------------------
3747  * Is polygon A overlapping or above polygon B? i.e. is
3748  * the lower most point of A at or above the lower most point
3749  * of B?
3750  *-------------------------------------------------------*/
3751 Datum
3753 {
3754  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3755  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3756  bool result;
3757 
3758  result = polya->boundbox.low.y >= polyb->boundbox.low.y;
3759 
3760  /*
3761  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3762  */
3763  PG_FREE_IF_COPY(polya, 0);
3764  PG_FREE_IF_COPY(polyb, 1);
3765 
3766  PG_RETURN_BOOL(result);
3767 }
3768 
3769 
3770 /*-------------------------------------------------------
3771  * Is polygon A the same as polygon B? i.e. are all the
3772  * points the same?
3773  * Check all points for matches in both forward and reverse
3774  * direction since polygons are non-directional and are
3775  * closed shapes.
3776  *-------------------------------------------------------*/
3777 Datum
3779 {
3780  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3781  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3782  bool result;
3783 
3784  if (polya->npts != polyb->npts)
3785  result = false;
3786  else
3787  result = plist_same(polya->npts, polya->p, polyb->p);
3788 
3789  /*
3790  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3791  */
3792  PG_FREE_IF_COPY(polya, 0);
3793  PG_FREE_IF_COPY(polyb, 1);
3794 
3795  PG_RETURN_BOOL(result);
3796 }
3797 
3798 /*-----------------------------------------------------------------
3799  * Determine if polygon A overlaps polygon B
3800  *-----------------------------------------------------------------*/
3801 Datum
3803 {
3804  POLYGON *polya = PG_GETARG_POLYGON_P(0);
3805  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
3806  bool result;
3807 
3808  Assert(polya->npts > 0 && polyb->npts > 0);
3809 
3810  /* Quick check by bounding box */
3811  result = box_ov(&polya->boundbox, &polyb->boundbox);
3812 
3813  /*
3814  * Brute-force algorithm - try to find intersected edges, if so then
3815  * polygons are overlapped else check is one polygon inside other or not
3816  * by testing single point of them.
3817  */
3818  if (result)
3819  {
3820  int ia,
3821  ib;
3822  LSEG sa,
3823  sb;
3824 
3825  /* Init first of polya's edge with last point */
3826  sa.p[0] = polya->p[polya->npts - 1];
3827  result = false;
3828 
3829  for (ia = 0; ia < polya->npts && !result; ia++)
3830  {
3831  /* Second point of polya's edge is a current one */
3832  sa.p[1] = polya->p[ia];
3833 
3834  /* Init first of polyb's edge with last point */
3835  sb.p[0] = polyb->p[polyb->npts - 1];
3836 
3837  for (ib = 0; ib < polyb->npts && !result; ib++)
3838  {
3839  sb.p[1] = polyb->p[ib];
3840  result = lseg_interpt_lseg(NULL, &sa, &sb);
3841  sb.p[0] = sb.p[1];
3842  }
3843 
3844  /*
3845  * move current endpoint to the first point of next edge
3846  */
3847  sa.p[0] = sa.p[1];
3848  }
3849 
3850  if (!result)
3851  {
3852  result = (point_inside(polya->p, polyb->npts, polyb->p) ||
3853  point_inside(polyb->p, polya->npts, polya->p));
3854  }
3855  }
3856 
3857  /*
3858  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
3859  */
3860  PG_FREE_IF_COPY(polya, 0);
3861  PG_FREE_IF_COPY(polyb, 1);
3862 
3863  PG_RETURN_BOOL(result);
3864 }
3865 
3866 /*
3867  * Tests special kind of segment for in/out of polygon.
3868  * Special kind means:
3869  * - point a should be on segment s
3870  * - segment (a,b) should not be contained by s
3871  * Returns true if:
3872  * - segment (a,b) is collinear to s and (a,b) is in polygon
3873  * - segment (a,b) s not collinear to s. Note: that doesn't
3874  * mean that segment is in polygon!
3875  */
3876 
3877 static bool
3878 touched_lseg_inside_poly(Point *a, Point *b, LSEG *s, POLYGON *poly, int start)
3879 {
3880  /* point a is on s, b is not */
3881  LSEG t;
3882 
3883  t.p[0] = *a;
3884  t.p[1] = *b;
3885 
3886  if (point_eq_point(a, s->p))
3887  {
3888  if (lseg_contain_point(&t, s->p + 1))
3889  return lseg_inside_poly(b, s->p + 1, poly, start);
3890  }
3891  else if (point_eq_point(a, s->p + 1))
3892  {
3893  if (lseg_contain_point(&t, s->p))
3894  return lseg_inside_poly(b, s->p, poly, start);
3895  }
3896  else if (lseg_contain_point(&t, s->p))
3897  {
3898  return lseg_inside_poly(b, s->p, poly, start);
3899  }
3900  else if (lseg_contain_point(&t, s->p + 1))
3901  {
3902  return lseg_inside_poly(b, s->p + 1, poly, start);
3903  }
3904 
3905  return true; /* may be not true, but that will check later */
3906 }
3907 
3908 /*
3909  * Returns true if segment (a,b) is in polygon, option
3910  * start is used for optimization - function checks
3911  * polygon's edges starting from start
3912  */
3913 static bool
3914 lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start)
3915 {
3916  LSEG s,
3917  t;
3918  int i;
3919  bool res = true,
3920  intersection = false;
3921 
3922  t.p[0] = *a;
3923  t.p[1] = *b;
3924  s.p[0] = poly->p[(start == 0) ? (poly->npts - 1) : (start - 1)];
3925 
3926  for (i = start; i < poly->npts && res; i++)
3927  {
3928  Point interpt;
3929 
3931 
3932  s.p[1] = poly->p[i];
3933 
3934  if (lseg_contain_point(&s, t.p))
3935  {
3936  if (lseg_contain_point(&s, t.p + 1))
3937  return true; /* t is contained by s */
3938 
3939  /* Y-cross */
3940  res = touched_lseg_inside_poly(t.p, t.p + 1, &s, poly, i + 1);
3941  }
3942  else if (lseg_contain_point(&s, t.p + 1))
3943  {
3944  /* Y-cross */
3945  res = touched_lseg_inside_poly(t.p + 1, t.p, &s, poly, i + 1);
3946  }
3947  else if (lseg_interpt_lseg(&interpt, &t, &s))
3948  {
3949  /*
3950  * segments are X-crossing, go to check each subsegment
3951  */
3952 
3953  intersection = true;
3954  res = lseg_inside_poly(t.p, &interpt, poly, i + 1);
3955  if (res)
3956  res = lseg_inside_poly(t.p + 1, &interpt, poly, i + 1);
3957  }
3958 
3959  s.p[0] = s.p[1];
3960  }
3961 
3962  if (res && !intersection)
3963  {
3964  Point p;
3965 
3966  /*
3967  * if X-intersection wasn't found then check central point of tested
3968  * segment. In opposite case we already check all subsegments
3969  */
3970  p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0);
3971  p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0);
3972 
3973  res = point_inside(&p, poly->npts, poly->p);
3974  }
3975 
3976  return res;
3977 }
3978 
3979 /*
3980  * Check whether the first polygon contains the second
3981  */
3982 static bool
3983 poly_contain_poly(POLYGON *contains_poly, POLYGON *contained_poly)
3984 {
3985  int i;
3986  LSEG s;
3987 
3988  Assert(contains_poly->npts > 0 && contained_poly->npts > 0);
3989 
3990  /*
3991  * Quick check to see if contained's bounding box is contained in
3992  * contains' bb.
3993  */
3994  if (!box_contain_box(&contains_poly->boundbox, &contained_poly->boundbox))
3995  return false;
3996 
3997  s.p[0] = contained_poly->p[contained_poly->npts - 1];
3998 
3999  for (i = 0; i < contained_poly->npts; i++)
4000  {
4001  s.p[1] = contained_poly->p[i];
4002  if (!lseg_inside_poly(s.p, s.p + 1, contains_poly, 0))
4003  return false;
4004  s.p[0] = s.p[1];
4005  }
4006 
4007  return true;
4008 }
4009 
4010 Datum
4012 {
4013  POLYGON *polya = PG_GETARG_POLYGON_P(0);
4014  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
4015  bool result;
4016 
4017  result = poly_contain_poly(polya, polyb);
4018 
4019  /*
4020  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
4021  */
4022  PG_FREE_IF_COPY(polya, 0);
4023  PG_FREE_IF_COPY(polyb, 1);
4024 
4025  PG_RETURN_BOOL(result);
4026 }
4027 
4028 
4029 /*-----------------------------------------------------------------
4030  * Determine if polygon A is contained by polygon B
4031  *-----------------------------------------------------------------*/
4032 Datum
4034 {
4035  POLYGON *polya = PG_GETARG_POLYGON_P(0);
4036  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
4037  bool result;
4038 
4039  /* Just switch the arguments and pass it off to poly_contain */
4040  result = poly_contain_poly(polyb, polya);
4041 
4042  /*
4043  * Avoid leaking memory for toasted inputs ... needed for rtree indexes
4044  */
4045  PG_FREE_IF_COPY(polya, 0);
4046  PG_FREE_IF_COPY(polyb, 1);
4047 
4048  PG_RETURN_BOOL(result);
4049 }
4050 
4051 
4052 Datum
4054 {
4055  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4056  Point *p = PG_GETARG_POINT_P(1);
4057 
4058  PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
4059 }
4060 
4061 Datum
4063 {
4064  Point *p = PG_GETARG_POINT_P(0);
4065  POLYGON *poly = PG_GETARG_POLYGON_P(1);
4066 
4067  PG_RETURN_BOOL(point_inside(p, poly->npts, poly->p) != 0);
4068 }
4069 
4070 
4071 Datum
4073 {
4074 #ifdef NOT_USED
4075  POLYGON *polya = PG_GETARG_POLYGON_P(0);
4076  POLYGON *polyb = PG_GETARG_POLYGON_P(1);
4077 #endif
4078 
4079  ereport(ERROR,
4080  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4081  errmsg("function \"poly_distance\" not implemented")));
4082 
4083  PG_RETURN_NULL();
4084 }
4085 
4086 
4087 /***********************************************************************
4088  **
4089  ** Routines for 2D points.
4090  **
4091  ***********************************************************************/
4092 
4093 Datum
4095 {
4096  float8 x = PG_GETARG_FLOAT8(0);
4097  float8 y = PG_GETARG_FLOAT8(1);
4098  Point *result;
4099 
4100  result = (Point *) palloc(sizeof(Point));
4101 
4102  point_construct(result, x, y);
4103 
4104  PG_RETURN_POINT_P(result);
4105 }
4106 
4107 
4108 static inline void
4109 point_add_point(Point *result, Point *pt1, Point *pt2)
4110 {
4111  point_construct(result,
4112  float8_pl(pt1->x, pt2->x),
4113  float8_pl(pt1->y, pt2->y));
4114 }
4115 
4116 Datum
4118 {
4119  Point *p1 = PG_GETARG_POINT_P(0);
4120  Point *p2 = PG_GETARG_POINT_P(1);
4121  Point *result;
4122 
4123  result = (Point *) palloc(sizeof(Point));
4124 
4125  point_add_point(result, p1, p2);
4126 
4127  PG_RETURN_POINT_P(result);
4128 }
4129 
4130 
4131 static inline void
4132 point_sub_point(Point *result, Point *pt1, Point *pt2)
4133 {
4134  point_construct(result,
4135  float8_mi(pt1->x, pt2->x),
4136  float8_mi(pt1->y, pt2->y));
4137 }
4138 
4139 Datum
4141 {
4142  Point *p1 = PG_GETARG_POINT_P(0);
4143  Point *p2 = PG_GETARG_POINT_P(1);
4144  Point *result;
4145 
4146  result = (Point *) palloc(sizeof(Point));
4147 
4148  point_sub_point(result, p1, p2);
4149 
4150  PG_RETURN_POINT_P(result);
4151 }
4152 
4153 
4154 static inline void
4155 point_mul_point(Point *result, Point *pt1, Point *pt2)
4156 {
4157  point_construct(result,
4158  float8_mi(float8_mul(pt1->x, pt2->x),
4159  float8_mul(pt1->y, pt2->y)),
4160  float8_pl(float8_mul(pt1->x, pt2->y),
4161  float8_mul(pt1->y, pt2->x)));
4162 }
4163 
4164 Datum
4166 {
4167  Point *p1 = PG_GETARG_POINT_P(0);
4168  Point *p2 = PG_GETARG_POINT_P(1);
4169  Point *result;
4170 
4171  result = (Point *) palloc(sizeof(Point));
4172 
4173  point_mul_point(result, p1, p2);
4174 
4175  PG_RETURN_POINT_P(result);
4176 }
4177 
4178 
4179 static inline void
4180 point_div_point(Point *result, Point *pt1, Point *pt2)
4181 {
4182  float8 div;
4183 
4184  div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y));
4185 
4186  point_construct(result,
4187  float8_div(float8_pl(float8_mul(pt1->x, pt2->x),
4188  float8_mul(pt1->y, pt2->y)), div),
4189  float8_div(float8_mi(float8_mul(pt1->y, pt2->x),
4190  float8_mul(pt1->x, pt2->y)), div));
4191 }
4192 
4193 Datum
4195 {
4196  Point *p1 = PG_GETARG_POINT_P(0);
4197  Point *p2 = PG_GETARG_POINT_P(1);
4198  Point *result;
4199 
4200  result = (Point *) palloc(sizeof(Point));
4201 
4202  point_div_point(result, p1, p2);
4203 
4204  PG_RETURN_POINT_P(result);
4205 }
4206 
4207 
4208 /***********************************************************************
4209  **
4210  ** Routines for 2D boxes.
4211  **
4212  ***********************************************************************/
4213 
4214 Datum
4216 {
4217  Point *p1 = PG_GETARG_POINT_P(0);
4218  Point *p2 = PG_GETARG_POINT_P(1);
4219  BOX *result;
4220 
4221  result = (BOX *) palloc(sizeof(BOX));
4222 
4223  box_construct(result, p1, p2);
4224 
4225  PG_RETURN_BOX_P(result);
4226 }
4227 
4228 Datum
4230 {
4231  BOX *box = PG_GETARG_BOX_P(0);
4232  Point *p = PG_GETARG_POINT_P(1);
4233  BOX *result;
4234 
4235  result = (BOX *) palloc(sizeof(BOX));
4236 
4237  point_add_point(&result->high, &box->high, p);
4238  point_add_point(&result->low, &box->low, p);
4239 
4240  PG_RETURN_BOX_P(result);
4241 }
4242 
4243 Datum
4245 {
4246  BOX *box = PG_GETARG_BOX_P(0);
4247  Point *p = PG_GETARG_POINT_P(1);
4248  BOX *result;
4249 
4250  result = (BOX *) palloc(sizeof(BOX));
4251 
4252  point_sub_point(&result->high, &box->high, p);
4253  point_sub_point(&result->low, &box->low, p);
4254 
4255  PG_RETURN_BOX_P(result);
4256 }
4257 
4258 Datum
4260 {
4261  BOX *box = PG_GETARG_BOX_P(0);
4262  Point *p = PG_GETARG_POINT_P(1);
4263  BOX *result;
4264  Point high,
4265  low;
4266 
4267  result = (BOX *) palloc(sizeof(BOX));
4268 
4269  point_mul_point(&high, &box->high, p);
4270  point_mul_point(&low, &box->low, p);
4271 
4272  box_construct(result, &high, &low);
4273 
4274  PG_RETURN_BOX_P(result);
4275 }
4276 
4277 Datum
4279 {
4280  BOX *box = PG_GETARG_BOX_P(0);
4281  Point *p = PG_GETARG_POINT_P(1);
4282  BOX *result;
4283  Point high,
4284  low;
4285 
4286  result = (BOX *) palloc(sizeof(BOX));
4287 
4288  point_div_point(&high, &box->high, p);
4289  point_div_point(&low, &box->low, p);
4290 
4291  box_construct(result, &high, &low);
4292 
4293  PG_RETURN_BOX_P(result);
4294 }
4295 
4296 /*
4297  * Convert point to empty box
4298  */
4299 Datum
4301 {
4302  Point *pt = PG_GETARG_POINT_P(0);
4303  BOX *box;
4304 
4305  box = (BOX *) palloc(sizeof(BOX));
4306 
4307  box->high.x = pt->x;
4308  box->low.x = pt->x;
4309  box->high.y = pt->y;
4310  box->low.y = pt->y;
4311 
4312  PG_RETURN_BOX_P(box);
4313 }
4314 
4315 /*
4316  * Smallest bounding box that includes both of the given boxes
4317  */
4318 Datum
4320 {
4321  BOX *box1 = PG_GETARG_BOX_P(0),
4322  *box2 = PG_GETARG_BOX_P(1),
4323  *container;
4324 
4325  container = (BOX *) palloc(sizeof(BOX));
4326 
4327  container->high.x = float8_max(box1->high.x, box2->high.x);
4328  container->low.x = float8_min(box1->low.x, box2->low.x);
4329  container->high.y = float8_max(box1->high.y, box2->high.y);
4330  container->low.y = float8_min(box1->low.y, box2->low.y);
4331 
4332  PG_RETURN_BOX_P(container);
4333 }
4334 
4335 
4336 /***********************************************************************
4337  **
4338  ** Routines for 2D paths.
4339  **
4340  ***********************************************************************/
4341 
4342 /* path_add()
4343  * Concatenate two paths (only if they are both open).
4344  */
4345 Datum
4347 {
4348  PATH *p1 = PG_GETARG_PATH_P(0);
4349  PATH *p2 = PG_GETARG_PATH_P(1);
4350  PATH *result;
4351  int size,
4352  base_size;
4353  int i;
4354 
4355  if (p1->closed || p2->closed)
4356  PG_RETURN_NULL();
4357 
4358  base_size = sizeof(p1->p[0]) * (p1->npts + p2->npts);
4359  size = offsetof(PATH, p) + base_size;
4360 
4361  /* Check for integer overflow */
4362  if (base_size / sizeof(p1->p[0]) != (p1->npts + p2->npts) ||
4363  size <= base_size)
4364  ereport(ERROR,
4365  (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
4366  errmsg("too many points requested")));
4367 
4368  result = (PATH *) palloc(size);
4369 
4370  SET_VARSIZE(result, size);
4371  result->npts = (p1->npts + p2->npts);
4372  result->closed = p1->closed;
4373  /* prevent instability in unused pad bytes */
4374  result->dummy = 0;
4375 
4376  for (i = 0; i < p1->npts; i++)
4377  {
4378  result->p[i].x = p1->p[i].x;
4379  result->p[i].y = p1->p[i].y;
4380  }
4381  for (i = 0; i < p2->npts; i++)
4382  {
4383  result->p[i + p1->npts].x = p2->p[i].x;
4384  result->p[i + p1->npts].y = p2->p[i].y;
4385  }
4386 
4387  PG_RETURN_PATH_P(result);
4388 }
4389 
4390 /* path_add_pt()
4391  * Translation operators.
4392  */
4393 Datum
4395 {
4396  PATH *path = PG_GETARG_PATH_P_COPY(0);
4397  Point *point = PG_GETARG_POINT_P(1);
4398  int i;
4399 
4400  for (i = 0; i < path->npts; i++)
4401  point_add_point(&path->p[i], &path->p[i], point);
4402 
4403  PG_RETURN_PATH_P(path);
4404 }
4405 
4406 Datum
4408 {
4409  PATH *path = PG_GETARG_PATH_P_COPY(0);
4410  Point *point = PG_GETARG_POINT_P(1);
4411  int i;
4412 
4413  for (i = 0; i < path->npts; i++)
4414  point_sub_point(&path->p[i], &path->p[i], point);
4415 
4416  PG_RETURN_PATH_P(path);
4417 }
4418 
4419 /* path_mul_pt()
4420  * Rotation and scaling operators.
4421  */
4422 Datum
4424 {
4425  PATH *path = PG_GETARG_PATH_P_COPY(0);
4426  Point *point = PG_GETARG_POINT_P(1);
4427  int i;
4428 
4429  for (i = 0; i < path->npts; i++)
4430  point_mul_point(&path->p[i], &path->p[i], point);
4431 
4432  PG_RETURN_PATH_P(path);
4433 }
4434 
4435 Datum
4437 {
4438  PATH *path = PG_GETARG_PATH_P_COPY(0);
4439  Point *point = PG_GETARG_POINT_P(1);
4440  int i;
4441 
4442  for (i = 0; i < path->npts; i++)
4443  point_div_point(&path->p[i], &path->p[i], point);
4444 
4445  PG_RETURN_PATH_P(path);
4446 }
4447 
4448 
4449 Datum
4451 {
4452 #ifdef NOT_USED
4453  PATH *path = PG_GETARG_PATH_P(0);
4454 #endif
4455 
4456  ereport(ERROR,
4457  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
4458  errmsg("function \"path_center\" not implemented")));
4459 
4460  PG_RETURN_NULL();
4461 }
4462 
4463 Datum
4465 {
4466  PATH *path = PG_GETARG_PATH_P(0);
4467  POLYGON *poly;
4468  int size;
4469  int i;
4470 
4471  /* This is not very consistent --- other similar cases return NULL ... */
4472  if (!path->closed)
4473  ereport(ERROR,
4474  (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
4475  errmsg("open path cannot be converted to polygon")));
4476 
4477  /*
4478  * Never overflows: the old size fit in MaxAllocSize, and the new size is
4479  * just a small constant larger.
4480  */
4481  size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * path->npts;
4482  poly = (POLYGON *) palloc(size);
4483 
4484  SET_VARSIZE(poly, size);
4485  poly->npts = path->npts;
4486 
4487  for (i = 0; i < path->npts; i++)
4488  {
4489  poly->p[i].x = path->p[i].x;
4490  poly->p[i].y = path->p[i].y;
4491  }
4492 
4493  make_bound_box(poly);
4494 
4495  PG_RETURN_POLYGON_P(poly);
4496 }
4497 
4498 
4499 /***********************************************************************
4500  **
4501  ** Routines for 2D polygons.
4502  **
4503  ***********************************************************************/
4504 
4505 Datum
4507 {
4508  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4509 
4510  PG_RETURN_INT32(poly->npts);
4511 }
4512 
4513 
4514 Datum
4516 {
4517  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4518  Point *result;
4519  CIRCLE circle;
4520 
4521  result = (Point *) palloc(sizeof(Point));
4522 
4523  poly_to_circle(&circle, poly);
4524  *result = circle.center;
4525 
4526  PG_RETURN_POINT_P(result);
4527 }
4528 
4529 
4530 Datum
4532 {
4533  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4534  BOX *box;
4535 
4536  box = (BOX *) palloc(sizeof(BOX));
4537  *box = poly->boundbox;
4538 
4539  PG_RETURN_BOX_P(box);
4540 }
4541 
4542 
4543 /* box_poly()
4544  * Convert a box to a polygon.
4545  */
4546 Datum
4548 {
4549  BOX *box = PG_GETARG_BOX_P(0);
4550  POLYGON *poly;
4551  int size;
4552 
4553  /* map four corners of the box to a polygon */
4554  size = offsetof(POLYGON, p) + sizeof(poly->p[0]) * 4;
4555  poly = (POLYGON *) palloc(size);
4556 
4557  SET_VARSIZE(poly, size);
4558  poly->npts = 4;
4559 
4560  poly->p[0].x = box->low.x;
4561  poly->p[0].y = box->low.y;
4562  poly->p[1].x = box->low.x;
4563  poly->p[1].y = box->high.y;
4564  poly->p[2].x = box->high.x;
4565  poly->p[2].y = box->high.y;
4566  poly->p[3].x = box->high.x;
4567  poly->p[3].y = box->low.y;
4568 
4569  box_construct(&poly->boundbox, &box->high, &box->low);
4570 
4571  PG_RETURN_POLYGON_P(poly);
4572 }
4573 
4574 
4575 Datum
4577 {
4578  POLYGON *poly = PG_GETARG_POLYGON_P(0);
4579  PATH *path;
4580  int size;
4581  int i;
4582 
4583  /*
4584  * Never overflows: the old size fit in MaxAllocSize, and the new size is
4585  * smaller by a small constant.
4586  */
4587  size = offsetof(PATH, p) + sizeof(path->p[0]) * poly->npts;
4588  path = (PATH *) palloc(size);
4589 
4590  SET_VARSIZE(path, size);
4591  path->npts = poly->npts;
4592  path->closed = true;
4593  /* prevent instability in unused pad bytes */
4594  path->dummy = 0;
4595 
4596  for (i = 0; i < poly->npts; i++)
4597  {
4598  path->p[i].x = poly->p[i].x;
4599  path->p[i].y = poly->p[i].y;
4600  }
4601 
4602  PG_RETURN_PATH_P(path);
4603 }
4604 
4605 
4606 /***********************************************************************
4607  **
4608  ** Routines for circles.
4609  **
4610  ***********************************************************************/
4611 
4612 /*----------------------------------------------------------
4613  * Formatting and conversion routines.
4614  *---------------------------------------------------------*/
4615 
4616 /* circle_in - convert a string to internal form.
4617  *
4618  * External format: (center and radius of circle)
4619  * "<(f8,f8),f8>"
4620  * also supports quick entry style "f8,f8,f8"
4621  */
4622 Datum
4624 {
4625  char *str = PG_GETARG_CSTRING(0);
4626  CIRCLE *circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4627  char *s,
4628  *cp;
4629  int depth = 0;
4630 
4631  s = str;
4632  while (isspace((unsigned char) *s))
4633  s++;
4634  if (*s == LDELIM_C)
4635  depth++, s++;
4636  else if (*s == LDELIM)
4637  {
4638  /* If there are two left parens, consume the first one */
4639  cp = (s + 1);
4640  while (isspace((unsigned char) *cp))
4641  cp++;
4642  if (*cp == LDELIM)
4643  depth++, s = cp;
4644  }
4645 
4646  /* pair_decode will consume parens around the pair, if any */
4647  pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str);
4648 
4649  if (*s == DELIM)
4650  s++;
4651 
4652  circle->radius = single_decode(s, &s, "circle", str);
4653  /* We have to accept NaN. */
4654  if (circle->radius < 0.0)
4655  ereport(ERROR,
4656  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4657  errmsg("invalid input syntax for type %s: \"%s\"",
4658  "circle", str)));
4659 
4660  while (depth > 0)
4661  {
4662  if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1)))
4663  {
4664  depth--;
4665  s++;
4666  while (isspace((unsigned char) *s))
4667  s++;
4668  }
4669  else
4670  ereport(ERROR,
4671  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4672  errmsg("invalid input syntax for type %s: \"%s\"",
4673  "circle", str)));
4674  }
4675 
4676  if (*s != '\0')
4677  ereport(ERROR,
4678  (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
4679  errmsg("invalid input syntax for type %s: \"%s\"",
4680  "circle", str)));
4681 
4682  PG_RETURN_CIRCLE_P(circle);
4683 }
4684 
4685 /* circle_out - convert a circle to external form.
4686  */
4687 Datum
4689 {
4690  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4692 
4693  initStringInfo(&str);
4694 
4697  pair_encode(circle->center.x, circle->center.y, &str);
4699  appendStringInfoChar(&str, DELIM);
4700  single_encode(circle->radius, &str);
4702 
4703  PG_RETURN_CSTRING(str.data);
4704 }
4705 
4706 /*
4707  * circle_recv - converts external binary format to circle
4708  */
4709 Datum
4711 {
4713  CIRCLE *circle;
4714 
4715  circle = (CIRCLE *) palloc(sizeof(CIRCLE));
4716 
4717  circle->center.x = pq_getmsgfloat8(buf);
4718  circle->center.y = pq_getmsgfloat8(buf);
4719  circle->radius = pq_getmsgfloat8(buf);
4720 
4721  /* We have to accept NaN. */
4722  if (circle->radius < 0.0)
4723  ereport(ERROR,
4724  (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION),
4725  errmsg("invalid radius in external \"circle\" value")));
4726 
4727  PG_RETURN_CIRCLE_P(circle);
4728 }
4729 
4730 /*
4731  * circle_send - converts circle to binary format
4732  */
4733 Datum
4735 {
4736  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4738 
4739  pq_begintypsend(&buf);
4740  pq_sendfloat8(&buf, circle->center.x);
4741  pq_sendfloat8(&buf, circle->center.y);
4742  pq_sendfloat8(&buf, circle->radius);
4744 }
4745 
4746 
4747 /*----------------------------------------------------------
4748  * Relational operators for CIRCLEs.
4749  * <, >, <=, >=, and == are based on circle area.
4750  *---------------------------------------------------------*/
4751 
4752 /* circles identical?
4753  *
4754  * We consider NaNs values to be equal to each other to let those circles
4755  * to be found.
4756  */
4757 Datum
4759 {
4760  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4761  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4762 
4763  PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) ||
4764  FPeq(circle1->radius, circle2->radius)) &&
4765  point_eq_point(&circle1->center, &circle2->center));
4766 }
4767 
4768 /* circle_overlap - does circle1 overlap circle2?
4769  */
4770 Datum
4772 {
4773  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4774  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4775 
4776  PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4777  float8_pl(circle1->radius, circle2->radius)));
4778 }
4779 
4780 /* circle_overleft - is the right edge of circle1 at or left of
4781  * the right edge of circle2?
4782  */
4783 Datum
4785 {
4786  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4787  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4788 
4789  PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius),
4790  float8_pl(circle2->center.x, circle2->radius)));
4791 }
4792 
4793 /* circle_left - is circle1 strictly left of circle2?
4794  */
4795 Datum
4797 {
4798  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4799  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4800 
4801  PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius),
4802  float8_mi(circle2->center.x, circle2->radius)));
4803 }
4804 
4805 /* circle_right - is circle1 strictly right of circle2?
4806  */
4807 Datum
4809 {
4810  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4811  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4812 
4813  PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius),
4814  float8_pl(circle2->center.x, circle2->radius)));
4815 }
4816 
4817 /* circle_overright - is the left edge of circle1 at or right of
4818  * the left edge of circle2?
4819  */
4820 Datum
4822 {
4823  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4824  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4825 
4826  PG_RETURN_BOOL(FPge(float8_mi(circle1->center.x, circle1->radius),
4827  float8_mi(circle2->center.x, circle2->radius)));
4828 }
4829 
4830 /* circle_contained - is circle1 contained by circle2?
4831  */
4832 Datum
4834 {
4835  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4836  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4837 
4838  PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4839  float8_mi(circle2->radius, circle1->radius)));
4840 }
4841 
4842 /* circle_contain - does circle1 contain circle2?
4843  */
4844 Datum
4846 {
4847  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4848  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4849 
4850  PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center),
4851  float8_mi(circle1->radius, circle2->radius)));
4852 }
4853 
4854 
4855 /* circle_below - is circle1 strictly below circle2?
4856  */
4857 Datum
4859 {
4860  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4861  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4862 
4863  PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius),
4864  float8_mi(circle2->center.y, circle2->radius)));
4865 }
4866 
4867 /* circle_above - is circle1 strictly above circle2?
4868  */
4869 Datum
4871 {
4872  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4873  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4874 
4875  PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius),
4876  float8_pl(circle2->center.y, circle2->radius)));
4877 }
4878 
4879 /* circle_overbelow - is the upper edge of circle1 at or below
4880  * the upper edge of circle2?
4881  */
4882 Datum
4884 {
4885  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4886  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4887 
4888  PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius),
4889  float8_pl(circle2->center.y, circle2->radius)));
4890 }
4891 
4892 /* circle_overabove - is the lower edge of circle1 at or above
4893  * the lower edge of circle2?
4894  */
4895 Datum
4897 {
4898  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4899  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4900 
4901  PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius),
4902  float8_mi(circle2->center.y, circle2->radius)));
4903 }
4904 
4905 
4906 /* circle_relop - is area(circle1) relop area(circle2), within
4907  * our accuracy constraint?
4908  */
4909 Datum
4911 {
4912  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4913  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4914 
4915  PG_RETURN_BOOL(FPeq(circle_ar(circle1), circle_ar(circle2)));
4916 }
4917 
4918 Datum
4920 {
4921  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4922  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4923 
4924  PG_RETURN_BOOL(FPne(circle_ar(circle1), circle_ar(circle2)));
4925 }
4926 
4927 Datum
4929 {
4930  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4931  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4932 
4933  PG_RETURN_BOOL(FPlt(circle_ar(circle1), circle_ar(circle2)));
4934 }
4935 
4936 Datum
4938 {
4939  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4940  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4941 
4942  PG_RETURN_BOOL(FPgt(circle_ar(circle1), circle_ar(circle2)));
4943 }
4944 
4945 Datum
4947 {
4948  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4949  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4950 
4951  PG_RETURN_BOOL(FPle(circle_ar(circle1), circle_ar(circle2)));
4952 }
4953 
4954 Datum
4956 {
4957  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
4958  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
4959 
4960  PG_RETURN_BOOL(FPge(circle_ar(circle1), circle_ar(circle2)));
4961 }
4962 
4963 
4964 /*----------------------------------------------------------
4965  * "Arithmetic" operators on circles.
4966  *---------------------------------------------------------*/
4967 
4968 /* circle_add_pt()
4969  * Translation operator.
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_add_point(&result->center, &circle->center, point);
4981  result->radius = circle->radius;
4982 
4983  PG_RETURN_CIRCLE_P(result);
4984 }
4985 
4986 Datum
4988 {
4989  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
4990  Point *point = PG_GETARG_POINT_P(1);
4991  CIRCLE *result;
4992 
4993  result = (CIRCLE *) palloc(sizeof(CIRCLE));
4994 
4995  point_sub_point(&result->center, &circle->center, point);
4996  result->radius = circle->radius;
4997 
4998  PG_RETURN_CIRCLE_P(result);
4999 }
5000 
5001 
5002 /* circle_mul_pt()
5003  * Rotation and scaling operators.
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_mul_point(&result->center, &circle->center, point);
5015  result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y));
5016 
5017  PG_RETURN_CIRCLE_P(result);
5018 }
5019 
5020 Datum
5022 {
5023  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5024  Point *point = PG_GETARG_POINT_P(1);
5025  CIRCLE *result;
5026 
5027  result = (CIRCLE *) palloc(sizeof(CIRCLE));
5028 
5029  point_div_point(&result->center, &circle->center, point);
5030  result->radius = float8_div(circle->radius, HYPOT(point->x, point->y));
5031 
5032  PG_RETURN_CIRCLE_P(result);
5033 }
5034 
5035 
5036 /* circle_area - returns the area of the circle.
5037  */
5038 Datum
5040 {
5041  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5042 
5043  PG_RETURN_FLOAT8(circle_ar(circle));
5044 }
5045 
5046 
5047 /* circle_diameter - returns the diameter of the circle.
5048  */
5049 Datum
5051 {
5052  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5053 
5054  PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0));
5055 }
5056 
5057 
5058 /* circle_radius - returns the radius of the circle.
5059  */
5060 Datum
5062 {
5063  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5064 
5065  PG_RETURN_FLOAT8(circle->radius);
5066 }
5067 
5068 
5069 /* circle_distance - returns the distance between
5070  * two circles.
5071  */
5072 Datum
5074 {
5075  CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0);
5076  CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1);
5077  float8 result;
5078 
5079  result = float8_mi(point_dt(&circle1->center, &circle2->center),
5080  float8_pl(circle1->radius, circle2->radius));
5081  if (result < 0.0)
5082  result = 0.0;
5083 
5084  PG_RETURN_FLOAT8(result);
5085 }
5086 
5087 
5088 Datum
5090 {
5091  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5092  Point *point = PG_GETARG_POINT_P(1);
5093  float8 d;
5094 
5095  d = point_dt(&circle->center, point);
5096  PG_RETURN_BOOL(d <= circle->radius);
5097 }
5098 
5099 
5100 Datum
5102 {
5103  Point *point = PG_GETARG_POINT_P(0);
5104  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5105  float8 d;
5106 
5107  d = point_dt(&circle->center, point);
5108  PG_RETURN_BOOL(d <= circle->radius);
5109 }
5110 
5111 
5112 /* dist_pc - returns the distance between
5113  * a point and a circle.
5114  */
5115 Datum
5117 {
5118  Point *point = PG_GETARG_POINT_P(0);
5119  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5120  float8 result;
5121 
5122  result = float8_mi(point_dt(point, &circle->center),
5123  circle->radius);
5124  if (result < 0.0)
5125  result = 0.0;
5126 
5127  PG_RETURN_FLOAT8(result);
5128 }
5129 
5130 /*
5131  * Distance from a circle to a point
5132  */
5133 Datum
5135 {
5136  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5137  Point *point = PG_GETARG_POINT_P(1);
5138  float8 result;
5139 
5140  result = float8_mi(point_dt(point, &circle->center), circle->radius);
5141  if (result < 0.0)
5142  result = 0.0;
5143 
5144  PG_RETURN_FLOAT8(result);
5145 }
5146 
5147 /* circle_center - returns the center point of the circle.
5148  */
5149 Datum
5151 {
5152  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5153  Point *result;
5154 
5155  result = (Point *) palloc(sizeof(Point));
5156  result->x = circle->center.x;
5157  result->y = circle->center.y;
5158 
5159  PG_RETURN_POINT_P(result);
5160 }
5161 
5162 
5163 /* circle_ar - returns the area of the circle.
5164  */
5165 static float8
5167 {
5168  return float8_mul(float8_mul(circle->radius, circle->radius), M_PI);
5169 }
5170 
5171 
5172 /*----------------------------------------------------------
5173  * Conversion operators.
5174  *---------------------------------------------------------*/
5175 
5176 Datum
5178 {
5179  Point *center = PG_GETARG_POINT_P(0);
5180  float8 radius = PG_GETARG_FLOAT8(1);
5181  CIRCLE *result;
5182 
5183  result = (CIRCLE *) palloc(sizeof(CIRCLE));
5184 
5185  result->center.x = center->x;
5186  result->center.y = center->y;
5187  result->radius = radius;
5188 
5189  PG_RETURN_CIRCLE_P(result);
5190 }
5191 
5192 Datum
5194 {
5195  CIRCLE *circle = PG_GETARG_CIRCLE_P(0);
5196  BOX *box;
5197  float8 delta;
5198 
5199  box = (BOX *) palloc(sizeof(BOX));
5200 
5201  delta = float8_div(circle->radius, sqrt(2.0));
5202 
5203  box->high.x = float8_pl(circle->center.x, delta);
5204  box->low.x = float8_mi(circle->center.x, delta);
5205  box->high.y = float8_pl(circle->center.y, delta);
5206  box->low.y = float8_mi(circle->center.y, delta);
5207 
5208  PG_RETURN_BOX_P(box);
5209 }
5210 
5211 /* box_circle()
5212  * Convert a box to a circle.
5213  */
5214 Datum
5216 {
5217  BOX *box = PG_GETARG_BOX_P(0);
5218  CIRCLE *circle;
5219 
5220  circle = (CIRCLE *) palloc(sizeof(CIRCLE));
5221 
5222  circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0);
5223  circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0);
5224 
5225  circle->radius = point_dt(&circle->center, &box->high);
5226 
5227  PG_RETURN_CIRCLE_P(circle);
5228 }
5229 
5230 
5231 Datum
5233 {
5234  int32 npts = PG_GETARG_INT32(0);
5235  CIRCLE *circle = PG_GETARG_CIRCLE_P(1);
5236  POLYGON *poly;
5237  int base_size,
5238  size;
5239  int i;
5240  float8 angle;
5241  float8 anglestep;
5242 
5243  if (FPzero(circle->radius))
5244  ereport(ERROR,
5245  (errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
5246  errmsg("cannot convert circle with radius zero to polygon")));
5247 
5248  if (npts < 2)
5249  ereport(ERROR,
5250  (errcode(ERRCODE_INVALID_PARAMETER_VALUE),
5251  errmsg("must request at least 2 points")));
5252 
5253  base_size = sizeof(poly->p[0]) * npts;
5254  size = offsetof(POLYGON, p) + base_size;
5255 
5256  /* Check for integer overflow */
5257  if (base_size / npts != sizeof(poly->p[0]) || size <= base_size)
5258  ereport(ERROR,
5259  (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED),
5260  errmsg("too many points requested")));
5261 
5262  poly = (POLYGON *) palloc0(size); /* zero any holes */
5263  SET_VARSIZE(poly, size);
5264  poly->npts = npts;
5265 
5266  anglestep = float8_div(2.0 * M_PI, npts);
5267 
5268  for (i = 0; i < npts; i++)
5269  {
5270  angle = float8_mul(anglestep, i);
5271 
5272  poly->p[i].x = float8_mi(circle->center.x,
5273  float8_mul(circle->radius, cos(angle)));
5274  poly->p[i].y = float8_pl(circle->center.y,
5275  float8_mul(circle->radius, sin(angle)));
5276  }
5277 
5278  make_bound_box(poly);
5279 
5280  PG_RETURN_POLYGON_P(poly);
5281 }
5282 
5283 /*
5284  * Convert polygon to circle
5285  *
5286  * The result must be preallocated.
5287  *
5288  * XXX This algorithm should use weighted means of line segments
5289  * rather than straight average values of points - tgl 97/01/21.
5290  */
5291 static void
5293 {
5294  int i;
5295 
5296  Assert(poly->npts > 0);
5297 
5298  result->center.x = 0;
5299  result->center.y = 0;
5300  result->radius = 0;
5301 
5302  for (i = 0; i < poly->npts; i++)
5303  point_add_point(&result->center, &result->center, &poly->p[i]);
5304  result->center.x = float8_div(result->center.x, poly->npts);
5305  result->center.y = float8_div(result->center.y, poly->npts);
5306 
5307  for (i = 0; i < poly->npts; i++)
5308  result->radius = float8_pl(result->radius,
5309  point_dt(&poly->p[i], &result->center));
5310  result->radius = float8_div(result->radius, poly->npts);
5311 }
5312 
5313 Datum
5315 {
5316  POLYGON *poly = PG_GETARG_POLYGON_P(0);
5317  CIRCLE *result;
5318 
5319  result = (CIRCLE *) palloc(sizeof(CIRCLE));
5320 
5321  poly_to_circle(result, poly);
5322 
5323  PG_RETURN_CIRCLE_P(result);
5324 }
5325 
5326 
5327 /***********************************************************************
5328  **
5329  ** Private routines for multiple types.
5330  **
5331  ***********************************************************************/
5332 
5333 /*
5334  * Test to see if the point is inside the polygon, returns 1/0, or 2 if
5335  * the point is on the polygon.
5336  * Code adapted but not copied from integer-based routines in WN: A
5337  * Server for the HTTP
5338  * version 1.15.1, file wn/image.c
5339  * http://hopf.math.northwestern.edu/index.html
5340  * Description of algorithm: http://www.linuxjournal.com/article/2197
5341  * http://www.linuxjournal.com/article/2029
5342  */
5343 
5344 #define POINT_ON_POLYGON INT_MAX
5345 
5346 static int
5347 point_inside(Point *p, int npts, Point *plist)
5348 {
5349  float8 x0,
5350  y0;
5351  float8 prev_x,
5352  prev_y;
5353  int i = 0;
5354  float8 x,
5355  y;
5356  int cross,
5357  total_cross = 0;
5358 
5359  Assert(npts > 0);
5360 
5361  /* compute first polygon point relative to single point */
5362  x0 = float8_mi(plist[0].x, p->x);
5363  y0 = float8_mi(plist[0].y, p->y);
5364 
5365  prev_x = x0;
5366  prev_y = y0;
5367  /* loop over polygon points and aggregate total_cross */
5368  for (i = 1; i < npts; i++)
5369  {
5370  /* compute next polygon point relative to single point */
5371  x = float8_mi(plist[i].x, p->x);
5372  y = float8_mi(plist[i].y, p->y);
5373 
5374  /* compute previous to current point crossing */
5375  if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON)
5376  return 2;
5377  total_cross += cross;
5378 
5379  prev_x = x;
5380  prev_y = y;
5381  }
5382 
5383  /* now do the first point */
5384  if ((cross = lseg_crossing(x0, y0, prev_x, prev_y)) == POINT_ON_POLYGON)
5385  return 2;
5386  total_cross += cross;
5387 
5388  if (total_cross != 0)
5389  return 1;
5390  return 0;
5391 }
5392 
5393 
5394 /* lseg_crossing()
5395  * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction.
5396  * Returns +/-1 if one point is on the positive X-axis.
5397  * Returns 0 if both points are on the positive X-axis, or there is no crossing.
5398  * Returns POINT_ON_POLYGON if the segment contains (0,0).
5399  * Wow, that is one confusing API, but it is used above, and when summed,
5400  * can tell is if a point is in a polygon.
5401  */
5402 
5403 static int
5404 lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y)
5405 {
5406  float8 z;
5407  int y_sign;
5408 
5409  if (FPzero(y))
5410  { /* y == 0, on X axis */
5411  if (FPzero(x)) /* (x,y) is (0,0)? */
5412  return POINT_ON_POLYGON;
5413  else if (FPgt(x, 0))
5414  { /* x > 0 */
5415  if (FPzero(prev_y)) /* y and prev_y are zero */
5416  /* prev_x > 0? */
5417  return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
5418  return FPlt(prev_y, 0.0) ? 1 : -1;
5419  }
5420  else
5421  { /* x < 0, x not on positive X axis */
5422  if (FPzero(prev_y))
5423  /* prev_x < 0? */
5424  return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON;
5425  return 0;
5426  }
5427  }
5428  else
5429  { /* y != 0 */
5430  /* compute y crossing direction from previous point */
5431  y_sign = FPgt(y, 0.0) ? 1 : -1;
5432 
5433  if (FPzero(prev_y))
5434  /* previous point was on X axis, so new point is either off or on */
5435  return FPlt(prev_x, 0.0) ? 0 : y_sign;
5436  else if ((y_sign < 0 && FPlt(prev_y, 0.0)) ||
5437  (y_sign > 0 && FPgt(prev_y, 0.0)))
5438  /* both above or below X axis */
5439  return 0; /* same sign */
5440  else
5441  { /* y and prev_y cross X-axis */
5442  if (FPge(x, 0.0) && FPgt(prev_x, 0.0))
5443  /* both non-negative so cross positive X-axis */
5444  return 2 * y_sign;
5445  if (FPlt(x, 0.0) && FPle(prev_x, 0.0))
5446  /* both non-positive so do not cross positive X-axis */
5447  return 0;
5448 
5449  /* x and y cross axes, see URL above point_inside() */
5450  z = float8_mi(float8_mul(float8_mi(x, prev_x), y),
5451  float8_mul(float8_mi(y, prev_y), x));
5452  if (FPzero(z))
5453  return POINT_ON_POLYGON;
5454  if ((y_sign < 0 && FPlt(z, 0.0)) ||
5455  (y_sign > 0 && FPgt(z, 0.0)))
5456  return 0;
5457  return 2 * y_sign;
5458  }
5459  }
5460 }
5461 
5462 
5463 static bool
5464 plist_same(int npts, Point *p1, Point *p2)
5465 {
5466  int i,
5467  ii,
5468  j;
5469 
5470  /* find match for first point */
5471  for (i = 0; i < npts; i++)
5472  {
5473  if (point_eq_point(&p2[i], &p1[0]))
5474  {
5475 
5476  /* match found? then look forward through remaining points */
5477  for (ii = 1, j = i + 1; ii < npts; ii++, j++)
5478  {
5479  if (j >= npts)
5480  j = 0;
5481  if (!point_eq_point(&p2[j], &p1[ii]))
5482  break;
5483  }
5484  if (ii == npts)
5485  return true;
5486 
5487  /* match not found forwards? then look backwards */
5488  for (ii = 1, j = i - 1; ii < npts; ii++, j--)
5489  {
5490  if (j < 0)
5491  j = (npts - 1);
5492  if (!point_eq_point(&p2[j], &p1[ii]))
5493  break;
5494  }
5495  if (ii == npts)
5496  return true;
5497  }
5498  }
5499 
5500  return false;
5501 }
5502 
5503 
5504 /*-------------------------------------------------------------------------
5505  * Determine the hypotenuse.
5506  *
5507  * If required, x and y are swapped to make x the larger number. The
5508  * traditional formula of x^2+y^2 is rearranged to factor x outside the
5509  * sqrt. This allows computation of the hypotenuse for significantly
5510  * larger values, and with a higher precision than when using the naive
5511  * formula. In particular, this cannot overflow unless the final result
5512  * would be out-of-range.
5513  *
5514  * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
5515  * = x * sqrt( 1 + y^2/x^2 )
5516  * = x * sqrt( 1 + y/x * y/x )
5517  *
5518  * It is expected that this routine will eventually be replaced with the
5519  * C99 hypot() function.
5520  *
5521  * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
5522  * case of hypot(inf,nan) results in INF, and not NAN.
5523  *-----------------------------------------------------------------------
5524  */
5525 float8
5527 {
5528  float8 yx,
5529  result;
5530 
5531  /* Handle INF and NaN properly */
5532  if (isinf(x) || isinf(y))
5533  return get_float8_infinity();
5534 
5535  if (isnan(x) || isnan(y))
5536  return get_float8_nan();
5537 
5538  /* Else, drop any minus signs */
5539  x = fabs(x);
5540  y = fabs(y);
5541 
5542  /* Swap x and y if needed to make x the larger one */
5543  if (x < y)
5544  {
5545  float8 temp = x;
5546 
5547  x = y;
5548  y = temp;
5549  }
5550 
5551  /*
5552  * If y is zero, the hypotenuse is x. This test saves a few cycles in
5553  * such cases, but more importantly it also protects against
5554  * divide-by-zero errors, since now x >= y.
5555  */
5556  if (y == 0.0)
5557  return x;
5558 
5559  /* Determine the hypotenuse */
5560  yx = y / x;
5561  result = x * sqrt(1.0 + (yx * yx));
5562 
5563  if (unlikely(isinf(result)))
5565  if (unlikely(result == 0.0))
5567 
5568  return result;
5569 }
Datum boxes