Context Navigation

Back to Ticket #63

Ticket #63: max_distance.patch

File max_distance.patch, 56.3 KB (added by nicklas, 3 years ago)

liblwgeom/liblwgeom.h


9	9	*
10	10	* This is free software; you can redistribute and/or modify it under
11	11	* the terms of the GNU General Public Licence. See the COPYING file.
12		*
	12	*
13	13	**********************************************************************/
14	14
15	15	#ifndef _LIBLWGEOM_H
…	…
21	21
22	22	/**
23	23	* @file liblwgeom.h
24		*
	24	*
25	25	* This library is the generic geometry handling section of PostGIS. The geometry
26	26	* objects, constructors, destructors, and a set of spatial processing functions,
27		* are implemented here.
28		*
29		* The library is designed for use in non-PostGIS applications if necessary. The
30		* units tests at cunit/cu_tester.c and the loader/dumper programs at
	27	* are implemented here.
	28	*
	29	* The library is designed for use in non-PostGIS applications if necessary. The
	30	* units tests at cunit/cu_tester.c and the loader/dumper programs at
31	31	* ../loader/shp2pgsql.c are examples of non-PostGIS applications using liblwgeom.
32		*
	32	*
33	33	* Programs using this library should set up the default memory managers and error
34		* handlers by implementing an lwgeom_init_allocators() function, which can be as
	34	* handlers by implementing an lwgeom_init_allocators() function, which can be as
35	35	* a wrapper around the lwgeom_install_default_allocators() function if you want
36	36	* no special handling for memory management and error reporting.
37	37	*/
…	…
71	71
72	72	#endif
73	73
74		/**
75		* Global functions for memory/logging handlers.
	74	/**
	75	* Global functions for memory/logging handlers.
76	76	*/
77	77	typedef void* (*lwallocator)(size_t size);
78	78	typedef void* (lwreallocator)(void mem, size_t size);
…	…
113	113	void lwerror(const char *fmt, ...);
114	114
115	115	/**
116		* The default memory/logging handlers installed by lwgeom_install_default_allocators()
	116	* The default memory/logging handlers installed by lwgeom_install_default_allocators()
117	117	*/
118	118	void *default_allocator(size_t size);
119	119	void default_reallocator(void mem, size_t size);
…	…
167	167
168	168	typedef struct
169	169	{
170		double xmin, ymin, zmin;
171		double xmax, ymax, zmax;
	170	double xmin, ymin, zmin;
	171	double xmax, ymax, zmax;
172	172	} BOX3D;
173	173
174	174	typedef struct chiptag
…	…
241	241	* NOTE: for GEOS integration, we'll probably set z=NaN
242	242	* so look out - z might be NaN for 2d geometries!
243	243	*/
244		typedef struct { double x,y,z; } POINT3DZ;
245		typedef struct { double x,y,z; } POINT3D; /* alias for POINT3DZ */
246		typedef struct { double x,y,m; } POINT3DM;
	244	typedef struct
	245	{
	246	double x,y,z;
	247	} POINT3DZ;
	248	typedef struct
	249	{
	250	double x,y,z;
	251	} POINT3D; /* alias for POINT3DZ */
	252	typedef struct
	253	{
	254	double x,y,m;
	255	} POINT3DM;
247	256
248	257
249	258	/*
…	…
252	261	*/
253	262	typedef struct
254	263	{
255		double x;
256		double y;
	264	double x;
	265	double y;
257	266	} POINT2D;
258	267
259	268	typedef struct
260	269	{
261		double x;
262		double y;
263		double z;
264		double m;
	270	double x;
	271	double y;
	272	double z;
	273	double m;
265	274	} POINT4D;
266	275
	276	typedef struct
	277	{
	278	double d; /the distance between p1 and p2/
	279	POINT2D p1;
	280	POINT2D p2;
	281	int thedir; /the direction of looking, if thedir = -1 then we look for maxdistance and if it is 1 then we look for mindistance/
	282	int twisted; /To preserve the order of incoming points to match the first and secon point in shortest and longest line/
	283	} DISTPTS;
267	284	/******************************************************************/
268	285
269	286	/*
…	…
274	291	*/
275	292	typedef struct
276	293	{
277		/* array of POINT 2D, 3D or 4D. probably missaligned. */
278		uchar *serialized_pointlist;
	294	/* array of POINT 2D, 3D or 4D. probably missaligned. */
	295	uchar *serialized_pointlist;
279	296
280		/* use TYPE_* macros to handle */
281		uchar dims;
	297	/* use TYPE_* macros to handle */
	298	uchar dims;
282	299
283		uint32 npoints;
	300	uint32 npoints;
284	301	} POINTARRAY;
285	302
286	303
287	304	/*
288	305	* Use the following to build pointarrays
289		* when number of points in output is not
	306	* when number of points in output is not
290	307	* known in advance
291	308	*/
292		typedef struct {
	309	typedef struct
	310	{
293	311	POINTARRAY *pa;
294	312	size_t ptsize;
295	313	size_t capacity; /* given in points */
…	…
311	329	* respectively.
312	330	*/
313	331	extern int dynptarray_addPoint4d(DYNPTARRAY dpa, POINT4D p4d,
314		int allow_duplicates);
	332	int allow_duplicates);
315	333
316	334	/******************************************************************
317	335	*
…	…
321	339
322	340	typedef struct
323	341	{
324		uchar type;
	342	uchar type;
325	343	BOX2DFLOAT4 *bbox;
326	344	uint32 SRID; /* -1 == unneeded */
327	345	void *data;
…	…
332	350	{
333	351	uchar type; /* POINTTYPE */
334	352	BOX2DFLOAT4 *bbox;
335		uint32 SRID;
336		POINTARRAY point; / hide 2d/3d (this will be an array of 1 point) */
	353	uint32 SRID;
	354	POINTARRAY point; / hide 2d/3d (this will be an array of 1 point) */
337	355	} LWPOINT; /* "light-weight point" */
338	356
339	357	/* LINETYPE */
…	…
341	359	{
342	360	uchar type; /* LINETYPE */
343	361	BOX2DFLOAT4 *bbox;
344		uint32 SRID;
	362	uint32 SRID;
345	363	POINTARRAY points; / array of POINT3D */
346	364	} LWLINE; /* "light-weight line" */
347	365
…	…
350	368	{
351	369	uchar type; /* POLYGONTYPE */
352	370	BOX2DFLOAT4 *bbox;
353		uint32 SRID;
	371	uint32 SRID;
354	372	int nrings;
355	373	POINTARRAY *rings; / list of rings (list of points) */
356	374	} LWPOLY; /* "light-weight polygon" */
…	…
358	376	/* MULTIPOINTTYPE */
359	377	typedef struct
360	378	{
361		uchar type;
	379	uchar type;
362	380	BOX2DFLOAT4 *bbox;
363		uint32 SRID;
	381	uint32 SRID;
364	382	int ngeoms;
365	383	LWPOINT **geoms;
366		} LWMPOINT;
	384	} LWMPOINT;
367	385
368	386	/* MULTILINETYPE */
369	387	typedef struct
370		{
371		uchar type;
	388	{
	389	uchar type;
372	390	BOX2DFLOAT4 *bbox;
373		uint32 SRID;
	391	uint32 SRID;
374	392	int ngeoms;
375	393	LWLINE **geoms;
376		} LWMLINE;
	394	} LWMLINE;
377	395
378	396	/* MULTIPOLYGONTYPE */
379	397	typedef struct
380		{
381		uchar type;
	398	{
	399	uchar type;
382	400	BOX2DFLOAT4 *bbox;
383		uint32 SRID;
	401	uint32 SRID;
384	402	int ngeoms;
385	403	LWPOLY **geoms;
386		} LWMPOLY;
	404	} LWMPOLY;
387	405
388	406	/* COLLECTIONTYPE */
389	407	typedef struct
390		{
391		uchar type;
	408	{
	409	uchar type;
392	410	BOX2DFLOAT4 *bbox;
393		uint32 SRID;
	411	uint32 SRID;
394	412	int ngeoms;
395	413	LWGEOM **geoms;
396		} LWCOLLECTION;
	414	} LWCOLLECTION;
397	415
398	416	/* CIRCSTRINGTYPE */
399	417	typedef struct
400	418	{
401		uchar type; /* CIRCSTRINGTYPE */
402		BOX2DFLOAT4 *bbox;
403		uint32 SRID;
404		POINTARRAY points; / array of POINT(3D/3DM) */
	419	uchar type; /* CIRCSTRINGTYPE */
	420	BOX2DFLOAT4 *bbox;
	421	uint32 SRID;
	422	POINTARRAY points; / array of POINT(3D/3DM) */
405	423	} LWCIRCSTRING; /* "light-weight circularstring" */
406	424
407	425	/* COMPOUNDTYPE */
408	426	typedef struct
409	427	{
410		uchar type; /* COMPOUNDTYPE */
411		BOX2DFLOAT4 *bbox;
412		uint32 SRID;
413		int ngeoms;
414		LWGEOM **geoms;
	428	uchar type; /* COMPOUNDTYPE */
	429	BOX2DFLOAT4 *bbox;
	430	uint32 SRID;
	431	int ngeoms;
	432	LWGEOM **geoms;
415	433	} LWCOMPOUND; /* "light-weight compound line" */
416	434
417	435	/* CURVEPOLYTYPE */
418	436	typedef struct
419	437	{
420		uchar type; /* CURVEPOLYTYPE */
421		BOX2DFLOAT4 *bbox;
422		uint32 SRID;
423		int nrings;
424		LWGEOM *rings; / list of rings (list of points) */
	438	uchar type; /* CURVEPOLYTYPE */
	439	BOX2DFLOAT4 *bbox;
	440	uint32 SRID;
	441	int nrings;
	442	LWGEOM *rings; / list of rings (list of points) */
425	443	} LWCURVEPOLY; /* "light-weight polygon" */
426	444
427	445	/* MULTICURVE */
428	446	typedef struct
429	447	{
430		uchar type;
431		BOX2DFLOAT4 *bbox;
432		uint32 SRID;
433		int ngeoms;
434		LWGEOM **geoms;
	448	uchar type;
	449	BOX2DFLOAT4 *bbox;
	450	uint32 SRID;
	451	int ngeoms;
	452	LWGEOM **geoms;
435	453	} LWMCURVE;
436	454
437	455	/* MULTISURFACETYPE */
438	456	typedef struct
439	457	{
440		uchar type;
441		BOX2DFLOAT4 *bbox;
442		uint32 SRID;
443		int ngeoms;
444		LWGEOM **geoms;
	458	uchar type;
	459	BOX2DFLOAT4 *bbox;
	460	uint32 SRID;
	461	int ngeoms;
	462	LWGEOM **geoms;
445	463	} LWMSURFACE;
446	464
447	465	/* Casts LWGEOM->LW* (return NULL if cast is illegal) */
…	…
466	484	extern LWGEOM lwpoint_as_lwgeom(LWPOINT obj);
467	485
468	486	/*
469		* Call this function everytime LWGEOM coordinates
	487	* Call this function everytime LWGEOM coordinates
470	488	* change so to invalidate bounding box
471	489	*/
472	490	extern void lwgeom_changed(LWGEOM *lwgeom);
…	…
552	570	* possibly corrupt the POINTARRAY.
553	571	*
554	572	* WARNING: Don't cast this to a POINT !
555		* it would not be reliable due to memory alignment constraints
	573	* it would not be reliable due to memory alignment constraints
