[postgis-tickets] r14612 - #3433, ST_ClusterIntersecting gives incorrect result for MultiPoints

[Daniel Baston]

Author: dbaston
Date: 2016-01-20 04:18:14 -0800 (Wed, 20 Jan 2016)
New Revision: 14612

Modified:
   branches/2.2/NEWS
   branches/2.2/liblwgeom/cunit/cu_geos_cluster.c
   branches/2.2/liblwgeom/lwgeom_geos_cluster.c
Log:
#3433, ST_ClusterIntersecting gives incorrect result for MultiPoints

Modified: branches/2.2/NEWS
===================================================================
--- branches/2.2/NEWS	2016-01-20 11:57:21 UTC (rev 14611)
+++ branches/2.2/NEWS	2016-01-20 12:18:14 UTC (rev 14612)
@@ -6,6 +6,7 @@
   - #3422, Improve ST_Split robustness on standard precision double
            systems (arm64, ppc64el, s390c, powerpc, ...)
   - #3427, Update spatial_ref_sys to EPSG version 8.8
+  - #3433, ST_ClusterIntersecting incorrect for MultiPoints
   - #3436, memory handling mistake in ptarray_clone_deep
 
 PostGIS 2.2.1

Modified: branches/2.2/liblwgeom/cunit/cu_geos_cluster.c
===================================================================
--- branches/2.2/liblwgeom/cunit/cu_geos_cluster.c	2016-01-20 11:57:21 UTC (rev 14611)
+++ branches/2.2/liblwgeom/cunit/cu_geos_cluster.c	2016-01-20 12:18:14 UTC (rev 14612)
@@ -238,6 +238,23 @@
 	perform_cluster_within_distance_test(1, wkt_inputs, 2, expected_outputs, 2);
 }
 
+static void multipoint_test(void)
+{
+	/* See #3433 */
+	char* wkt_inputs_mp[] = { "MULTIPOINT ((0 0), (0 1))", "POINT (0 0)"};
+	char* expected_outputs_mp[] = { "GEOMETRYCOLLECTION(MULTIPOINT ((0 0), (0 1)), POINT (0 0))"};
+
+	char* wkt_inputs_gc[] = { "GEOMETRYCOLLECTION (POINT (0 0), POINT (0 1))", "POINT (0 0)"};
+	char* expected_outputs_gc[] = { "GEOMETRYCOLLECTION(GEOMETRYCOLLECTION (POINT (0 0), POINT (0 1)), POINT (0 0))"};
+
+	char* wkt_inputs_pt[] = { "POINT (3 3)", "POINT (3 3)"};
+	char* expected_outputs_pt[] = { "GEOMETRYCOLLECTION(POINT (3 3), POINT (3 3))"};
+
+	perform_cluster_intersecting_test(wkt_inputs_mp, 2, expected_outputs_mp, 1);
+	perform_cluster_intersecting_test(wkt_inputs_gc, 2, expected_outputs_gc, 1);
+	perform_cluster_intersecting_test(wkt_inputs_pt, 2, expected_outputs_pt, 1);
+}
+
 void geos_cluster_suite_setup(void);
 void geos_cluster_suite_setup(void)
 {
@@ -247,4 +264,5 @@
 	PG_ADD_TEST(suite, basic_distance_test);
 	PG_ADD_TEST(suite, single_input_test);
 	PG_ADD_TEST(suite, empty_inputs_test);
+	PG_ADD_TEST(suite, multipoint_test);
 }

Modified: branches/2.2/liblwgeom/lwgeom_geos_cluster.c
===================================================================
--- branches/2.2/liblwgeom/lwgeom_geos_cluster.c	2016-01-20 11:57:21 UTC (rev 14611)
+++ branches/2.2/liblwgeom/lwgeom_geos_cluster.c	2016-01-20 12:18:14 UTC (rev 14612)
@@ -124,12 +124,26 @@
 
 	if (p != q && UF_find(cxt->uf, p) != UF_find(cxt->uf, q))
 	{
-		/* Lazy initialize prepared geometry */
-		if (cxt->prep == NULL)
+		int geos_type = GEOSGeomTypeId(cxt->geoms[p]);
+		int geos_result;
+
+		/* Don't build prepared a geometry around a Point or MultiPoint -
+		 * there are some problems in the implementation, and it's not clear
+		 * there would be a performance benefit in any case.  (See #3433)
+		 */
+		if (geos_type != GEOS_POINT && geos_type != GEOS_MULTIPOINT)
 		{
-			cxt->prep = GEOSPrepare(cxt->geoms[p]);
+			/* Lazy initialize prepared geometry */
+			if (cxt->prep == NULL)
+			{
+				cxt->prep = GEOSPrepare(cxt->geoms[p]);
+			}
+			geos_result = GEOSPreparedIntersects(cxt->prep, cxt->geoms[q]);
 		}
-		int geos_result = GEOSPreparedIntersects(cxt->prep, cxt->geoms[q]);
+		else
+		{
+			geos_result = GEOSIntersects(cxt->geoms[p], cxt->geoms[q]);
+		}
 		if (geos_result > 1)
 		{
 			cxt->error = geos_result;