#3074, first piece of infrastructure, count points in a piece

git-svn-id: http://svn.osgeo.org/postgis/trunk@13322 b70326c6-7e19-0410-871a-916f4a2858ee

liblwgeom/cunit/cu_misc.c

@@ -139,6 +139,24 @@ static void test_grid(void)
 	lwgeom_free(geomgrid);		
 }
 
+static void test_rect_count(void)
+{
+	GBOX box;
+	int n;
+	static char *wkt = "MULTIPOLYGON(((0 0, 10 0, 10 10, 0 10, 0 0)))";
+	LWGEOM *geom = lwgeom_from_wkt(wkt, LW_PARSER_CHECK_ALL);
+
+	box.xmin = -5;  box.ymin = -5;
+	box.xmax =  5;  box.ymax =  5;
+	n = lwgeom_npoints_in_rect(geom, &box);
+	CU_ASSERT_EQUAL(2, n);
+
+	box.xmin = -5;  box.ymin = -5;
+	box.xmax = 15;  box.ymax = 15; 
+	n = lwgeom_npoints_in_rect(geom, &box);
+	CU_ASSERT_EQUAL(5, n);
+}
+
 
 /*
 ** Used by the test harness to register the tests in this file.
@@ -153,4 +171,5 @@ void misc_suite_setup(void)
 	PG_ADD_TEST(suite, test_misc_area);
 	PG_ADD_TEST(suite, test_misc_wkb);
 	PG_ADD_TEST(suite, test_grid);
+	PG_ADD_TEST(suite, test_rect_count);
 }

liblwgeom/g_box.c

@@ -307,6 +307,17 @@ gbox_contains_2d(const GBOX *g1, const GBOX *g2)
 	return LW_TRUE;
 }
 
+int 
+gbox_contains_point2d(const GBOX *g, const POINT2D *p)
+{
+	if ( ( g->xmin <= p->x ) && ( g->xmax >= p->x ) &&
+	     ( g->ymin <= p->y ) && ( g->ymax >= p->y ) )
+	{
+		return LW_TRUE;
+	}
+	return LW_FALSE;
+}
+
 /**
 * Warning, this function is only good for x/y/z boxes, used
 * in unit testing of geodetic box generation.

liblwgeom/liblwgeom_internal.h

@@ -435,3 +435,7 @@ extern int _lwgeom_interrupt_requested;
 }
 
 #endif /* _LIBLWGEOM_INTERNAL_H */
+
+int ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox);
+int gbox_contains_point2d(const GBOX *g, const POINT2D *p);
+int lwgeom_npoints_in_rect(const LWGEOM *geom, const GBOX *gbox);

liblwgeom/lwgeom.c

@@ -1818,3 +1818,57 @@ lwgeom_grid(const LWGEOM *lwgeom, const gridspec *grid)
 			return NULL;
 	}
 }
+
+
+int
+lwgeom_npoints_in_rect(const LWGEOM *geom, const GBOX *gbox)
+{
+	const GBOX *geombox = lwgeom_get_bbox(geom);
+
+	if ( ! gbox_overlaps_2d(geombox, gbox) )
+		return 0;
+
+	switch ( geom->type )
+	{
+		case POINTTYPE:
+			return ptarray_npoints_in_rect(((LWPOINT*)geom)->point, gbox);
+
+		case CIRCSTRINGTYPE:
+		case LINETYPE:
+			return ptarray_npoints_in_rect(((LWLINE*)geom)->points, gbox);
+
+		case POLYGONTYPE:
+		{
+			int i, n = 0;
+			LWPOLY *poly = (LWPOLY*)geom;
+			for ( i = 0; i < poly->nrings; i++ )
+			{
+				n += ptarray_npoints_in_rect(poly->rings[i], gbox);
+			}
+			return n;
+		}
+
+		case MULTIPOINTTYPE:
+		case MULTILINETYPE:
+		case MULTIPOLYGONTYPE:
+		case COMPOUNDTYPE:
+		case COLLECTIONTYPE:
+		{
+			int i, n = 0;
+			LWCOLLECTION *col = (LWCOLLECTION*)geom;
+			for ( i = 0; i < col->ngeoms; i++ )
+			{
+				n += lwgeom_npoints_in_rect(col->geoms[i], gbox);
+			}
+			return n;
+		}
+		default:
+		{
+			lwerror("lwgeom_npoints_in_rect: Unsupported geometry type: %s",
+			        lwtype_name(geom->type));
+			return 0;
+		}
+	}
+	return 0;
+}
+

liblwgeom/ptarray.c

@@ -1818,3 +1818,19 @@ ptarray_grid(const POINTARRAY *pa, const gridspec *grid)
 	return dpa;
 }
 
+int 
+ptarray_npoints_in_rect(const POINTARRAY *pa, const GBOX *gbox)
+{
+	const POINT2D *pt;
+	int n = 0;
+	int i;
+	for ( i = 0; i < pa->npoints; i++ )
+	{
+		pt = getPoint2d_cp(pa, i);
+		if ( gbox_contains_point2d(gbox, pt) )
+			n++;
+	}
+	return n;
+}
+
+