556	574	*/
557	575	extern uchar getPoint_internal(const POINTARRAY pa, int n);
558	576
…	…
568	586	* 'points' points to. No data conversion is done.
569	587	*/
570	588	extern POINTARRAY pointArray_construct(uchar points, char hasz, char hasm,
571		uint32 npoints);
	589	uint32 npoints);
572	590
573		/*
	591	/*
574	592	* Calculate the (BOX3D) bounding box of a set of points.
575	593	* Returns an alloced BOX3D or NULL (for empty geom) in the first form.
576	594	* Write result in user-provided BOX3D in second form (return 0 if untouched).
…	…
639	657	* This is the binary representation of lwgeom compatible
640	658	* with postgresql varlena struct
641	659	*/
642		typedef struct {
	660	typedef struct
	661	{
643	662	uint32 size; /* varlena header (do not touch directly!) */
644	663	uchar type; /* encodes ndims, type, bbox presence,
645	664	srid presence */
…	…
656	675	* from the serialized form.
657	676	*/
658	677	extern PG_LWGEOM PG_LWGEOM_construct(uchar serialized, int SRID,
659		int wantbbox);
	678	int wantbbox);
660	679
661	680	/*
662	681	* Compute bbox of serialized geom
…	…
715	734
716	735	/*
717	736	* convert this point into its serialize form
718		* result's first char will be the 8bit type.
	737	* result's first char will be the 8bit type.
719	738	* See serialized form doc
720	739	*/
721	740	extern uchar lwpoint_serialize(LWPOINT point);
…	…
724	743	extern void lwpoint_serialize_buf(LWPOINT point, uchar buf, size_t *size);
725	744
726	745	/*
727		* find bounding box (standard one)
	746	* find bounding box (standard one)
728	747	* zmin=zmax=0 if 2d (might change to NaN)
729	748	*/
730	749	extern BOX3D lwpoint_compute_box3d(LWPOINT point);
…	…
1060	1079	* MEMORY MANAGEMENT
1061	1080	****************************************************************/
1062	1081
1063		/*
1064		* The lwfree_* family of functions frees all memory associated
	1082	/*
	1083	* The lwfree_* family of functions frees all memory associated
1065	1084	* with the pointer, including the serialized__pointlist in the
1066	1085	* point arrays. Do not use these on LWGEOMs de-serialized from
1067	1086	* PG_LWGEOMs or they will try to free an underlying structure
1068		* managed by PgSQL. Only use these on LWGEOMs you have
	1087	* managed by PgSQL. Only use these on LWGEOMs you have
1069	1088	* constructed yourself.
1070	1089	*/
1071
	1090
1072	1091	extern void ptarray_free(POINTARRAY *pa);
1073	1092	extern void lwpoint_free(LWPOINT *pt);
1074	1093	extern void lwline_free(LWLINE *line);
…	…
1084	1103
1085	1104	/*
1086	1105	* The *_release family of functions frees the LWGEOM structures
1087		* surrounding the POINTARRAYs but leaves the POINTARRAYs
	1106	* surrounding the POINTARRAYs but leaves the POINTARRAYs
1088	1107	* intact. Use these on LWGEOMs that have been de-serialized
1089		* from PG_LWGEOMs. Do not use these on LWGEOMs you have
	1108	* from PG_LWGEOMs. Do not use these on LWGEOMs you have
1090	1109	* constructed yourself, or you will leak lots of memory.
1091	1110	*/
1092	1111
…	…
1145	1164	extern void lwgeom_force3dz_recursive(uchar serialized, uchar optr, size_t *retsize);
1146	1165	extern void lwgeom_force3dm_recursive(uchar serialized, uchar optr, size_t *retsize);
1147	1166	extern void lwgeom_force4d_recursive(uchar serialized, uchar optr, size_t *retsize);
	1167
1148	1168	extern double distance2d_pt_pt(POINT2D p1, POINT2D p2);
	1169	extern void lw_dist2d_comp_pt_pt(POINT2D p1, POINT2D p2, DISTPTS *dl);
1149	1170	extern double distance2d_pt_seg(POINT2D p, POINT2D A, POINT2D *B);
1150		extern double distance2d_seg_seg(POINT2D A, POINT2D B, POINT2D C, POINT2D D);
1151		extern double distance2d_pt_ptarray(POINT2D p, POINTARRAY pa);
1152		extern double distance2d_ptarray_ptarray(POINTARRAY l1, POINTARRAY l2);
	1171	extern void lw_dist2d_comp_pt_seg(POINT2D p, POINT2D A, POINT2D B, DISTPTS dl);
	1172	extern void lw_dist2d_comp_seg_seg(POINT2D A, POINT2D B, POINT2D C, POINT2D D, DISTPTS *dl);
	1173	extern void lw_dist2d_comp_pt_ptarray(POINT2D p, POINTARRAY pa, DISTPTS *dl);
	1174	extern void lw_dist2d_comp_ptarray_ptarray(POINTARRAY l1, POINTARRAY l2, DISTPTS *dl);
1153	1175	extern int pt_in_ring_2d(POINT2D p, POINTARRAY ring);
1154	1176	extern int pt_in_poly_2d(POINT2D p, LWPOLY poly);
1155		extern ~~double distance2d_ptarray_poly(POINTARRAY pa, LWPOLY poly~~);
1156		extern ~~double distance2d_point_point(LWPOINT point1, LWPOINT point2~~);
1157		extern ~~double distance2d_point_line(LWPOINT point, LWLINE line~~);
1158		extern ~~double distance2d_line_line(LWLINE line1, LWLINE line2~~);
1159		extern ~~double distance2d_point_poly(LWPOINT point, LWPOLY poly~~);
1160		extern ~~double distance2d_poly_poly(LWPOLY poly1, LWPOLY poly2~~);
1161		extern ~~double distance2d_line_poly(LWLINE line, LWPOLY poly~~);
	1177	extern void lw_dist2d_comp_ptarray_poly(POINTARRAY pa, LWPOLY poly, DISTPTS *dl);
	1178	extern void lw_dist2d_comp_point_point(LWPOINT point1, LWPOINT point2, DISTPTS *dl);
	1179	extern void lw_dist2d_comp_point_line(LWPOINT point, LWLINE line, DISTPTS *dl);
	1180	extern void lw_dist2d_comp_line_line(LWLINE line1, LWLINE line2, DISTPTS *dl);
	1181	extern void lw_dist2d_comp_point_poly(LWPOINT point, LWPOLY poly, DISTPTS *dl);
	1182	extern void lw_dist2d_comp_poly_poly(LWPOLY poly1, LWPOLY poly2, DISTPTS *dl);
	1183	extern void lw_dist2d_comp_line_poly(LWLINE line, LWPOLY poly, DISTPTS *dl);
1162	1184	extern int azimuth_pt_pt(POINT2D p1, POINT2D p2, double *ret);
1163	1185	extern double lwgeom_mindistance2d_recursive(uchar lw1, uchar lw2);
1164	1186	extern double lwgeom_mindistance2d_recursive_tolerance(uchar lw1, uchar lw2, double tolerance);
	1187	extern double lwgeom_maxdistance2d_recursive(uchar lw1, uchar lw2);
	1188	extern double lwgeom_maxdistance2d_recursive_tolerance(uchar lw1, uchar lw2, double tolerance);
	1189	extern void lw_dist2d_comp(uchar lw1, uchar lw2, double tolerance, DISTPTS *dl);
1165	1190	extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad);
1166	1191	extern int32 lwgeom_npoints(uchar *serialized);
1167	1192	extern char ptarray_isccw(const POINTARRAY *pa);
…	…
1224	1249	/*
1225	1250	* Clone an LWGEOM
1226	1251	* pointarray are not copied.
1227		* BBOXes are copied
	1252	* BBOXes are copied
1228	1253	*/
1229	1254	extern LWGEOM lwgeom_clone(const LWGEOM lwgeom);
1230	1255	extern LWPOINT lwpoint_clone(const LWPOINT lwgeom);
…	…
1240	1265	* Take ownership of arguments
1241	1266	*/
1242	1267	extern LWPOINT lwpoint_construct(int SRID, BOX2DFLOAT4 bbox,
1243		POINTARRAY *point);
	1268	POINTARRAY *point);
1244	1269	extern LWLINE lwline_construct(int SRID, BOX2DFLOAT4 bbox,
1245		POINTARRAY *points);
	1270	POINTARRAY *points);
1246	1271
1247	1272	/*
1248	1273	* Construct a new LWPOLY. arrays (points/points per ring) will NOT be copied
1249	1274	* use SRID=-1 for unknown SRID (will have 8bit type's S = 0)
1250	1275	*/
1251	1276	extern LWPOLY lwpoly_construct(int SRID, BOX2DFLOAT4 bbox,
1252		unsigned int nrings, POINTARRAY **points);
	1277	unsigned int nrings, POINTARRAY **points);
1253	1278
1254	1279	extern LWCOLLECTION *lwcollection_construct(unsigned int type, int SRID,
1255		BOX2DFLOAT4 bbox, unsigned int ngeoms, LWGEOM *geoms);
	1280	BOX2DFLOAT4 bbox, unsigned int ngeoms, LWGEOM *geoms);
1256	1281	extern LWCOLLECTION *lwcollection_construct_empty(int SRID,
1257		char hasZ, char hasM);
	1282	char hasZ, char hasM);
1258	1283
1259	1284	/*
1260	1285	* Construct a new LWCIRCSTRING. arrays (points/points per ring) will NOT be copied
…	…
1288	1313	*/
1289	1314
1290	1315	extern POINTARRAY ptarray_addPoint(POINTARRAY pa, uchar *p, size_t pdims,
1291		unsigned int where);
	1316	unsigned int where);
