geo_ops.c source code [PostgreSQL/src/backend/utils/adt/geo_ops.c]

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

Browse the source code of PostgreSQL/src/backend/utils/adt/geo_ops.c