Context Navigation

geo_ops.c

Revision 47, 114.4 kB (checked in by cs186, 4 years ago)
Make MyStuff? consistent with the source tree!

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

Note: See TracBrowser for help on using the browser.

Context Navigation

root / hw / hw2 / MyStuff / geo_ops.c

Download in other formats: