xyzm (rebased)

.gitignore

@@ -15,3 +15,5 @@ Makefile
 */*/Makefile
 pgsql/pointcloud.control
 pgsql/pointcloud--*
+pgsql_postgis/pointcloud_postgis--*.sql
+pgsql_postgis/pointcloud_postgis.control

README.md

@@ -276,12 +276,12 @@ Now that you have created two tables, you'll see entries for them in the `pointc
 > 
 >     1     
 
-**PC_Envelope(p pcpatch)** returns **bytea**
+**PC_Envelope_AsBinary(p pcpatch)** returns **bytea**
 
 > Return the OGC "well-known binary" format for *bounds* of the patch.
 > Useful for performing intersection tests with geometries.
 > 
->     SELECT PC_Envelope(pa) FROM patches LIMIT 1;
+>     SELECT PC_Envelope_AsBinary(pa) FROM patches LIMIT 1;
 >
 >     \x0103000000010000000500000090c2f5285cbf5fc0e17a
 >     14ae4781464090c2f5285cbf5fc0ec51b81e858b46400ad7
@@ -519,6 +519,11 @@ The `pointcloud_postgis` extension adds functions that allow you to use PostgreS
 > 
 >     POINT Z (-127 45 124)
 
+**PC_Envelope(pcpatch)** returns **geometry**<br/>
+**pcpatch::geometry** returns **geometry**
+
+> Get the PcPatch bounds as a PostGIS geometry
+
 ## Compressions ##
 
 One of the issues with LIDAR data is that there is a lot of it. To deal with data volumes, PostgreSQL Pointcloud allows schemas to declare their preferred compression method in the `<pc:metadata>` block of the schema document. In the example schema, we declared our compression as follows:

lib/cunit/cu_pc_patch.c

@@ -454,8 +454,11 @@ test_patch_wkb()
 	PCPOINTLIST *pl1;
 	PCPATCH_UNCOMPRESSED *pu1, *pu2;
 	PCPATCH *pa1, *pa2, *pa3, *pa4;
-	size_t z1, z2;
-	uint8_t *wkb1, *wkb2;
+	size_t z1, z2, z3;
+	uint8_t *wkb1, *wkb2, *wkb3, *hexwkb;
+
+	static char *hexresult_ndr = "01030000000100000005000000000000000000000000000000000000000000000000000000CDCCCCCCCC8C4B40EC51B81E852B4440CDCCCCCCCC8C4B40EC51B81E852B4440000000000000000000000000000000000000000000000000";
+	static char *hexresult_xdr = "00000000030000000100000005000000000000000000000000000000000000000000000000404B8CCCCCCCCCCD40442B851EB851EC404B8CCCCCCCCCCD40442B851EB851EC000000000000000000000000000000000000000000000000";
 
 	pl1 = pc_pointlist_make(npts);
 
@@ -474,6 +477,7 @@ test_patch_wkb()
 	// str = hexbytes_from_bytes(wkb1, z1);
 	// printf("str\n%s\n",str);
 	pa2 = pc_patch_from_wkb(simpleschema, wkb1, z1);
+	pcfree(wkb1);
 
 	// printf("pa2\n%s\n",pc_patch_to_string(pa2));
 
@@ -497,6 +501,18 @@ test_patch_wkb()
 	CU_ASSERT_EQUAL(pu1->npoints, pu2->npoints);
 	CU_ASSERT(memcmp(pu1->data, pu2->data, pu1->datasize) == 0);
 
+	wkb3 = pc_bounds_to_geometry_wkb(&pa1->bounds,simpleschema->srid,&z3);
+	hexwkb = hexbytes_from_bytes(wkb3,z3);
+	if ( machine_endian() == PC_NDR )
+	{
+		CU_ASSERT_STRING_EQUAL(hexwkb, hexresult_ndr);
+	}
+	else
+	{
+		CU_ASSERT_STRING_EQUAL(hexwkb, hexresult_xdr);
+	}
+	pcfree(hexwkb);
+	pcfree(wkb3);
 
 	pc_pointlist_free(pl1);
 	pc_patch_free(pa1);
@@ -505,7 +521,6 @@ test_patch_wkb()
 	pc_patch_free(pa4);
 	pc_patch_free((PCPATCH*)pu1);
 	pc_patch_free((PCPATCH*)pu2);
-	pcfree(wkb1);
 }
 
 

lib/cunit/cu_pc_schema.c

@@ -205,9 +205,9 @@ test_schema_clone(void)
 	pc_schema_free(clone);
 
 	/* See https://github.com/pgpointcloud/pointcloud/issues/66 */