1292	1317	extern POINTARRAY ptarray_removePoint(POINTARRAY pa, unsigned int where);
1293	1318
1294	1319	extern int ptarray_isclosed2d(const POINTARRAY *pa);
…	…
1323	1348	#define PARSER_CHECK_ALL (PARSER_CHECK_MINPOINTS \| PARSER_CHECK_ODD \| PARSER_CHECK_CLOSURE)
1324	1349
1325	1350	/*
1326		* Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM
	1351	* Parser result structure: returns the result of attempting to convert (E)WKT/(E)WKB to LWGEOM
1327	1352	*/
1328	1353	typedef struct struct_lwgeom_parser_result
1329	1354	{
…	…
1338	1363	* Parser error messages (these must match the message array in lwgparse.c)
1339	1364	*/
1340	1365	#define PARSER_ERROR_MOREPOINTS 1
1341		#define PARSER_ERROR_ODDPOINTS 2
1342		#define PARSER_ERROR_UNCLOSED 3
1343		#define PARSER_ERROR_MIXDIMS 4
	1366	#define PARSER_ERROR_ODDPOINTS 2
	1367	#define PARSER_ERROR_UNCLOSED 3
	1368	#define PARSER_ERROR_MIXDIMS 4
1344	1369	#define PARSER_ERROR_INVALIDGEOM 5
1345	1370	#define PARSER_ERROR_INVALIDWKBTYPE 6
1346	1371	#define PARSER_ERROR_INCONTINUOUS 7
1347	1372
1348	1373
1349	1374	/*
1350		* Unparser result structure: returns the result of attempting to convert LWGEOM to (E)WKT/(E)WKB
	1375	* Unparser result structure: returns the result of attempting to convert LWGEOM to (E)WKT/(E)WKB
1351	1376	*/
1352	1377	typedef struct struct_lwgeom_unparser_result
1353	1378	{
…	…
1362	1387	* Unparser error messages (these must match the message array in lwgunparse.c)
1363	1388	*/
1364	1389	#define UNPARSER_ERROR_MOREPOINTS 1
1365		#define UNPARSER_ERROR_ODDPOINTS 2
1366		#define UNPARSER_ERROR_UNCLOSED 3
	1390	#define UNPARSER_ERROR_ODDPOINTS 2
	1391	#define UNPARSER_ERROR_UNCLOSED 3
1367	1392
1368	1393
1369	1394	/* Parser access routines */

liblwgeom/measures.c


7	7	*
8	8	* This is free software; you can redistribute and/or modify it under
9	9	* the terms of the GNU General Public Licence. See the COPYING file.
10		*
	10	*
11	11	**********************************************************************/
12	12
13	13	#include <math.h>
…	…
39	39	if ( memcmp(&first, &last, sizeof(POINT2D)) )
40	40	{
41	41	lwerror("pt_in_ring_2d: V[n] != V[0] (%g %g != %g %g)",
42		first.x, first.y, last.x, last.y);
43
	42	first.x, first.y, last.x, last.y);
	43
44	44	}
45	45	#endif
46	46
…	…
49	49
50	50	/* loop through all edges of the polygon */
51	51	getPoint2d_p(ring, 0, &v1);
52		for (i=0; i<ring->npoints-1; i++)
53		{
	52	for (i=0; i<ring->npoints-1; i++)
	53	{
54	54	double vt;
55	55	getPoint2d_p(ring, i+1, &v2);
56	56
57	57	/* edge from vertex i to vertex i+1 */
58		if
	58	if
59	59	(
60		/* an upward crossing */
61		((v1.y <= p->y) && (v2.y > p->y))
62		/* a downward crossing */
63		\|\| ((v1.y > p->y) && (v2.y <= p->y))
	60	/* an upward crossing */
	61	((v1.y <= p->y) && (v2.y > p->y))
	62	/* a downward crossing */
	63	\|\| ((v1.y > p->y) && (v2.y <= p->y))
64	64	)
65		{
	65	{
66	66
67	67	vt = (double)(p->y - v1.y) / (v2.y - v1.y);
68	68
…	…
70	70	if (p->x < v1.x + vt * (v2.x - v1.x))
71	71	{
72	72	/* a valid crossing of y=p.y right of p.x */
73		++cn;
	73	++cn;
74	74	}
75	75	}
76	76	v1 = v2;
…	…
81	81	return (cn&1); /* 0 if even (out), and 1 if odd (in) */
82	82	}
83	83
	84	/Compares incomming points and stores the points closest to each other or most far away from each other depending on dl->thedir (thedirection) /
	85	void lw_dist2d_comp_pt_pt(POINT2D thep1, POINT2D thep2,DISTPTS *dl)
	86	{
	87	double hside = thep2->x - thep1->x;
	88	double vside = thep2->y - thep1->y;
	89	double dist = sqrt ( hsidehside + vsidevside );
	90	if (((dl->d - dist)*(dl->thedir))>0)
	91	{
	92	dl->d = dist;
	93	/* lwnotice("lw_dist2d_comp_pt_pt %e",dl->twisted);*/
	94	if (dl->twisted>0)
	95	{
	96	dl->p1 = *thep1;
	97	dl->p2 = *thep2;
	98	}
	99	else
	100	{
	101	dl->p1 = *thep2;
	102	dl->p2 = *thep1;
	103	}
	104	}
	105	return;
	106	}
	107
	108	/The old function nessecary for ptarray_segmentize2d in ptarray.c/
84	109	double distance2d_pt_pt(POINT2D p1, POINT2D p2)
85	110	{
86	111	double hside = p2->x - p1->x;
…	…
94	119	); */
95	120	}
96	121
97		/distance2d from p to line A->B /
	122	/*lw_dist2d_comp from p to line A->B
	123	This one is now sending every occation to lw_dist2d_comp_pt_pt
	124	Berfore it was handling occations where r was between 0 and 1 internally and just returning the distance without identifying the points.
	125	To get this points it was nessecary to change and it also showed to be about 10%faster. */
	126	void lw_dist2d_comp_pt_seg(POINT2D p, POINT2D A, POINT2D B, DISTPTS dl)
	127	{
	128	/lwnotice("lw_dist2d_comp_pt_seg %e",dl->d);/
	129	double r;
	130
	131	/if start==end, then use pt distance /
	132	if ( ( A->x == B->x) && (A->y == B->y) )
	133	{
	134	lw_dist2d_comp_pt_pt(p,A,dl);
	135	return;
	136	}
	137	/*
	138	* otherwise, we use comp.graphics.algorithms
	139	* Frequently Asked Questions method
	140	*
	141	* (1) AC dot AB
	142	* r = ---------
	143	* \|\|AB\|\|^2
	144	* r has the following meaning:
	145	* r=0 P = A
	146	* r=1 P = B
	147	* r<0 P is on the backward extension of AB
	148	* r>1 P is on the forward extension of AB
	149	* 0<r<1 P is interior to AB
	150	*/
	151
	152	r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( (B->x-A->x)(B->x-A->x) +(B->y-A->y)(B->y-A->y) );
	153
	154	/*This is for finding the maxdistance.
	155	the maxdistance have to be between two vertexes as I understand it,
	156	compared to mindistance which can be ebtween some point between the vertices and a vertex.*/
	157	if (dl->thedir < 1)
	158	{
	159	if (r>=0.5)
	160	{
	161	lw_dist2d_comp_pt_pt(p,A,dl);
	162	return;
	163	}
	164	if (r<0.5)
	165	{
	166	lw_dist2d_comp_pt_pt(p,B,dl);
	167	return;
	168	}
	169	}
	170
	171	if (r<0) /If the first vertex A is closest to the point p/
	172	{
	173	lw_dist2d_comp_pt_pt(p,A,dl);
	174	return;
	175	}
	176	if (r>1) /If the second vertex B is closest to the point p/
	177	{
	178
	179	lw_dist2d_comp_pt_pt(p,B,dl);
	180	return;
	181	}
	182
	183
	184	/else if the point p is closer to some point between a and b then we find that point and send it to lw_dist2d_comp_pt_pt/
	185	POINT2D c;
	186	c.x=A->x + r * (B->x-A->x);
	187	c.y=A->y + r * (B->y-A->y);
	188
	189	lw_dist2d_comp_pt_pt(p,&c,dl);
	190	return;
	191	}
	192
	193
	194
	195
	196	/The old function nessecary for ptarray_segmentize2d in ptarray.c/
98	197	double distance2d_pt_seg(POINT2D p, POINT2D A, POINT2D *B)
99	198	{
100	199	double r,s;
…	…
108	207	* Frequently Asked Questions method
109	208	*
110	209	* (1) AC dot AB
111		* r = ---------
112		* \|\|AB\|\|^2
	210	* r = ---------
	211	* \|\|AB\|\|^2
113	212	* r has the following meaning:
114	213	* r=0 P = A
115	214	* r=1 P = B
…	…
135	234	*/
136	235
137	236	s = ( (A->y-p->y)(B->x-A->x)- (A->x-p->x)(B->y-A->y) ) /
138		( (B->x-A->x)(B->x-A->x) +(B->y-A->y)(B->y-A->y) );
	237	( (B->x-A->x)(B->x-A->x) +(B->y-A->y)(B->y-A->y) );
139	238
140	239	return LW_ABS(s) * sqrt(
141		(B->x-A->x)(B->x-A->x) + (B->y-A->y)(B->y-A->y)
142		);
	240	(B->x-A->x)(B->x-A->x) + (B->y-A->y)(B->y-A->y)
	241	);
143	242	}
144	243
145		/* find the minimum 2d distance from AB to CD */
146		double distance2d_seg_seg(POINT2D A, POINT2D B, POINT2D C, POINT2D D)
	244
	245
	246
	247
	248	/*This function is changed so it is not doing any comparasion of distance
	249	but just sending every possible combination further to lw_dist2d_comp_pt_seg*/
	250	void lw_dist2d_comp_seg_seg(POINT2D A, POINT2D B, POINT2D C, POINT2D D, DISTPTS *dl)
147	251	{
148
	252	/lwnotice("lw_dist2d_comp_seg_seg %e \n",dl->d);/
149	253	double s_top, s_bot,s;
150	254	double r_top, r_bot,r;
151	255
152		~~LWDEBUGF(2, "distance2d_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",~~
153		~~A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);~~
154	256
	257	LWDEBUGF(2, "lw_dist2d_comp_seg_seg [%g,%g]->[%g,%g] by [%g,%g]->[%g,%g]",
	258	A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
155	259
156	260	/A and B are the same point /
157	261	if ( ( A->x == B->x) && (A->y == B->y) )
158		return distance2d_pt_seg(A,C,D);
	262	{
	263	lw_dist2d_comp_pt_seg(A,C,D,dl);
	264	return;
	265	}
	266	/U and V are the same point /
159	267
160		~~/U and V are the same point /~~
161
162	268	if ( ( C->x == D->x) && (C->y == D->y) )
163		return distance2d_pt_seg(D,A,B);
164
	269	{
	270	dl->twisted= ((dl->twisted) * (-1));
	271	lw_dist2d_comp_pt_seg(D,A,B,dl);
	272	return;
	273	}
165	274	/* AB and CD are line segments */
166	275	/* from comp.graphics.algo
167	276
…	…
190	299	s_top = (A->y-C->y)(B->x-A->x) - (A->x-C->x)(B->y-A->y);
191	300	s_bot = (B->x-A->x)(D->y-C->y) - (B->y-A->y)(D->x-C->x);
192	301
193		~~if ( (r_bot==0) \|\| (s_bot == 0) )~~
194		{
195		~~return (~~
196		~~LW_MIN(distance2d_pt_seg(A,C,D),~~
197		~~LW_MIN(distance2d_pt_seg(B,C,D),~~
198		~~LW_MIN(distance2d_pt_seg(C,A,B),~~
199		~~distance2d_pt_seg(D,A,B))~~
200		)
201		)
202		);
203		}
204	302	s = s_top/s_bot;
205	303	r= r_top/r_bot;
206	304
207		if ((r<0) \|\| (r>1) \|\| (s<0) \|\| (s>1) )
	305
	306	if ((r<0)\|\|(r>1)\|\|(s<0)\|\|(s>1)\|\|(dl->thedir <0))
208	307	{
209		/no intersection /
210		return (
211		LW_MIN(distance2d_pt_seg(A,C,D),
212		LW_MIN(distance2d_pt_seg(B,C,D),
213		LW_MIN(distance2d_pt_seg(C,A,B),
214		distance2d_pt_seg(D,A,B))
215		)
216		)
217		);
218
	308	lw_dist2d_comp_pt_seg(A,C,D,dl);
	309	lw_dist2d_comp_pt_seg(B,C,D,dl);
	310
	311	dl->twisted= ((dl->twisted) * (-1)); /here we change the order of inputted geometrys and that we notice by changing sign on dl->twisted/
	312	lw_dist2d_comp_pt_seg(C,A,B,dl);
	313	lw_dist2d_comp_pt_seg(D,A,B,dl);
	314	return;
219	315	}
220	316	else
221		return -0; /intersection exists /
	317	{
	318	if (dl->thedir >0) /If there is intersection we identify the intersection point and return it but only if we are looking for mindistance/
	319	{ /lwnotice("intersection lines: %e %e %e %e",A->x,A->y,B->x,B->y);/
	320	POINT2D theP;
	321
	322	if(((A->x==C->x)&&(A->y==C->y))\|\|((A->x==D->x)&&(A->y==D->y)))
	323	{
	324	theP.x = A->x;
	325	theP.y = A->y;
	326	}
	327	else if(((B->x==C->x)&&(B->y==C->y))\|\|((B->x==D->x)&&(B->y==D->y)))
	328	{
	329	theP.x = B->x;
	330	theP.y = B->y;
	331	}
	332	else
	333	{
	334	theP.x = A->x+r*(B->x-A->x);
	335	theP.y = A->y+r*(B->y-A->y);
	336	}
	337	/lwnotice("intersection: %e %e ",theP.x, theP.y);/
	338	dl->d=0;
	339	dl->p1=theP;
	340	dl->p2=theP;
	341	}
	342	return;
222	343
	344	}
	345
223	346	}
224	347
225		/*
226		* search all the segments of pointarray to see which one is closest to p1
227		* Returns minimum distance between point and pointarray
228		*/
229		double distance2d_pt_ptarray(POINT2D p, POINTARRAY pa)
	348	void lw_dist2d_comp_pt_ptarray(POINT2D p, POINTARRAY pa,DISTPTS *dl)
230	349	{
231		double result = 0;
	350
	351
232	352	int t;
233	353	POINT2D start, end;
234
	354	int twist = dl->twisted;
235	355	getPoint2d_p(pa, 0, &start);
236	356
237	357	for (t=1; t<pa->npoints; t++)
238	358	{
239		d~~ouble d~~ist;
	359	dl->twisted=twist;
240	360	getPoint2d_p(pa, t, &end);
241		dist = distance2d_pt_seg(p, &start, &end);
242		if (t==1) result = dist;
243		else result = LW_MIN(result, dist);
	361	lw_dist2d_comp_pt_seg(p, &start, &end,dl);
	362	start = end;
244	363
245		~~if ( result == 0 ) return 0;~~
246
247		~~start = end;~~
248	364	}
249	365
250		return ~~result~~;
	366	return;
251	367	}
252	368
253	369	/* test each segment of l1 against each segment of l2. Return min */
254		~~double distance2d_ptarray_ptarray(POINTARRAY l1, POINTARRAY l2~~)
	370	void lw_dist2d_comp_ptarray_ptarray(POINTARRAY l1, POINTARRAY l2,DISTPTS *dl)
255	371	{
256		double result = 99999999999.9;
257		char result_okay = 0; /result is a valid min /
	372	/lwnotice("lw_dist2d_comp_ptarray_ptarray");/
258	373	int t,u;
259	374	POINT2D start, end;
260	375	POINT2D start2, end2;
	376	int twist = dl->twisted;
	377	LWDEBUGF(2, "lw_dist2d_comp_ptarray_ptarray called (points: %d-%d)",
	378	l1->npoints, l2->npoints);
261	379
262		~~LWDEBUGF(2, "distance2d_ptarray_ptarray called (points: %d-%d)",~~
263		~~l1->npoints, l2->npoints);~~
264
265	380	getPoint2d_p(l1, 0, &start);
266	381	for (t=1; t<l1->npoints; t++) /for each segment in L1 /
267	382	{
268	383	getPoint2d_p(l1, t, &end);
269
270	384	getPoint2d_p(l2, 0, &start2);
271	385	for (u=1; u<l2->npoints; u++) /for each segment in L2 /
272	386	{
273		~~double dist;~~
274
275	387	getPoint2d_p(l2, u, &end2);
	388	dl->twisted=twist;
	389	lw_dist2d_comp_seg_seg(&start, &end, &start2, &end2,dl);
	390	#if PGIS_DEBUG > 1
	391	printf("line_line; seg %i * seg %i, dist = %g\n",t,u,dist_this);
	392	#endif
276	393
277		dist = distance2d_seg_seg(&start, &end, &start2, &end2);
	394	#ifdef PGIS_DEBUG
	395	lwnotice(" seg%d-seg%d dist: %f, mindist: %f",
	396	t, u, dist, result);
	397	#endif
278	398
279		~~LWDEBUGF(4, "line_line; seg %i * seg %i, dist = %g\n",t,u,dist);~~
280
281		~~if (result_okay)~~
282		~~result = LW_MIN(result,dist);~~
283		~~else~~
284		{
285		~~result_okay = 1;~~
286		~~result = dist;~~
287		}
288
289		~~LWDEBUGF(3, " seg%d-seg%d dist: %f, mindist: %f",~~
290		~~t, u, dist, result);~~
291
292		~~if (result <= 0) return 0; /intersection /~~
293
294	399	start2 = end2;
295	400	}
296	401	start = end;
297	402	}
298
299		return result;
	403	return ;
300	404	}
301	405
	406
	407
	408
	409
302	410	/* true if point is in poly (and not in its holes) */
303	411	int pt_in_poly_2d(POINT2D p, LWPOLY poly)
304	412	{
…	…
327	435	* otherwise return min distance to a ring (could be outside
328	436	* polygon or inside a hole)
329	437	*/
330		double distance2d_ptarray_poly(POINTARRAY pa, LWPOLY poly)
	438
	439
	440
	441
	442
	443
	444	void lw_dist2d_comp_ptarray_poly(POINTARRAY pa, LWPOLY poly, DISTPTS *dl)
331	445	{
332	446	POINT2D pt;
333	447	int i;
334		~~double mindist = 0;~~
335	448
336		~~LWDEBUGF(2, "distance2d_ptarray_poly called (%d rings)", poly->nrings);~~
337	449
	450	LWDEBUGF(2, "lw_dist2d_comp_ptarray_poly called (%d rings)", poly->nrings);
	451
	452
338	453	for (i=0; i<poly->nrings; i++)
339	454	{
340		double dist = distance2d_ptarray_ptarray(pa, poly->rings[i]);
341		if (i) mindist = LW_MIN(mindist, dist);
342		else mindist = dist;
	455	lw_dist2d_comp_ptarray_ptarray(pa, poly->rings[i], dl);
343	456
344	457	LWDEBUGF(3, " distance from ring %d: %f, mindist: %f",
345		i, dist, mindist);
	458	i, dist, mindist);
346	459
347		~~if ( mindist <= 0 ) return 0.0; /* intersection */~~
348	460	}
349	461
350	462	/*
…	…
357	469	* Outside outer ring, so min distance to a ring
358	470	* is the actual min distance
359	471	*/
360		if ( ! pt_in_ring_2d(&pt, poly->rings[0]) ) return mindist;
	472	if ( ! pt_in_ring_2d(&pt, poly->rings[0]) )
	473	{
	474	return ;
	475	}
361	476
362
363	477	/*
364	478	* Its in the outer ring.
365	479	* Have to check if its inside a hole
…	…
372	486	* Its inside a hole, then the actual
373	487	* distance is the min ring distance
374	488	*/
375		return ~~mindist~~;
	489	return;
376	490	}
377	491	}
	492	if (dl->thedir >0)
	493	{
	494	dl->d=0;
	495	dl->p1.x=pt.x;
	496	dl->p1.y=pt.y;
	497	dl->p2.x=pt.x;
	498	dl->p2.y=pt.y;
	499	}
	500	return ; /* Not in hole, so inside polygon */
378	501
379		~~return 0.0; /* Not in hole, so inside polygon */~~
380	502	}
381	503
382		double distance2d_point_point(LWPOINT point1, LWPOINT point2)
	504
	505
	506
	507
	508	void lw_dist2d_comp_point_point(LWPOINT point1, LWPOINT point2, DISTPTS *dl)
383	509	{
384	510	POINT2D p1;
385	511	POINT2D p2;
…	…
387	513	getPoint2d_p(point1->point, 0, &p1);
388	514	getPoint2d_p(point2->point, 0, &p2);
389	515
390		return distance2d_pt_pt(&p1, &p2);
	516	lw_dist2d_comp_pt_pt(&p1, &p2,dl);
	517	return;
391	518	}
392	519
393		~~double distance2d_point_line(LWPOINT point, LWLINE line~~)
	520	void lw_dist2d_comp_point_line(LWPOINT point, LWLINE line, DISTPTS *dl)
394	521	{
395	522	POINT2D p;
396	523	POINTARRAY *pa = line->points;
397	524	getPoint2d_p(point->point, 0, &p);
398		return distance2d_pt_ptarray(&p, pa);
	525	lw_dist2d_comp_pt_ptarray(&p, pa, dl);
	526	return;
399	527	}
400	528
401		~~double distance2d_line_line(LWLINE line1, LWLINE line2~~)
	529	void lw_dist2d_comp_line_line(LWLINE line1, LWLINE line2, DISTPTS *dl)
402	530	{
403	531	POINTARRAY *pa1 = line1->points;
404	532	POINTARRAY *pa2 = line2->points;
405		return distance2d_ptarray_ptarray(pa1, pa2);
	533	lw_dist2d_comp_ptarray_ptarray(pa1, pa2, dl);
	534	return;
406	535	}
407	536
408	537	/*
…	…
410	539	* 2. if in the boundary, test to see if its in a hole.
411	540	* if so, then return dist to hole, else return 0 (point in polygon)
412	541	*/
413		double distance2d_point_poly(LWPOINT point, LWPOLY poly)
	542
	543
	544
	545
	546
	547	void lw_dist2d_comp_point_poly(LWPOINT point, LWPOLY poly, DISTPTS *dl)
414	548	{
415	549	POINT2D p;
416	550	int i;
417	551
418	552	getPoint2d_p(point->point, 0, &p);
419	553
420		LWDEBUG(2, "~~distance2d~~_point_poly called");
	554	LWDEBUG(2, "lw_dist2d_comp_point_poly called");
421	555
	556
422	557	/* Return distance to outer ring if not inside it */
423	558	if ( ! pt_in_ring_2d(&p, poly->rings[0]) )
424	559	{
	560
425	561	LWDEBUG(3, " not inside outer-ring");
426	562
427		return distance2d_pt_ptarray(&p, poly->rings[0]);
	563
	564	lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl);
	565	return;
428	566	}
429	567
	568	if (dl->thedir < 0)
	569	{
	570	LWDEBUG(3, "inside outer-ring but looking for maxdistance");
	571	lw_dist2d_comp_pt_ptarray(&p, poly->rings[0], dl);
	572	return;
	573	}
430	574	/*
431	575	* Inside the outer ring.
432	576	* Scan though each of the inner rings looking to
433	577	* see if its inside. If not, distance==0.
434	578	* Otherwise, distance = pt to ring distance
435	579	*/
436		for (i=1; i<poly->nrings; i++)
	580	for (i=1; i<poly->nrings; i++)
437	581	{
438	582	/* Inside a hole. Distance = pt -> ring */
439	583	if ( pt_in_ring_2d(&p, poly->rings[i]) )
440	584	{
441	585	LWDEBUG(3, " inside an hole");
442	586
443		return distance2d_pt_ptarray(&p, poly->rings[i]);
	587	lw_dist2d_comp_pt_ptarray(&p, poly->rings[i], dl);
	588	return;
444	589	}
445	590	}
446	591
447	592	LWDEBUG(3, " inside the polygon");
448
449		return 0.0; /* Is inside the polygon */
	593	if (dl->thedir >0)
	594	{
	595	dl->d=0;
	596	dl->p1.x=p.x;
	597	dl->p1.y=p.y;
	598	dl->p2.x=p.x;
	599	dl->p2.y=p.y;
	600	}
	601	return ; /* Is inside the polygon */
450	602	}
451	603
	604
	605
	606
