From 4f033847fea4176bde6585bb628f6d2de25b8fd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Lemoine?= Date: Wed, 29 Mar 2017 17:20:54 +0200 Subject: [PATCH 1/3] Add a PC_BoundingDiagonalAsBinary(p PCPATCH) function --- README.md | 12 +++ lib/cunit/CMakeLists.txt | 1 + lib/cunit/Makefile | 3 +- lib/cunit/cu_pc_util.c | 109 ++++++++++++++++++++++++ lib/cunit/cu_tester.c | 2 + lib/pc_api.h | 9 ++ lib/pc_api_internal.h | 9 ++ lib/pc_stats.c | 2 +- lib/pc_util.c | 154 ++++++++++++++++++++++++++++++++++ pgsql/expected/pointcloud.out | 11 +++ pgsql/pc_inout.c | 45 ++++++++++ pgsql/pointcloud.sql.in | 4 + pgsql/sql/pointcloud.sql | 7 ++ 13 files changed, 366 insertions(+), 2 deletions(-) create mode 100644 lib/cunit/cu_pc_util.c diff --git a/README.md b/README.md index f35ea19..c6a49fa 100644 --- a/README.md +++ b/README.md @@ -492,6 +492,18 @@ Now that you have created two tables, you'll see entries for them in the `pointc > a3703dba5fc0ec51b81e858b46400ad7a3703dba5fc0e17a > 14ae4781464090c2f5285cbf5fc0e17a14ae47814640 +**PC_BoundingDiagonalAsBinary(p pcpatch)** returns **bytea** + +> Return the OGC "well-known binary" format for the bounding diagonal of the patch. +> +> SELECT PC_BoundingDiagonalAsBinary( +> PC_Patch(ARRAY[ +> PC_MakePoint(1, ARRAY[0.,0.,0.,10.]), +> PC_MakePoint(1, ARRAY[1.,1.,1.,10.]), +> PC_MakePoint(1, ARRAY[10.,10.,10.,10.])])); +> +> \x01020000a0e610000002000000000000000000000000000000000000000000000000000000000000000000244000000000000024400000000000002440 + ## PostGIS Integration ## The `pointcloud_postgis` extension adds functions that allow you to use PostgreSQL Pointcloud with PostGIS, converting PcPoint and PcPatch to Geometry and doing spatial filtering on point cloud data. The `pointcloud_postgis` extension depends on both the `postgis` and `pointcloud` extensions, so they must be installed first: diff --git a/lib/cunit/CMakeLists.txt b/lib/cunit/CMakeLists.txt index 66541f6..f9c2fca 100644 --- a/lib/cunit/CMakeLists.txt +++ b/lib/cunit/CMakeLists.txt @@ -10,6 +10,7 @@ set (PC_TEST_SOURCES cu_pc_patch.c cu_pc_patch_lazperf.c cu_pc_sort.c + cu_pc_util.c cu_tester.c ) diff --git a/lib/cunit/Makefile b/lib/cunit/Makefile index 0222317..2b5307a 100644 --- a/lib/cunit/Makefile +++ b/lib/cunit/Makefile @@ -15,7 +15,8 @@ OBJS = \ cu_pc_patch.o \ cu_pc_patch_ght.o \ cu_pc_patch_lazperf.o \ - cu_pc_sort.o + cu_pc_sort.o \ + cu_pc_util.o ifeq ($(CUNIT_LDFLAGS),) # No cunit? Emit message and continue diff --git a/lib/cunit/cu_pc_util.c b/lib/cunit/cu_pc_util.c new file mode 100644 index 0000000..16038dd --- /dev/null +++ b/lib/cunit/cu_pc_util.c @@ -0,0 +1,109 @@ +/*********************************************************************** +* cu_pc_util.c +* +* Testing for the util functions +* +* Portions Copyright (c) 2017, Oslandia +* +***********************************************************************/ + +#include "CUnit/Basic.h" +#include "cu_tester.h" + + +/* GLOBALS ************************************************************/ + +static PCSCHEMA *schema = NULL; +static const char *xmlfile = "data/pdal-schema.xml"; + +/* Setup/teardown for this suite */ +static int +init_suite(void) +{ + char *xmlstr = file_to_str(xmlfile); + schema = pc_schema_from_xml(xmlstr); + pcfree(xmlstr); + if ( !schema ) return 1; + + return 0; +} + +static int +clean_suite(void) +{ + pc_schema_free(schema); + return 0; +} + + +/* TESTS **************************************************************/ + +static void +test_bounding_diagonal_wkb_from_bounds() +{ + PCBOUNDS bounds; + size_t wkbsize; + uint8_t *wkb; + char *wkbhex; + + bounds.xmin = -10; + bounds.xmax = 10; + bounds.ymin = -10; + bounds.ymax = 10; + + wkb = pc_bounding_diagonal_wkb_from_bounds(&bounds, schema, &wkbsize); + CU_ASSERT(wkb != NULL); + CU_ASSERT(wkbsize == 41); + + wkbhex = hexbytes_from_bytes(wkb, wkbsize); + CU_ASSERT(wkbhex != NULL); + CU_ASSERT_STRING_EQUAL(wkbhex, "01020000000200000000000000000024C000000000000024C000000000000024400000000000002440"); + + pcfree(wkb); + pcfree(wkbhex); +} + +static void +test_bounding_diagonal_wkb_from_stats() +{ + PCSTATS *stats; + size_t wkbsize; + uint8_t *wkb; + char *wkbhex; + + stats = pc_stats_new(schema); + + pc_point_set_x(&stats->min, -10); + pc_point_set_x(&stats->max, 10); + pc_point_set_y(&stats->min, -10); + pc_point_set_y(&stats->max, 10); + pc_point_set_z(&stats->min, -10); + pc_point_set_z(&stats->max, 10); + + wkb = pc_bounding_diagonal_wkb_from_stats(stats, &wkbsize); + CU_ASSERT(wkb != NULL); + CU_ASSERT(wkbsize == 73); + + wkbhex = hexbytes_from_bytes(wkb, wkbsize); + CU_ASSERT(wkbhex != NULL); + CU_ASSERT_STRING_EQUAL(wkbhex, "01020000C00200000000000000000024C000000000000024C000000000000024C000000000000000000000000000002440000000000000244000000000000024400000000000000000"); + + pc_stats_free(stats); + pcfree(wkb); + pcfree(wkbhex); +} + +/* REGISTER ***********************************************************/ + +CU_TestInfo util_tests[] = { + PC_TEST(test_bounding_diagonal_wkb_from_bounds), + PC_TEST(test_bounding_diagonal_wkb_from_stats), + CU_TEST_INFO_NULL +}; + +CU_SuiteInfo util_suite = { + .pName = "util", + .pInitFunc = init_suite, + .pCleanupFunc = clean_suite, + .pTests = util_tests +}; diff --git a/lib/cunit/cu_tester.c b/lib/cunit/cu_tester.c index bee4167..306a78f 100644 --- a/lib/cunit/cu_tester.c +++ b/lib/cunit/cu_tester.c @@ -20,6 +20,7 @@ extern CU_SuiteInfo ght_suite; extern CU_SuiteInfo bytes_suite; extern CU_SuiteInfo lazperf_suite; extern CU_SuiteInfo sort_suite; +extern CU_SuiteInfo util_suite; /** * CUnit error handler @@ -56,6 +57,7 @@ int main(int argc, char *argv[]) bytes_suite, lazperf_suite, sort_suite, + util_suite, CU_SUITE_INFO_NULL }; diff --git a/lib/pc_api.h b/lib/pc_api.h index 94c3d1a..1d38218 100644 --- a/lib/pc_api.h +++ b/lib/pc_api.h @@ -431,6 +431,9 @@ int pc_bytes_deserialize(const uint8_t *buf, const PCDIMENSION *dim, PCBYTES *pc /** Wrap serialized stats in a new stats objects */ PCSTATS* pc_stats_new_from_data(const PCSCHEMA *schema, const uint8_t *mindata, const uint8_t *maxdata, const uint8_t *avgdata); +/** Allocate a stats object */ +PCSTATS* pc_stats_new(const PCSCHEMA *schema); + /** Free a stats object */ void pc_stats_free(PCSTATS *stats); @@ -446,6 +449,12 @@ int pc_patch_compute_extent(PCPATCH *patch); /** True/false if bounds intersect */ int pc_bounds_intersects(const PCBOUNDS *b1, const PCBOUNDS *b2); +/** Returns OGC WKB of the bounding diagonal of XY bounds */ +uint8_t* pc_bounding_diagonal_wkb_from_bounds(const PCBOUNDS *bounds, const PCSCHEMA *schema, size_t *wkbsize); + +/** Returns OGC WKB of the bounding diagonal of XY, XYZ, XYM or XYZM bounds */ +uint8_t* pc_bounding_diagonal_wkb_from_stats(const PCSTATS *stats, size_t *wkbsize); + /** Subset batch based on less-than condition on dimension */ PCPATCH* pc_patch_filter_lt_by_name(const PCPATCH *pa, const char *name, double val); diff --git a/lib/pc_api_internal.h b/lib/pc_api_internal.h index 21b5433..f23dacc 100644 --- a/lib/pc_api_internal.h +++ b/lib/pc_api_internal.h @@ -106,6 +106,15 @@ int16_t wkb_get_int16(const uint8_t *wkb, int flip_endian); /** Read the number of points from a wkb */ uint32_t wkb_get_npoints(const uint8_t *wkb); +/** Write a double into a byte array */ +uint8_t *wkb_set_double(uint8_t *wkb, double d); + +/** Write a uint32 into a byte array */ +uint8_t *wkb_set_uint32(uint8_t *wkb, uint32_t i); + +/** Write a char into a byte array */ +uint8_t *wkb_set_char(uint8_t *wkb, char c); + /** Force a byte array into the machine endianness */ uint8_t* uncompressed_bytes_flip_endian(const uint8_t *bytebuf, const PCSCHEMA *schema, uint32_t npoints); diff --git a/lib/pc_stats.c b/lib/pc_stats.c index ed01d37..1bc827f 100644 --- a/lib/pc_stats.c +++ b/lib/pc_stats.c @@ -89,7 +89,7 @@ pc_stats_new_from_data(const PCSCHEMA *schema, const uint8_t *mindata, const uin * point shells and the data areas underneath. Used for initial calcution * of patch stats, when objects first created. */ -static PCSTATS * +PCSTATS * pc_stats_new(const PCSCHEMA *schema) { /*size_t sz = schema->size;*/ diff --git a/lib/pc_util.c b/lib/pc_util.c index 1519e8a..6d803c5 100644 --- a/lib/pc_util.c +++ b/lib/pc_util.c @@ -156,6 +156,30 @@ wkb_get_int16(const uint8_t *wkb, int flip_endian) return i; } +uint8_t * +wkb_set_double(uint8_t *wkb, double d) +{ + memcpy(wkb, &d, 8); + wkb += 8; + return wkb; +} + +uint8_t * +wkb_set_uint32(uint8_t *wkb, uint32_t i) +{ + memcpy(wkb, &i, 4); + wkb += 4; + return wkb; +} + +uint8_t * +wkb_set_char(uint8_t *wkb, char c) +{ + memcpy(wkb, &c, 1); + wkb += 1; + return wkb; +} + uint32_t wkb_get_pcid(const uint8_t *wkb) { @@ -206,6 +230,7 @@ wkb_get_npoints(const uint8_t *wkb) } return npoints; } + uint8_t* uncompressed_bytes_flip_endian(const uint8_t *bytebuf, const PCSCHEMA *schema, uint32_t npoints) { @@ -263,3 +288,132 @@ void pc_bounds_merge(PCBOUNDS *b1, const PCBOUNDS *b2) if ( b2->ymax > b1->ymax ) b1->ymax = b2->ymax; } +static uint32_t srid_mask = 0x20000000; +static uint32_t m_mask = 0x40000000; +static uint32_t z_mask = 0x80000000; + +uint8_t * +pc_bounding_diagonal_wkb_from_bounds( + const PCBOUNDS *bounds, const PCSCHEMA *schema, size_t *wkbsize) +{ + uint8_t *wkb, *ptr; + uint32_t wkbtype; + size_t size; + + wkbtype = 2; /* WKB LINESTRING */ + size = 1 + 4 + 4 + (2 * 2 * 8); /* endian + type + npoints + 2 dbl pts */ + + if ( schema->srid != 0 ) + { + wkbtype |= srid_mask; + size += 4; + } + + wkb = pcalloc(size); + ptr = wkb; + + ptr = wkb_set_char(ptr, machine_endian()); /* Endian flag */ + ptr = wkb_set_uint32(ptr, wkbtype); /* Geometry type */ + + if ( schema->srid != 0 ) + { + ptr = wkb_set_uint32(ptr, schema->srid); /* SRID */ + } + + ptr = wkb_set_uint32(ptr, 2); /* NPOINTS = 2 */ + + // point 1 + ptr = wkb_set_double(ptr, bounds->xmin); + ptr = wkb_set_double(ptr, bounds->ymin); + + // point 2 + ptr = wkb_set_double(ptr, bounds->xmax); + ptr = wkb_set_double(ptr, bounds->ymax); + + if ( wkbsize ) + *wkbsize = size; + + return wkb; +} + +uint8_t * +pc_bounding_diagonal_wkb_from_stats(const PCSTATS *stats, size_t *wkbsize) +{ + const PCSCHEMA *schema = stats->min.schema; + const PCPOINT *stat; + uint8_t *wkb, *ptr; + uint32_t wkbtype; + size_t size; + double val; + + wkbtype = 2; /* WKB LINESTRING */ + size = 1 + 4 + 4 + (2 * 2 * 8); /* endian + type + npoints + 2 dbl pts */ + + if ( schema->srid != 0 ) + { + wkbtype |= srid_mask; + size += 4; + } + if ( schema->zdim ) + { + wkbtype |= z_mask; + size += 2 * 8; + } + if ( schema->mdim ) + { + wkbtype |= m_mask; + size += 2 * 8; + } + + wkb = pcalloc(size); + ptr = wkb; + + ptr = wkb_set_char(ptr, machine_endian()); /* Endian flag */ + ptr = wkb_set_uint32(ptr, wkbtype); /* Geometry type */ + + if ( schema->srid != 0 ) + { + ptr = wkb_set_uint32(ptr, schema->srid); /* SRID */ + } + + ptr = wkb_set_uint32(ptr, 2); /* NPOINTS = 2 */ + + // point 1 + stat = &stats->min; + pc_point_get_x(stat, &val); + ptr = wkb_set_double(ptr, val); + pc_point_get_y(stat, &val); + ptr = wkb_set_double(ptr, val); + if ( schema->zdim ) + { + pc_point_get_z(stat, &val); + ptr = wkb_set_double(ptr, val); + } + if ( schema->mdim ) + { + pc_point_get_m(stat, &val); + ptr = wkb_set_double(ptr, val); + } + + // point 2 + stat = &stats->max; + pc_point_get_x(stat, &val); + ptr = wkb_set_double(ptr, val); + pc_point_get_y(stat, &val); + ptr = wkb_set_double(ptr, val); + if ( schema->zdim ) + { + pc_point_get_z(stat, &val); + ptr = wkb_set_double(ptr, val); + } + if ( schema->mdim ) + { + pc_point_get_m(stat, &val); + ptr = wkb_set_double(ptr, val); + } + + if ( wkbsize ) + *wkbsize = size; + + return wkb; +} diff --git a/pgsql/expected/pointcloud.out b/pgsql/expected/pointcloud.out index 3b5acfe..1d32e3e 100644 --- a/pgsql/expected/pointcloud.out +++ b/pgsql/expected/pointcloud.out @@ -601,4 +601,15 @@ FROM ( SELECT #78 | -1 | -1 | 0 | 0 | 4862413 | 4862413 | 1 | 1 (1 row) +-- test for PC_BoundingDiagonalAsBinary +SELECT PC_BoundingDiagonalAsBinary( + PC_Patch(ARRAY[ + PC_MakePoint(1, ARRAY[0.,0.,0.,10.]), + PC_MakePoint(1, ARRAY[1.,1.,1.,10.]), + PC_MakePoint(1, ARRAY[10.,10.,10.,10.])])); + pc_boundingdiagonalasbinary +---------------------------------------------------------------------------------------------------------------------- + \x010200008002000000000000000000000000000000000000000000000000000000000000000000244000000000000024400000000000002440 +(1 row) + TRUNCATE pointcloud_formats; diff --git a/pgsql/pc_inout.c b/pgsql/pc_inout.c index fd734cc..a6e1be3 100644 --- a/pgsql/pc_inout.c +++ b/pgsql/pc_inout.c @@ -32,6 +32,7 @@ Datum pcpoint_as_text(PG_FUNCTION_ARGS); Datum pcpatch_as_text(PG_FUNCTION_ARGS); Datum pcpoint_as_bytea(PG_FUNCTION_ARGS); Datum pcpatch_bytea_envelope(PG_FUNCTION_ARGS); +Datum pcpatch_bounding_diagonal_as_bytea(PG_FUNCTION_ARGS); static void @@ -313,6 +314,50 @@ Datum pcpatch_envelope_as_bytea(PG_FUNCTION_ARGS) PG_RETURN_BYTEA_P(wkb); } +PG_FUNCTION_INFO_V1(pcpatch_bounding_diagonal_as_bytea); +Datum pcpatch_bounding_diagonal_as_bytea(PG_FUNCTION_ARGS) +{ + static const size_t stats_size_guess = 400; + SERIALIZED_PATCH *serpatch; + PCSCHEMA *schema; + PCSTATS *stats; + uint8 *bytes; + size_t bytes_size; + bytea *wkb; + size_t wkb_size; + + serpatch = PG_GETHEADERX_SERPATCH_P(0, stats_size_guess); + schema = pc_schema_from_pcid(serpatch->pcid, fcinfo); + + if ( schema->zdim || schema->mdim ) + { + if ( stats_size_guess < pc_stats_size(schema) ) + serpatch = PG_GETHEADERX_SERPATCH_P(0, pc_stats_size(schema)); + + stats = pc_patch_stats_deserialize(schema, serpatch->data); + if ( !stats ) + PG_RETURN_NULL(); + + bytes = pc_bounding_diagonal_wkb_from_stats(stats, &bytes_size); + + pc_stats_free(stats); + } + else + { + bytes = pc_bounding_diagonal_wkb_from_bounds( + &serpatch->bounds, schema, &bytes_size); + } + + wkb_size = VARHDRSZ + bytes_size; + wkb = palloc(wkb_size); + memcpy(VARDATA(wkb), bytes, bytes_size); + SET_VARSIZE(wkb, wkb_size); + + pcfree(bytes); + + PG_RETURN_BYTEA_P(wkb); +} + PG_FUNCTION_INFO_V1(pc_typmod_in); Datum pc_typmod_in(PG_FUNCTION_ARGS) { diff --git a/pgsql/pointcloud.sql.in b/pgsql/pointcloud.sql.in index 5909482..6bd4748 100644 --- a/pgsql/pointcloud.sql.in +++ b/pgsql/pointcloud.sql.in @@ -313,6 +313,10 @@ CREATE OR REPLACE FUNCTION PC_Range(p pcpatch, first int4, count int4) RETURNS pcpatch AS 'MODULE_PATHNAME', 'pcpatch_range' LANGUAGE 'c' IMMUTABLE STRICT; +CREATE OR REPLACE FUNCTION PC_BoundingDiagonalAsBinary(p pcpatch) + RETURNS bytea AS 'MODULE_PATHNAME', 'pcpatch_bounding_diagonal_as_bytea' + LANGUAGE 'c' IMMUTABLE STRICT; + ------------------------------------------------------------------- -- POINTCLOUD_COLUMNS ------------------------------------------------------------------- diff --git a/pgsql/sql/pointcloud.sql b/pgsql/sql/pointcloud.sql index 6d38e4d..ecfe477 100644 --- a/pgsql/sql/pointcloud.sql +++ b/pgsql/sql/pointcloud.sql @@ -352,4 +352,11 @@ FROM ( SELECT 'y',0) p ) foo; +-- test for PC_BoundingDiagonalAsBinary +SELECT PC_BoundingDiagonalAsBinary( + PC_Patch(ARRAY[ + PC_MakePoint(1, ARRAY[0.,0.,0.,10.]), + PC_MakePoint(1, ARRAY[1.,1.,1.,10.]), + PC_MakePoint(1, ARRAY[10.,10.,10.,10.])])); + TRUNCATE pointcloud_formats; From 7dd0bccfee96bd9052158c898bdedf4c08f6e47a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Lemoine?= Date: Mon, 3 Apr 2017 10:15:39 +0200 Subject: [PATCH 2/3] Add a PC_BoundingDiagonal(p pcpatch) function --- README.md | 23 +++++++++++++++++++++++ pgsql_postgis/pointcloud_postgis--1.0.sql | 10 ++++++++++ 2 files changed, 33 insertions(+) diff --git a/README.md b/README.md index c6a49fa..82343ad 100644 --- a/README.md +++ b/README.md @@ -559,6 +559,29 @@ The `pointcloud_postgis` extension adds functions that allow you to use PostgreS > > POLYGON((-126.99 45.01,-126.99 45.09,-126.91 45.09,-126.91 45.01,-126.99 45.01)) +**PC_BoundingDiagonal(pcpatch)** returns **geometry** + +> Returns the bounding diagonal of a patch. This is a LineString (2D), a LineString Z or a LineString M or a LineString ZM, based on the existence of the Z and M dimensions in the patch. This function is useful for creating an index on a patch column. +> +> SELECT ST_AsText(PC_BoundingDiagonal(pa)) FROM patches; +> st_astext +> ------------------------------------------------ +> LINESTRING Z (-126.99 45.01 1,-126.91 45.09 9) +> LINESTRING Z (-126 46 100,-126 46 100) +> LINESTRING Z (-126.2 45.8 80,-126.11 45.89 89) +> LINESTRING Z (-126.4 45.6 60,-126.31 45.69 69) +> LINESTRING Z (-126.3 45.7 70,-126.21 45.79 79) +> LINESTRING Z (-126.8 45.2 20,-126.71 45.29 29) +> LINESTRING Z (-126.5 45.5 50,-126.41 45.59 59) +> LINESTRING Z (-126.6 45.4 40,-126.51 45.49 49) +> LINESTRING Z (-126.9 45.1 10,-126.81 45.19 19) +> LINESTRING Z (-126.7 45.3 30,-126.61 45.39 39) +> LINESTRING Z (-126.1 45.9 90,-126.01 45.99 99) +> (11 rows)gq +> +> For example, this is how one may want to create an index: +> +> CREATE INDEX ON patches USING GIST(PC_BoundingDiagonal(patch) gist_geometry_ops_nd); ## Compressions ## diff --git a/pgsql_postgis/pointcloud_postgis--1.0.sql b/pgsql_postgis/pointcloud_postgis--1.0.sql index e37f2a0..ea96725 100644 --- a/pgsql_postgis/pointcloud_postgis--1.0.sql +++ b/pgsql_postgis/pointcloud_postgis--1.0.sql @@ -54,3 +54,13 @@ CREATE OR REPLACE FUNCTION PC_Intersects(geometry, pcpatch) SELECT PC_Intersects($2, $1) $$ LANGUAGE 'sql'; + +----------------------------------------------------------------------------- +-- Function from pcpatch to LineString +-- +CREATE OR REPLACE FUNCTION PC_BoundingDiagonal(pcpatch) + RETURNS geometry AS + $$ + SELECT ST_GeomFromEWKB(PC_BoundingDiagonalAsBinary($1)) + $$ + LANGUAGE 'sql'; From eb0c186758d603953569286802635d6918894d8a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Lemoine?= Date: Mon, 3 Apr 2017 10:15:49 +0200 Subject: [PATCH 3/3] Update the NEWS file --- NEWS | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS b/NEWS index a2a0272..8a4cf25 100644 --- a/NEWS +++ b/NEWS @@ -11,6 +11,7 @@ - PC_Sort(pcpatch,text[]) (#106) - PC_IsSorted(pcpatch,text[],boolean) (#106) - PC_Range(pcpatch, int, int) returns pcpatch (#152) + - PC_BoundingDiagonalAsBinary(pcpatch) and PC_BoundingDiagonal(pcpach) (#158) - Enhancements - Support sigbits encoding for 64bit integers (#61) - Warn about truncated values (#68)