-	  xmlstr = "<pc:PointCloudSchema xmlns:pc='x'><pc:dimension><pc:position>1</pc:position></pc:dimension></pc:PointCloudSchema>";
+	xmlstr = "<pc:PointCloudSchema xmlns:pc='x'><pc:dimension><pc:position>1</pc:position></pc:dimension></pc:PointCloudSchema>";
 	i = pc_schema_from_xml(xmlstr, &myschema);
-	  CU_ASSERT_EQUAL(i, PC_SUCCESS);
+	CU_ASSERT_EQUAL(i, PC_SUCCESS);
 	clone = pc_schema_clone(myschema);
 	CU_ASSERT_EQUAL(clone->ndims, myschema->ndims);
 	CU_ASSERT_EQUAL(clone->dims[0]->name, NULL);

lib/pc_api.h

@@ -150,8 +150,19 @@ typedef struct
 	double xmax;
 	double ymin;
 	double ymax;
+	double zmin;
+	double zmax;
+	double mmin;
+	double mmax;
 } PCBOUNDS;
 
+typedef struct
+{
+	double xmin, ymin, zmin;
+	double xmax, ymax, zmax;
+	int32_t srid;
+} PCBOX3D;
+
 /* Used for generic patch statistics */
 typedef struct
 {
@@ -437,6 +448,12 @@ int pc_patch_compute_extent(PCPATCH *patch);
 /** True/false if bounds intersect */
 int pc_bounds_intersects(const PCBOUNDS *b1, const PCBOUNDS *b2);
 
+/** Return the bounds as an OGC WKB geometry */
+uint8_t *pc_bounds_to_geometry_wkb(const PCBOUNDS *bounds, uint32_t srid, size_t *wkbsize);
+
+/** Return the bounds as a BOX3D */
+PCBOX3D *pc_bounds_to_box3d(const PCBOUNDS *bounds, uint32_t srid);
+
 /** Subset batch based on less-than condition on dimension */
 PCPATCH* pc_patch_filter_lt_by_name(const PCPATCH *pa, const char *name, double val);
 

lib/pc_api_internal.h

@@ -269,11 +269,13 @@ void pc_bytes_to_ptr(uint8_t *buf, PCBYTES pcb, int n);
 */
 
 /** Initialize with very large mins and very small maxes */
-void pc_bounds_init(PCBOUNDS *b);
+void pc_bounds_init(PCBOUNDS *b, const PCSCHEMA *schema);
 /** Copy a bounds */
 PCSTATS* pc_stats_clone(const PCSTATS *stats);
 /** Expand extents of b1 to encompass b2 */
 void pc_bounds_merge(PCBOUNDS *b1, const PCBOUNDS *b2);
+/** Expand extents of b to encompass p */
+void pc_bounds_expand(PCBOUNDS *b, const PCPOINT *p);
 
 /****************************************************************************
 * BITMAPS

lib/pc_filter.c

@@ -165,7 +165,7 @@ pc_patch_dimensional_filter(const PCPATCH_DIMENSIONAL *pdl, const PCBITMAP *map)
 		stats.max = pc_value_scale_offset(stats.max, dim);
 		stats.sum = pc_value_scale_offset(stats.sum, dim);
 
-		/* Save the X/Y stats for use in bounds later */
+		/* Save the XYZM stats for use in bounds later */
 		if ( i == pdl->schema->x_position )
 		{
 			fpdl->bounds.xmin = stats.min;
@@ -176,6 +176,16 @@ pc_patch_dimensional_filter(const PCPATCH_DIMENSIONAL *pdl, const PCBITMAP *map)
 			fpdl->bounds.ymin = stats.min;
 			fpdl->bounds.ymax = stats.max;
 		}
+		else if ( i == pdl->schema->z_position )
+		{
+			fpdl->bounds.zmin = stats.min;
+			fpdl->bounds.zmax = stats.max;
+		}
+		else if ( i == pdl->schema->m_position )
+		{
+			fpdl->bounds.mmin = stats.min;
+			fpdl->bounds.mmax = stats.max;
+		}
 
 		pc_point_set_double_by_index(&(fpdl->stats->min), i, stats.min);
 		pc_point_set_double_by_index(&(fpdl->stats->max), i, stats.max);

lib/pc_patch_dimensional.c

@@ -21,7 +21,7 @@ typedef struct
 	int8_t readonly;
 	const PCSCHEMA *schema;
 	uint32_t npoints;
-	double xmin, xmax, ymin, ymax;
+	double xmin, xmax, ymin, ymax, zmin, zmax, mmin, mmax;
 	PCSTATS *stats;
 	PCBYTES *bytes;
 } PCPATCH_DIMENSIONAL;