452	607	/*
453	608	* Brute force.
454	609	* Test to see if any rings intersect.
…	…
456	611	* Test to see if one inside the other and if they are inside holes.
457	612	* Find min distance ring-to-ring.
458	613	*/
459		double distance2d_poly_poly(LWPOLY poly1, LWPOLY poly2)
	614
	615
	616
	617	void lw_dist2d_comp_poly_poly(LWPOLY poly1, LWPOLY poly2, DISTPTS *dl)
460	618	{
461	619	POINT2D pt;
462		double mindist = -1;
	620
463	621	int i;
464	622
465		LWDEBUG(2, "~~distance2d~~_poly_poly called");
	623	LWDEBUG(2, "lw_dist2d_comp_poly_poly called");
466	624
467	625	/* if poly1 inside poly2 return 0 */
468	626	getPoint2d_p(poly1->rings[0], 0, &pt);
469		if ( pt_in_poly_2d(&pt, poly2) ) return 0.0;
	627	if (( pt_in_poly_2d(&pt, poly2)) && (dl->thedir >0))
	628	{
	629	dl->d=0;
	630	dl->p1.x=pt.x;
	631	dl->p1.y=pt.y;
	632	dl->p2.x=pt.x;
	633	dl->p2.y=pt.y;
	634	return ;
	635	}
470	636
471	637	/* if poly2 inside poly1 return 0 */
472	638	getPoint2d_p(poly2->rings[0], 0, &pt);
473		if ( pt_in_poly_2d(&pt, poly1) ) return 0.0;
	639	if ( (pt_in_poly_2d(&pt, poly1)) && (dl->thedir >0))
	640	{
	641	dl->d=0;
	642	dl->p1.x=pt.x;
	643	dl->p1.y=pt.y;
	644	dl->p2.x=pt.x;
	645	dl->p2.y=pt.y;
	646	return ;
	647	}
474	648
475	649	LWDEBUG(3, " polys not inside each other");
476	650
…	…
484	658	int j;
485	659	for (j=0; j<poly2->nrings; j++)
486	660	{
487		double d = distance2d_ptarray_ptarray(poly1->rings[i],
488		poly2->rings[j]);
489		if ( d <= 0 ) return 0.0;
490
	661	lw_dist2d_comp_ptarray_ptarray(poly1->rings[i], poly2->rings[j],dl);
	662	if ( dl->d <= 0 )
	663	{
	664	if (dl->thedir >0)
	665	{
	666	return;
	667	}
	668
	669	}
491	670	/* mindist is -1 when not yet set */
492		if (mindist > -1) mindist = LW_MIN(mindist, d);
493		else mindist = d;
494
495		LWDEBUGF(3, " ring%i-%i dist: %f, mindist: %f", i, j, d, mindist);
	671	LWDEBUGF(3, " ring%i-%i dist: %f", i, j, d);
496	672	}
497
498	673	}
499
500	674	/* otherwise return closest approach of rings (no intersection) */
501		return ~~mindist~~;
	675	return ;
