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