@@ -179,30 +179,54 @@ pc_patch_dimensional_free(PCPATCH_DIMENSIONAL *pdl)
 int
 pc_patch_dimensional_compute_extent(PCPATCH_DIMENSIONAL *pdl)
 {
-	double xmin, xmax, ymin, ymax, xavg, yavg;
+	double min, max, avg;
 	int rv;
 	PCBYTES *pcb;
 
 	assert(pdl);
 	assert(pdl->schema);
 
+	pc_bounds_init(&(pdl->bounds), pdl->schema);
+
 	/* Get x extremes */
-	pcb = &(pdl->bytes[pdl->schema->x_position]);
-	rv = pc_bytes_minmax(pcb, &xmin, &xmax, &xavg);
-	if ( PC_FAILURE == rv ) return PC_FAILURE;
-	xmin = pc_value_scale_offset(xmin, pdl->schema->dims[pdl->schema->x_position]);
-	xmax = pc_value_scale_offset(xmax, pdl->schema->dims[pdl->schema->x_position]);
-	pdl->bounds.xmin = xmin;
-	pdl->bounds.xmax = xmax;
+	if(pdl->schema->x_position!=-1)
+	{
+		pcb = &(pdl->bytes[pdl->schema->x_position]);
+		rv = pc_bytes_minmax(pcb, &min, &max, &avg);
+		if ( PC_FAILURE == rv ) return PC_FAILURE;
+		pdl->bounds.xmin = pc_value_scale_offset(min, pdl->schema->dims[pdl->schema->x_position]);
+		pdl->bounds.xmax = pc_value_scale_offset(max, pdl->schema->dims[pdl->schema->x_position]);
+	}
 
 	/* Get y extremes */
-	pcb = &(pdl->bytes[pdl->schema->y_position]);
-	rv = pc_bytes_minmax(pcb, &ymin, &ymax, &yavg);
-	if ( PC_FAILURE == rv ) return PC_FAILURE;
-	ymin = pc_value_scale_offset(ymin, pdl->schema->dims[pdl->schema->y_position]);
-	ymax = pc_value_scale_offset(ymax, pdl->schema->dims[pdl->schema->y_position]);
-	pdl->bounds.ymin = ymin;
-	pdl->bounds.ymax = ymax;
+	if(pdl->schema->y_position!=-1)
+	{
+		pcb = &(pdl->bytes[pdl->schema->y_position]);
+		rv = pc_bytes_minmax(pcb, &min, &max, &avg);
+		if ( PC_FAILURE == rv ) return PC_FAILURE;
+		pdl->bounds.ymin = pc_value_scale_offset(min, pdl->schema->dims[pdl->schema->y_position]);
+		pdl->bounds.ymax = pc_value_scale_offset(max, pdl->schema->dims[pdl->schema->y_position]);
+	}
+
+	/* Get z extremes */
+	if(pdl->schema->z_position!=-1)
+	{
+		pcb = &(pdl->bytes[pdl->schema->z_position]);
+		rv = pc_bytes_minmax(pcb, &min, &max, &avg);
+		if ( PC_FAILURE == rv ) return PC_FAILURE;
+		pdl->bounds.zmin = pc_value_scale_offset(min, pdl->schema->dims[pdl->schema->z_position]);
+		pdl->bounds.zmax = pc_value_scale_offset(max, pdl->schema->dims[pdl->schema->z_position]);
+	}
+
+	/* Get m extremes */
+	if(pdl->schema->m_position!=-1)
+	{
+		pcb = &(pdl->bytes[pdl->schema->m_position]);
+		rv = pc_bytes_minmax(pcb, &min, &max, &avg);
+		if ( PC_FAILURE == rv ) return PC_FAILURE;
+		pdl->bounds.mmin = pc_value_scale_offset(min, pdl->schema->dims[pdl->schema->m_position]);
+		pdl->bounds.mmax = pc_value_scale_offset(max, pdl->schema->dims[pdl->schema->m_position]);
+	}
 
 	return PC_SUCCESS;
 }