502	676
503	677	}
504	678
505		~~double distance2d_line_poly(LWLINE line, LWPOLY poly~~)
	679	void lw_dist2d_comp_line_poly(LWLINE line, LWPOLY poly, DISTPTS *dl)
506	680	{
507		return distance2d_ptarray_poly(line->points, poly);
	681	lw_dist2d_comp_ptarray_poly(line->points, poly, dl);
	682	return;
508	683	}
509	684
510
511	685	/find the 2d length of the given POINTARRAY (even if it's 3d) /
512	686	double lwgeom_pointarray_length2d(POINTARRAY *pts)
513	687	{
…	…
517	691	POINT2D to;
518	692
519	693	if ( pts->npoints < 2 ) return 0.0;
520		for (i=0; i<pts->npoints-1;i++)
	694	for (i=0; i<pts->npoints-1; i++)
521	695	{
522	696	getPoint2d_p(pts, i, &frm);
523	697	getPoint2d_p(pts, i+1, &to);
524	698	dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) ) +
525		((frm.y - to.y)*(frm.y - to.y) ) );
	699	((frm.y - to.y)*(frm.y - to.y) ) );
526	700	}
527	701	return dist;
528	702	}
…	…
544	718	/* compute 2d length if 3d is not available */
545	719	if ( ! TYPE_HASZ(pts->dims) ) return lwgeom_pointarray_length2d(pts);
546	720
547		for (i=0; i<pts->npoints-1;i++)
	721	for (i=0; i<pts->npoints-1; i++)
548	722	{
549	723	getPoint3dz_p(pts, i, &frm);
550	724	getPoint3dz_p(pts, i+1, &to);
551	725	dist += sqrt( ( (frm.x - to.x)*(frm.x - to.x) ) +
552		((frm.y - to.y)*(frm.y - to.y) ) +
553		((frm.z - to.z)*(frm.z - to.z) ) );
	726	((frm.y - to.y)*(frm.y - to.y) ) +
	727	((frm.z - to.z)*(frm.z - to.z) ) );
554	728	}
555
556	729	return dist;
557	730	}
558	731
…	…
562	735	double
563	736	lwgeom_curvepolygon_area(LWCURVEPOLY *curvepoly)
564	737	{
565		LWPOLY poly = (LWPOLY )lwgeom_segmentize((LWGEOM *)curvepoly, 32);
566		return lwgeom_polygon_area(poly);
	738	LWPOLY poly = (LWPOLY )lwgeom_segmentize((LWGEOM *)curvepoly, 32);
	739	return lwgeom_polygon_area(poly);
567	740	}
568	741
569	742	/*
…	…
589	762	LWDEBUGF(4, " rings %d has %d points", i, ring->npoints);
590	763
591	764	for (j=0; j<ring->npoints-1; j++)
592		{
	765	{
593	766	getPoint2d_p(ring, j, &p1);
594	767	getPoint2d_p(ring, j+1, &p2);
595	768	ringarea += ( p1.x * p2.y ) - ( p1.y * p2.x );
…	…
644	817	}
645	818
646	819	double
	820	lwgeom_maxdistance2d_recursive(uchar lw1, uchar lw2)
	821	{
	822	return lwgeom_maxdistance2d_recursive_tolerance( lw1, lw2, 99999999 );
	823	}
	824
	825
	826	double
	827	lwgeom_maxdistance2d_recursive_tolerance(uchar lw1, uchar lw2, double tolerance)
	828	{
	829	/double thedist;/
	830	DISTPTS thedl;
	831	thedl.thedir = (-1);
	832	thedl.d= thedl.thedir * 99999999;
	833	DISTPTS *dl = &thedl;
	834	lw_dist2d_comp( lw1, lw2,tolerance, dl );
	835	return dl->d;
	836	/lwgeom_release((DISTPTS )dl);
	837	return thedist;
	838	return 1212.00;*/
	839	}
	840
	841	double
