Merge pull request #158 from elemoine/pc-boundingdiagonal

Add PC_BoundingDiagonal(p PCPATCH)
This commit is contained in:
Éric Lemoine 2017-04-05 10:26:44 +02:00 committed by GitHub
commit 2c49df040f
15 changed files with 400 additions and 2 deletions

1
NEWS
View File

@ -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)

View File

@ -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:
@ -547,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 ##

View File

@ -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
)

View File

@ -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

109
lib/cunit/cu_pc_util.c Normal file
View File

@ -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
};

View File

@ -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
};

View File

@ -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);

View File

@ -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);

View File

@ -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;*/

View File

@ -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;
}

View File

@ -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;

View File

@ -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)
{

View File

@ -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
-------------------------------------------------------------------

View File

@ -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;

View File

@ -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';