lib/pc_patch_ght.c

@@ -402,6 +402,12 @@ pc_patch_ght_compute_extent(PCPATCH_GHT *patch)
 	patch->bounds.ymin = area.y.min;
 	patch->bounds.ymax = area.y.max;
 
+	// provide a conservative non-discriminative bounding box
+	patch->bounds.zmin = -DBLMAX;
+	patch->bounds.zmax =  DBLMAX;
+	patch->bounds.mmin = -DBLMAX;
+	patch->bounds.mmax =  DBLMAX;
+
 	// ght_tree_free(tree);
 
 	return PC_SUCCESS;

lib/pc_patch_uncompressed.c

@@ -180,7 +180,7 @@ pc_patch_uncompressed_make(const PCSCHEMA *s, uint32_t maxpoints)
 	{
 		pch->data = pcalloc(datasize);
 	}
-	pc_bounds_init(&(pch->bounds));
+	pc_bounds_init(&(pch->bounds), s);
 
 	return pch;
 }
@@ -190,24 +190,17 @@ pc_patch_uncompressed_compute_extent(PCPATCH_UNCOMPRESSED *patch)
 {
 	int i;
 	PCPOINT *pt = pc_point_from_data(patch->schema, patch->data);
-	PCBOUNDS b;
-	double x, y;
+	double v;
 
 	/* Calculate bounds */
-	pc_bounds_init(&b);
+	pc_bounds_init(&(patch->bounds), patch->schema);
 	for ( i = 0; i < patch->npoints; i++ )
 	{
 		/* Just push the data buffer forward by one point at a time */
 		pt->data = patch->data + i * patch->schema->size;
-		x = pc_point_get_x(pt);
-		y = pc_point_get_y(pt);
-		if ( b.xmin > x ) b.xmin = x;
-		if ( b.ymin > y ) b.ymin = y;
-		if ( b.xmax < x ) b.xmax = x;
-		if ( b.ymax < y ) b.ymax = y;
+		pc_bounds_expand(&(patch->bounds),pt);
 	}
 
-	patch->bounds = b;
 	pcfree(pt);
 	return PC_SUCCESS;
 }
@@ -282,7 +275,7 @@ pc_patch_uncompressed_from_pointlist(const PCPOINTLIST *pl)
 	ptr = pch->data;
 
 	/* Initialize bounds */
-	pc_bounds_init(&(pch->bounds));
+	pc_bounds_init(&(pch->bounds), s);
 
 	/* Set up basic info */
 	pch->readonly = PC_FALSE;
@@ -418,12 +411,7 @@ pc_patch_uncompressed_add_point(PCPATCH_UNCOMPRESSED *c, const PCPOINT *p)
 	c->npoints += 1;
 
 	/* Update bounding box */
-	x = pc_point_get_x(p);
-	y = pc_point_get_y(p);
-	if ( c->bounds.xmin > x ) c->bounds.xmin = x;
-	if ( c->bounds.ymin > y ) c->bounds.ymin = y;
-	if ( c->bounds.xmax < x ) c->bounds.xmax = x;
-	if ( c->bounds.ymax < y ) c->bounds.ymax = y;
+	pc_bounds_expand(&(c->bounds),p);
 
 	return PC_SUCCESS;
 }

lib/pc_point.c

@@ -329,6 +329,9 @@ pc_point_to_geometry_wkb(const PCPOINT *pt, size_t *wkbsize)
 	uint32_t srid = pt->schema->srid;
 	double x, y, z, m;
 
+	if (pt->schema->x_position < 0 || pt->schema->y_position < 0 )
+		return NULL;
+
 	if ( srid != 0 )
 	{
 		wkbtype |= srid_mask;
@@ -379,7 +382,7 @@ pc_point_to_geometry_wkb(const PCPOINT *pt, size_t *wkbsize)
 
 	if ( pt->schema->m_position > -1 )
 	{
-		m = pc_point_get_z(pt);
+		m = pc_point_get_m(pt);
 		memcpy(ptr, &m, 8); /* M */
 		ptr += 8;
 	}