647	842	lwgeom_mindistance2d_recursive(uchar lw1, uchar lw2)
648	843	{
649		return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 );
	844	return lwgeom_mindistance2d_recursive_tolerance( lw1, lw2, 0.0 );
650	845	}
651	846
	847
652	848	double
653	849	lwgeom_mindistance2d_recursive_tolerance(uchar lw1, uchar lw2, double tolerance)
654	850	{
	851	double thedist;
	852	DISTPTS thedl;
	853	thedl.thedir = (1);
	854	thedl.d= thedl.thedir * 99999999;
	855	DISTPTS *dl = &thedl;
	856	lw_dist2d_comp( lw1, lw2,tolerance, dl );
	857	thedist = dl->d;
	858	/lwgeom_release((DISTPTS )dl);*/
	859	return thedist;
	860	}
	861
	862	void
	863	lw_dist2d_comp(uchar lw1, uchar lw2, double tolerance, DISTPTS *dl)
	864	{
655	865	LWGEOM_INSPECTED in1, in2;
656	866	int i, j;
657		~~double mindist = -1;~~
658	867
659	868	in1 = lwgeom_inspect(lw1);
660	869	in2 = lwgeom_inspect(lw2);
…	…
663	872	{
664	873	uchar *g1 = lwgeom_getsubgeometry_inspected(in1, i);
665	874	int t1 = lwgeom_getType(g1[0]);
666		~~double dist=tolerance;~~
667	875
668		/* ~~Argument 1 i~~s a multitype... recurse */
	876	/* it's a multitype... recurse */
669	877	if ( lwgeom_contains_subgeoms(t1) )
670	878	{
671		dist = lwgeom_mindistance2d_recursive_tolerance(g1, lw2, tolerance);
672		if ( dist <= tolerance ) return tolerance; /* can't be closer */
673		if ( mindist == -1 ) mindist = dist;
674		else mindist = LW_MIN(dist, mindist);
	879	lw_dist2d_comp(g1, lw2, tolerance, dl);
	880	if ( ((tolerance - dl->d)*dl->thedir) >0 )
	881	{
	882	dl->d=tolerance;
	883	return; /* can't be closer */
	884	}
675	885	continue;
676	886	}
677	887
…	…
680	890	uchar *g2 = lwgeom_getsubgeometry_inspected(in2, j);
681	891	int t2 = lwgeom_getType(g2[0]);
682	892
683		~~/* Argument 2 is a multitype... recurse */~~
684	893	if ( lwgeom_contains_subgeoms(t2) )
685	894	{
686		dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance);
687		if ( dist <= tolerance ) return tolerance; /* can't be closer */
688		if ( mindist == -1 ) mindist = dist;
689		else mindist = LW_MIN(dist, mindist);
	895	lw_dist2d_comp(g1, g2, tolerance, dl);
	896	if ( ((tolerance - dl->d)*dl->thedir) >0 )
	897	{
	898	dl->d=tolerance;
	899	return; /* can't be closer */
	900	}
690	901	continue;
691	902	}
692
693	903	if ( t1 == POINTTYPE )
694	904	{
695	905	if ( t2 == POINTTYPE )
696	906	{
697		dist = distance2d_point_point(
698		lwpoint_deserialize(g1),
699		lwpoint_deserialize(g2)
700		);
	907	dl->twisted=1;
	908	lw_dist2d_comp_point_point(
	909	lwpoint_deserialize(g1),
	910	lwpoint_deserialize(g2),
	911	dl);
701	912	}
702	913	else if ( t2 == LINETYPE )
703	914	{
704		dist = distance2d_point_line(
705		lwpoint_deserialize(g1),
706		lwline_deserialize(g2)
707		);
	915	dl->twisted=1;
	916	lw_dist2d_comp_point_line(
	917	lwpoint_deserialize(g1),
	918	lwline_deserialize(g2),
	919	dl);
708	920	}
709	921	else if ( t2 == POLYGONTYPE )
710	922	{
711		dist = distance2d_point_poly(
712		lwpoint_deserialize(g1),
713		lwpoly_deserialize(g2)
714		);
	923	dl->twisted=1;
	924	lw_dist2d_comp_point_poly(
	925	lwpoint_deserialize(g1),
	926	lwpoly_deserialize(g2),
	927	dl);
715	928	}
716	929	else
717	930	{
718		lwerror("Unsupported geometry type: %s", lwgeom_typename(t2));
719		}
	931	lwerror("Unsupported geometry type1: %s", lwgeom_typename(t2));
	932	}
720	933	}
721	934	else if ( t1 == LINETYPE )
722	935	{
723	936	if ( t2 == POINTTYPE )
724	937	{
725		dist = distance2d_point_line(
726		lwpoint_deserialize(g2),
727		lwline_deserialize(g1)
728		);
	938	dl->twisted=(-1);
	939	lw_dist2d_comp_point_line(
	940	lwpoint_deserialize(g2),
	941	lwline_deserialize(g1),
	942	dl);
729	943	}
730	944	else if ( t2 == LINETYPE )
731	945	{
732		dist = distance2d_line_line(
733		lwline_deserialize(g1),
734		lwline_deserialize(g2)
735		);
	946	dl->twisted=1;
	947	/lwnotice("start line");/
	948	lw_dist2d_comp_line_line(
	949	lwline_deserialize(g1),
	950	lwline_deserialize(g2),
	951	dl);
736	952	}
737	953	else if ( t2 == POLYGONTYPE )
738	954	{
739		dist = distance2d_line_poly(
740		lwline_deserialize(g1),
741		lwpoly_deserialize(g2)
742		);
	955	dl->twisted=1;
	956	lw_dist2d_comp_line_poly(
	957	lwline_deserialize(g1),
	958	lwpoly_deserialize(g2),
	959	dl);
743	960	}
744	961	else
745	962	{
746		lwerror("Unsupported geometry type: %s", lwgeom_typename(t2));
747		}
	963	lwerror("Unsupported geometry type2: %s", lwgeom_typename(t2));
	964	}
748	965	}
749	966	else if ( t1 == POLYGONTYPE )
750	967	{
751	968	if ( t2 == POLYGONTYPE )
752	969	{
753		dist = distance2d_poly_poly(
754		lwpoly_deserialize(g2),
755		lwpoly_deserialize(g1)
756		);
	970	dl->twisted=(-1);
	971	lw_dist2d_comp_poly_poly(
	972	lwpoly_deserialize(g2),
	973	lwpoly_deserialize(g1),
	974	dl);
757	975	}
758	976	else if ( t2 == POINTTYPE )
759	977	{
760		dist = distance2d_point_poly(
761		lwpoint_deserialize(g2),
762		lwpoly_deserialize(g1)
763		);
	978	dl->twisted=(-1);
	979	lw_dist2d_comp_point_poly(
	980	lwpoint_deserialize(g2),
	981	lwpoly_deserialize(g1),
	982	dl);
764	983	}
765	984	else if ( t2 == LINETYPE )
766	985	{
767		dist = distance2d_line_poly(
768		lwline_deserialize(g2),
769		lwpoly_deserialize(g1)
770		);
	986	dl->twisted=(-1);
	987	lw_dist2d_comp_line_poly(
	988	lwline_deserialize(g2),
	989	lwpoly_deserialize(g1),
	990	dl);
771	991	}
772	992	else
773	993	{
774		lwerror("Unsupported geometry type: %s", lwgeom_typename(t2));
775		}
	994	lwerror("Unsupported geometry type3: %s", lwgeom_typename(t2));
	995	}
776	996	}
777		// else if (lwgeom_contains_subgeoms(t1)) /* it's a multitype... recurse */
778		// {
779		// dist = lwgeom_mindistance2d_recursive_tolerance(g1, g2, tolerance);
780		// }
	997	/* else if (lwgeom_contains_subgeoms(t1))
	998	{
	999	lw_dist2d_comp(g1, g2, tolerance,dl);
	1000
	1001	}*/
781	1002	else
782	1003	{
783		lwerror("Unsupported geometry type: %s", lwgeom_typename(t1));
	1004	lwerror("Unsupported geometry type4: %s", lwgeom_typename(t1));
784	1005	}
785	1006
786		~~if (mindist == -1 ) mindist = dist;~~
787		~~else mindist = LW_MIN(dist, mindist);~~
788	1007
789		~~LWDEBUGF(3, "dist %d-%d: %f - mindist: %f",~~
790		~~i, j, dist, mindist);~~
791	1008
792	1009
793		if (mindist <= tolerance) return tolerance; /* can't be closer */
	1010	if ( ((tolerance - dl->d)*dl->thedir) >0 )
	1011	{
	1012	dl->d=tolerance;
	1013	return; /* can't be closer */
	1014	}
794	1015
795	1016	}
796	1017
797	1018	}
798
799		if (mindist<0) mindist = 0;
800
801		return mindist;
	1019	return ;
802	1020	}
803	1021
804	1022

postgis/lwgeom_functions_basic.c


43	43	Datum LWGEOM_length_linestring(PG_FUNCTION_ARGS);
44	44	Datum LWGEOM_perimeter2d_poly(PG_FUNCTION_ARGS);
45	45	Datum LWGEOM_perimeter_poly(PG_FUNCTION_ARGS);
	46	Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS);
46	47	Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS);
47		Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS);
	48	Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS);
	49	Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS);
48	50	Datum LWGEOM_inside_circle_point(PG_FUNCTION_ARGS);
49	51	Datum LWGEOM_collect(PG_FUNCTION_ARGS);
50	52	Datum LWGEOM_accum(PG_FUNCTION_ARGS);
…	…
602	604	TYPE_SETZM(nring->dims, 0, 0);
603	605	nring->npoints = ring->npoints;
604	606	nring->serialized_pointlist =
605		lwalloc(ring->npoints*sizeof(POINT2D));
	607	lwalloc(ring->npoints*sizeof(POINT2D));
606	608	loc = nring->serialized_pointlist;
607	609	for (k=0; k<ring->npoints; k++)
608	610	{
…	…
624	626	}
625	627
626	628	if ( type != MULTIPOINTTYPE && type != MULTIPOLYGONTYPE &&
627		type != MULTILINETYPE && type != COLLECTIONTYPE &&
628		type != COMPOUNDTYPE && type != CURVEPOLYTYPE &&
629		type != MULTICURVETYPE && type != MULTISURFACETYPE)
	629	type != MULTILINETYPE && type != COLLECTIONTYPE &&
	630	type != COMPOUNDTYPE && type != CURVEPOLYTYPE &&
	631	type != MULTICURVETYPE && type != MULTISURFACETYPE)
630	632	{
631	633	lwerror("lwgeom_force2d_recursive: unknown geometry: %d",
632	634	type);
…	…
818	820	TYPE_SETZM(nring->dims, 1, 0);
819	821	nring->npoints = ring->npoints;
820	822	nring->serialized_pointlist =
821		lwalloc(ring->npoints*sizeof(POINT3DZ));
	823	lwalloc(ring->npoints*sizeof(POINT3DZ));
822	824	loc = nring->serialized_pointlist;
823	825	for (k=0; k<ring->npoints; k++)
824	826	{
…	…
1012	1014	TYPE_SETZM(nring->dims, 0, 1);
1013	1015	nring->npoints = ring->npoints;
1014	1016	nring->serialized_pointlist =
1015		lwalloc(ring->npoints*sizeof(POINT3DM));
	1017	lwalloc(ring->npoints*sizeof(POINT3DM));
1016	1018	loc = nring->serialized_pointlist;
1017	1019	for (k=0; k<ring->npoints; k++)
1018	1020	{
…	…
1034	1036	}
1035	1037
1036	1038	if ( type != MULTIPOINTTYPE && type != MULTIPOLYGONTYPE &&
1037		type != MULTILINETYPE && type != COLLECTIONTYPE &&
1038		type != COMPOUNDTYPE && type != CURVEPOLYTYPE &&
1039		type != MULTICURVETYPE && type != MULTISURFACETYPE)
	1039	type != MULTILINETYPE && type != COLLECTIONTYPE &&
	1040	type != COMPOUNDTYPE && type != CURVEPOLYTYPE &&
	1041	type != MULTICURVETYPE && type != MULTISURFACETYPE)
1040	1042	{
1041	1043	lwerror("lwgeom_force3dm_recursive: unknown geometry: %d",
1042	1044	type);
…	…
1224	1226	TYPE_SETZM(nring->dims, 1, 1);
1225	1227	nring->npoints = ring->npoints;
1226	1228	nring->serialized_pointlist =
1227		lwalloc(ring->npoints*sizeof(POINT4D));
	1229	lwalloc(ring->npoints*sizeof(POINT4D));
1228	1230	loc = nring->serialized_pointlist;
1229	1231	for (k=0; k<ring->npoints; k++)
1230	1232	{
…	…
1252	1254
1253	1255	/* Add type */
1254	1256	*optr = lwgeom_makeType_full(
1255		1, 1,
1256		lwgeom_hasSRID(serialized[0]),
1257		type, lwgeom_hasBBOX(serialized[0]));
	1257	1, 1,
	1258	lwgeom_hasSRID(serialized[0]),
	1259	type, lwgeom_hasBBOX(serialized[0]));
1258	1260	optr++;
1259	1261	totsize++;
1260	1262	loc=serialized+1;
…	…
1449	1451	* automatic bbox addition FOR_COMPLEX_GEOMS.
1450	1452	*/
1451	1453	if ( TYPE_GETTYPE(geom->type) == COLLECTIONTYPE &&
1452		TYPE_HASBBOX(geom->type) )
	1454	TYPE_HASBBOX(geom->type) )
1453	1455	{
1454	1456	PG_RETURN_POINTER(geom);
1455	1457	}
…	…
1473	1475	lwgeom->bbox = NULL;
1474	1476	lwgeoms[0] = lwgeom;
1475	1477	lwgeom = (LWGEOM *)lwcollection_construct(COLLECTIONTYPE,
1476		SRID, bbox, 1,
1477		lwgeoms);
	1478	SRID, bbox, 1,
	1479	lwgeoms);
1478	1480	}
1479	1481
1480	1482	result = pglwgeom_serialize(lwgeom);
…	…
1510	1512	/* deserialize into lwgeoms[0] */
1511	1513	lwgeom = lwgeom_deserialize(SERIALIZED_FORM(geom));
1512	1514	ogeom = lwgeom_as_multi(lwgeom);
1513		printf("ogeom %p\n",ogeom);
1514		printf("ogeom->type %d\n", ogeom->type);
	1515	printf("ogeom %p\n",ogeom);
	1516	printf("ogeom->type %d\n", ogeom->type);
1515	1517
1516	1518	result = pglwgeom_serialize(ogeom);
1517	1519
…	…
1520	1522	PG_RETURN_POINTER(result);
1521	1523	}
1522	1524
	1525
	1526	PG_FUNCTION_INFO_V1(LWGEOM_shortestline2d);
	1527	Datum LWGEOM_shortestline2d(PG_FUNCTION_ARGS)
	1528	{
	1529
	1530	PG_LWGEOM *geom1;
	1531	PG_LWGEOM *geom2;
	1532	LWPOINT *point1;
	1533	LWPOINT *point2;
	1534	int SRID;
	1535	double x1,x2,y1,y2;
	1536	PG_LWGEOM *result=NULL;
	1537	LWPOINT *lwpoints[2];
	1538	LWLINE *outline;
	1539
	1540	DISTPTS thedl;
	1541	DISTPTS *dl = &thedl;
	1542	dl->d=99999999;
	1543	dl->thedir=(1);
	1544
	1545	geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
	1546	geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
	1547
	1548	if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2))
	1549	{
	1550	elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
	1551	PG_RETURN_NULL();
	1552	}
	1553
	1554	SRID = pglwgeom_getSRID(geom1);
	1555
	1556	lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2),0.0, dl );
	1557
	1558	x1=dl->p1.x;
	1559	y1=dl->p1.y;
	1560	x2=dl->p2.x;
	1561	y2=dl->p2.y;
	1562	/*
	1563	if(x1==0)
	1564	{
	1565	if((x2==0)&&(y1==0)&&(y2==0))
	1566	{
	1567	PG_RETURN_NULL();
	1568	}
	1569	}
	1570	*/
	1571	point1 = make_lwpoint2d(SRID, x1, y1);
	1572	point2 = make_lwpoint2d(SRID, x2, y2);
	1573
	1574	lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1)));
	1575	lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2)));
	1576
	1577	outline = lwline_from_lwpointarray(SRID, 2, lwpoints);
	1578
	1579	result = pglwgeom_serialize((LWGEOM *)outline);
	1580	PG_RETURN_POINTER(result);
	1581
	1582	}
	1583
	1584	PG_FUNCTION_INFO_V1(LWGEOM_longestline2d);
	1585	Datum LWGEOM_longestline2d(PG_FUNCTION_ARGS)
	1586	{
	1587
	1588	PG_LWGEOM *geom1;
	1589	PG_LWGEOM *geom2;
	1590	LWPOINT *point1;
	1591	LWPOINT *point2;
	1592	int SRID;
	1593
	1594	double x1,x2,y1,y2;
	1595	PG_LWGEOM *result=NULL;
	1596	LWPOINT *lwpoints[2];
	1597	LWLINE *outline;
	1598	DISTPTS thedl;
	1599	DISTPTS *dl = &thedl;
	1600	dl->d=0;
	1601	dl->thedir=(-1);
	1602
	1603	geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
	1604	geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
	1605
	1606	if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2))
	1607	{
	1608	elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
	1609	PG_RETURN_NULL();
	1610	}
	1611
	1612	SRID = pglwgeom_getSRID(geom1);
	1613
	1614	lw_dist2d_comp( SERIALIZED_FORM(geom1),SERIALIZED_FORM(geom2),99999999, dl );
	1615
	1616	x1=dl->p1.x;
	1617	y1=dl->p1.y;
	1618	x2=dl->p2.x;
	1619	y2=dl->p2.y;
	1620
	1621	point1 = make_lwpoint2d(SRID, x1, y1);
	1622	point2 = make_lwpoint2d(SRID, x2, y2);
	1623
	1624	lwpoints[0] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point1)));
	1625	lwpoints[1] = lwpoint_deserialize(SERIALIZED_FORM(pglwgeom_serialize((LWGEOM *)point2)));
	1626
	1627	outline = lwline_from_lwpointarray(SRID, 2, lwpoints);
	1628
	1629	result = pglwgeom_serialize((LWGEOM *)outline);
	1630
	1631	PG_RETURN_POINTER(result);
	1632	}
	1633
1523	1634	/* Minimum 2d distance between objects in geom1 and geom2. */
1524	1635	PG_FUNCTION_INFO_V1(LWGEOM_mindistance2d);
1525	1636	Datum LWGEOM_mindistance2d(PG_FUNCTION_ARGS)
…	…
1540	1651	}
1541	1652
1542	1653	mindist = lwgeom_mindistance2d_recursive(SERIALIZED_FORM(geom1),
1543		SERIALIZED_FORM(geom2));
	1654	SERIALIZED_FORM(geom2));
1544	1655
1545	1656	PROFSTOP(PROF_QRUN);
1546	1657	PROFREPORT("dist",geom1, geom2, NULL);
…	…
1551	1662	PG_RETURN_FLOAT8(mindist);
1552	1663	}
1553	1664
1554		/* Minimum 2d distance between objects in geom1 and geom2. */
1555		PG_FUNCTION_INFO_V1(LWGEOM_dwithin);
1556		Datum LWGEOM_dwithin(PG_FUNCTION_ARGS)
	1665	PG_FUNCTION_INFO_V1(LWGEOM_maxdistance2d_linestring);
	1666	Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS)
1557	1667	{
1558	1668	PG_LWGEOM *geom1;
1559	1669	PG_LWGEOM *geom2;
1560		double m~~indist, tolerance~~;
	1670	double maxdist;
1561	1671
1562	1672	PROFSTART(PROF_QRUN);
1563	1673
1564	1674	geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1565	1675	geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1566		~~tolerance = PG_GETARG_FLOAT8(2);~~
1567	1676
1568		~~if ( tolerance < 0 )~~
1569		{
1570		~~elog(ERROR,"Tolerance cannot be less than zero\n");~~
1571		~~PG_RETURN_NULL();~~
1572		}
1573
1574	1677	if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2))
1575	1678	{
1576	1679	elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
1577	1680	PG_RETURN_NULL();
1578	1681	}
1579	1682
1580		mindist = lwgeom_mindistance2d_recursive_tolerance(
1581		SERIALIZED_FORM(geom1),
1582		SERIALIZED_FORM(geom2),
1583		tolerance
1584		);
	1683	maxdist = lwgeom_maxdistance2d_recursive(SERIALIZED_FORM(geom1),
	1684	SERIALIZED_FORM(geom2));
1585	1685
1586	1686	PROFSTOP(PROF_QRUN);
1587		PROFREPORT("dist",geom1, geom2, NULL);
	1687	PROFREPORT("maxdist",geom1, geom2, NULL);
1588	1688
1589	1689	PG_FREE_IF_COPY(geom1, 0);
1590	1690	PG_FREE_IF_COPY(geom2, 1);
1591	1691
1592		PG_RETURN_~~BOOL(tolerance >= min~~dist);
	1692	PG_RETURN_FLOAT8(maxdist);
1593	1693	}
1594	1694
1595		/*
1596		* Maximum 2d distance between linestrings.
1597		* Returns NULL if given geoms are not linestrings.
1598		* This is very bogus (or I'm missing its meaning)
1599		*/
1600		PG_FUNCTION_INFO_V1(LWGEOM_maxdistance2d_linestring);
1601		Datum LWGEOM_maxdistance2d_linestring(PG_FUNCTION_ARGS)
	1695	/* Minimum 2d distance between objects in geom1 and geom2. */
	1696	PG_FUNCTION_INFO_V1(LWGEOM_dwithin);
	1697	Datum LWGEOM_dwithin(PG_FUNCTION_ARGS)
1602	1698	{
1603
1604	1699	PG_LWGEOM *geom1;
1605	1700	PG_LWGEOM *geom2;
1606		LWLINE *line1;
1607		LWLINE *line2;
1608		double maxdist = 0;
1609		int i;
	1701	double mindist, tolerance;
1610	1702
1611		elog(ERROR, "This function is unimplemented yet");
1612		PG_RETURN_NULL();
	1703	PROFSTART(PROF_QRUN);
1613	1704
1614	1705	geom1 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
1615		~~line1 = lwline_deserialize(SERIALIZED_FORM(geom1));~~
1616		~~if ( line1 == NULL ) PG_RETURN_NULL(); /* not a linestring */~~
1617
1618	1706	geom2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
1619		line2 = lwline_deserialize(SERIALIZED_FORM(geom2));
1620		if ( line2 == NULL ) PG_RETURN_NULL(); /* not a linestring */
	1707	tolerance = PG_GETARG_FLOAT8(2);
1621	1708
	1709	if ( tolerance < 0 )
	1710	{
	1711	elog(ERROR,"Tolerance cannot be less than zero\n");
	1712	PG_RETURN_NULL();
	1713	}
	1714
1622	1715	if (pglwgeom_getSRID(geom1) != pglwgeom_getSRID(geom2))
1623	1716	{
1624	1717	elog(ERROR,"Operation on two GEOMETRIES with different SRIDs\n");
1625	1718	PG_RETURN_NULL();
1626	1719	}
1627	1720
1628		for (i=0; i<line1->points->npoints; i++)
1629		{
1630		POINT2D p;
1631		double dist;
	1721	mindist = lwgeom_mindistance2d_recursive_tolerance(
	1722	SERIALIZED_FORM(geom1),
	1723	SERIALIZED_FORM(geom2),
	1724	tolerance
	1725	);
1632	1726
1633		~~getPoint2d_p(line1->points, i, &p~~);
1634		~~dist = distance2d_pt_ptarray(&p, line2->points~~);
	1727	PROFSTOP(PROF_QRUN);
	1728	PROFREPORT("dist",geom1, geom2, NULL);
1635	1729
1636		~~if (dist > maxdist) maxdist = dist;~~
1637		}
1638
1639	1730	PG_FREE_IF_COPY(geom1, 0);
1640	1731	PG_FREE_IF_COPY(geom2, 1);
1641	1732
1642		PG_RETURN_~~FLOAT8(max~~dist);
	1733	PG_RETURN_BOOL(tolerance >= mindist);
1643	1734	}
1644	1735
1645	1736	/*
	1737	* LWGEOM_maxdistance2d_linestring is deleted and I have a new aproach to it /Nicklas Avén
	1738	*/
	1739	/*
1646	1740	* Longitude shift:
1647	1741	* Y remains the same
1648	1742	* X is converted:
…	…
1788	1882	lwgeom_dropSRID(lwgeoms[1]);
1789	1883
1790	1884	outlwg = (LWGEOM *)lwcollection_construct(
1791		outtype, SRID,
1792		box, 2, lwgeoms);
	1885	outtype, SRID,
	1886	box, 2, lwgeoms);
1793	1887
1794	1888	result = pglwgeom_serialize(outlwg);
1795	1889
…	…
2062	2156	POSTGIS_DEBUGF(3, "LWGEOM_collect_garray: outtype = %d", outtype);
2063	2157
2064	2158	outlwg = (LWGEOM *)lwcollection_construct(
2065		outtype, SRID,
2066		box, nelems, lwgeoms);
	2159	outtype, SRID,
	2160	box, nelems, lwgeoms);
2067	2161
2068	2162	result = pglwgeom_serialize(outlwg);
2069	2163
…	…
2175	2269	if ( TYPE_GETTYPE(geom->type) != POINTTYPE ) continue;
2176	2270
2177	2271	lwpoints[npoints++] =
2178		lwpoint_deserialize(SERIALIZED_FORM(geom));
	2272	lwpoint_deserialize(SERIALIZED_FORM(geom));
2179	2273
2180	2274	/* Check SRID homogeneity */
2181	2275	if ( npoints == 1 )
…	…
2232	2326	pglwg2 = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
2233	2327
2234	2328	if ( ! TYPE_GETTYPE(pglwg1->type) == POINTTYPE \|\|
2235		! TYPE_GETTYPE(pglwg2->type) == POINTTYPE )
	2329	! TYPE_GETTYPE(pglwg2->type) == POINTTYPE )
2236	2330	{
2237	2331	elog(ERROR, "Input geometries must be points");
2238	2332	PG_RETURN_NULL();
…	…
2436	2530
2437	2531
2438	2532	if (box.xmin == box.xmax &&
2439		box.ymin == box.ymax)
	2533	box.ymin == box.ymax)
2440	2534	{
2441	2535	/* Construct and serialize point */
2442	2536	LWPOINT *point = make_lwpoint2d(SRID, box.xmin, box.ymin);
2443	2537	ser = lwpoint_serialize(point);
2444	2538	}
2445	2539	else if (box.xmin == box.xmax \|\|
2446		box.ymin == box.ymax)
	2540	box.ymin == box.ymax)
2447	2541	{
2448	2542	LWLINE *line;
2449	2543	POINT2D pts = palloc(sizeof(POINT2D)2);
…	…
2533	2627
2534	2628	/* Avoid deserialize/serialize steps */
2535	2629	if ( (TYPE_GETTYPE(ingeom->type) == POINTTYPE) \|\|
2536		(TYPE_GETTYPE(ingeom->type) == MULTIPOINTTYPE) )
	2630	(TYPE_GETTYPE(ingeom->type) == MULTIPOINTTYPE) )
2537	2631	PG_RETURN_POINTER(ingeom);
2538	2632
2539	2633	inlwgeom = lwgeom_deserialize(SERIALIZED_FORM(ingeom));
…	…
3080	3174	g1_bvol.ymax = g1_bvol.ymax + dist;
3081	3175
3082	3176	if ( (g1_bvol.xmin > geom2->bbox->xmax) \|\|
3083		(g1_bvol.xmax < geom2->bbox->xmin) \|\|
3084		(g1_bvol.ymin > geom2->bbox->ymax) \|\|
3085		(g1_bvol.ymax < geom2->bbox->ymin)
	3177	(g1_bvol.xmax < geom2->bbox->xmin) \|\|
	3178	(g1_bvol.ymin > geom2->bbox->ymax) \|\|
	3179	(g1_bvol.ymax < geom2->bbox->ymin)
3086	3180	)
3087	3181	{
3088	3182	PG_RETURN_BOOL(FALSE); /bbox not overlap /
…	…
3216	3310	if (curve != NULL)
3217	3311	{
3218	3312	lwgeom_affine_ptarray(curve->points,
3219		afac, bfac, cfac,
3220		dfac, efac, ffac,
3221		gfac, hfac, ifac,
3222		xoff, yoff, zoff);
	3313	afac, bfac, cfac,
	3314	dfac, efac, ffac,
	3315	gfac, hfac, ifac,
	3316	xoff, yoff, zoff);
3223	3317	lwgeom_release((LWGEOM *)curve);
3224	3318	continue;
3225	3319	}
…	…
3287	3381	Datum ST_GeoHash(PG_FUNCTION_ARGS)
3288	3382	{
3289	3383
3290		PG_LWGEOM *geom = NULL;
3291		int precision = 0;
3292		int len = 0;
3293		char *geohash = NULL;
3294		char *result = NULL;
	3384	PG_LWGEOM *geom = NULL;
	3385	int precision = 0;
	3386	int len = 0;
	3387	char *geohash = NULL;
	3388	char *result = NULL;
3295	3389
3296		if( PG_ARGISNULL(0) )
3297		{
	3390	if ( PG_ARGISNULL(0) )
	3391	{
3298	3392	PG_RETURN_NULL();
3299		}
	3393	}
3300	3394
3301		geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
	3395	geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
3302	3396
3303		if( ! PG_ARGISNULL(1) )
3304		{
3305		precision = PG_GETARG_INT32(1);
3306		}
	3397	if ( ! PG_ARGISNULL(1) )
	3398	{
	3399	precision = PG_GETARG_INT32(1);
	3400	}
3307	3401
3308		geohash = lwgeom_geohash((LWGEOM*)(pglwgeom_deserialize(geom)), precision);
	3402	geohash = lwgeom_geohash((LWGEOM*)(pglwgeom_deserialize(geom)), precision);
3309	3403
3310		if( ! geohash )
3311		{
	3404	if ( ! geohash )
	3405	{
3312	3406	elog(ERROR,"ST_GeoHash: lwgeom_geohash returned NULL.\n");
3313	3407	PG_RETURN_NULL();
3314		}
	3408	}
3315	3409
3316	3410	len = strlen(geohash) + VARHDRSZ;
3317		result = palloc(len);
	3411	result = palloc(len);
3318	3412	SET_VARSIZE(result, len);
3319	3413	memcpy(VARDATA(result), geohash, len-VARHDRSZ);
3320	3414	pfree(geohash);
3321	3415	PG_RETURN_POINTER(result);
3322
	3416
3323	3417	}

Download in other formats:

Original Format