lwin_twkb.c

Go to the documentation of this file.

/**********************************************************************
 *
 * PostGIS - Spatial Types for PostgreSQL
 * http://postgis.net
 *
 * PostGIS is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 2 of the License, or
 * (at your option) any later version.
 *
 * PostGIS is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with PostGIS.  If not, see <http://www.gnu.org/licenses/>.
 *
 **********************************************************************
 *
 * Copyright (C) 2014 Nicklas Avén
 *
 **********************************************************************/
 
 
#include <math.h>
#include "liblwgeom_internal.h"
#include "lwgeom_log.h"
#include "varint.h"
 
#define TWKB_IN_MAXCOORDS 4
 
typedef struct
{
        /* Pointers to the bytes */
        const uint8_t *twkb; /* Points to start of TWKB */
        const uint8_t *twkb_end; /* Points to end of TWKB */
        const uint8_t *pos; /* Current read position */
 
        uint32_t check; /* Simple validity checks on geometries */
        uint32_t lwtype; /* Current type we are handling */
 
        uint8_t has_bbox;
        uint8_t has_size;
        uint8_t has_idlist;
        uint8_t has_z;
        uint8_t has_m;
        uint8_t is_empty;
 
        /* Precision factors to convert ints to double */
        double factor;
        double factor_z;
        double factor_m;
 
        uint64_t size;
 
        /* Info about current geometry */
        uint8_t magic_byte; /* the magic byte contain info about if twkb contain id, size info, bboxes and precision */
 
        int ndims; /* Number of dimensions */
 
        int64_t *coords; /* An array to keep delta values from 4 dimensions */
 
} twkb_parse_state;
 
 
LWGEOM* lwgeom_from_twkb_state(twkb_parse_state *s);
 
 
/**********************************************************************/
 
static inline void twkb_parse_state_advance(twkb_parse_state *s, size_t next)
{
        if( (s->pos + next) > s->twkb_end)
        {
                lwerror("%s: TWKB structure does not match expected size!", __func__);
                // lwnotice("TWKB structure does not match expected size!");
        }
 
        s->pos += next;
}
 
static inline int64_t twkb_parse_state_varint(twkb_parse_state *s)
{
        size_t size;
        int64_t val = varint_s64_decode(s->pos, s->twkb_end, &size);
        twkb_parse_state_advance(s, size);
        return val;
}
 
static inline uint64_t twkb_parse_state_uvarint(twkb_parse_state *s)
{
        size_t size;
        uint64_t val = varint_u64_decode(s->pos, s->twkb_end, &size);
        twkb_parse_state_advance(s, size);
        return val;
}
 
static inline double twkb_parse_state_double(twkb_parse_state *s, double factor)
{
        size_t size;
        int64_t val = varint_s64_decode(s->pos, s->twkb_end, &size);
        twkb_parse_state_advance(s, size);
        return val / factor;
}
 
static inline void twkb_parse_state_varint_skip(twkb_parse_state *s)
{
        size_t size = varint_size(s->pos, s->twkb_end);
 
        if ( ! size )
                lwerror("%s: no varint to skip", __func__);
 
        twkb_parse_state_advance(s, size);
        return;
}
 
 
 
static uint32_t lwtype_from_twkb_type(uint8_t twkb_type)
{
        switch (twkb_type)
        {
                case 1:
                        return POINTTYPE;
                case 2:
                        return LINETYPE;
                case 3:
                        return POLYGONTYPE;
                case 4:
                        return MULTIPOINTTYPE;
                case 5:
                        return MULTILINETYPE;
                case 6:
                        return MULTIPOLYGONTYPE;
                case 7:
                        return COLLECTIONTYPE;
 
                default: /* Error! */
                        lwerror("Unknown WKB type");
                        return 0;
        }
        return 0;
}
 
static uint8_t byte_from_twkb_state(twkb_parse_state *s)
{
        uint8_t val = *(s->pos);
        twkb_parse_state_advance(s, WKB_BYTE_SIZE);
        return val;
}
 
 
static POINTARRAY* ptarray_from_twkb_state(twkb_parse_state *s, uint32_t npoints)
{
        POINTARRAY *pa = NULL;
        uint32_t ndims = s->ndims;
        uint32_t i;
        double *dlist;
 
        LWDEBUG(2,"Entering ptarray_from_twkb_state");
        LWDEBUGF(4,"Pointarray has %d points", npoints);
 
        /* Empty! */
        if( npoints == 0 )
                return ptarray_construct_empty(s->has_z, s->has_m, 0);
 
        pa = ptarray_construct(s->has_z, s->has_m, npoints);
        dlist = (double*)(pa->serialized_pointlist);
        for( i = 0; i < npoints; i++ )
        {
                int j = 0;
                /* X */
                s->coords[j] += twkb_parse_state_varint(s);
                dlist[ndims*i + j] = s->coords[j] / s->factor;
                j++;
                /* Y */
                s->coords[j] += twkb_parse_state_varint(s);
                dlist[ndims*i + j] = s->coords[j] / s->factor;
                j++;
                /* Z */
                if ( s->has_z )
                {
                        s->coords[j] += twkb_parse_state_varint(s);
                        dlist[ndims*i + j] = s->coords[j] / s->factor_z;
                        j++;
                }
                /* M */
                if ( s->has_m )
                {
                        s->coords[j] += twkb_parse_state_varint(s);
                        dlist[ndims*i + j] = s->coords[j] / s->factor_m;
                        j++;
                }
        }
 
        return pa;
}
 
static LWPOINT* lwpoint_from_twkb_state(twkb_parse_state *s)
{
        static uint32_t npoints = 1;
        POINTARRAY *pa;
 
        LWDEBUG(2,"Entering lwpoint_from_twkb_state");
 
        if ( s->is_empty )
                return lwpoint_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        pa = ptarray_from_twkb_state(s, npoints);
        return lwpoint_construct(SRID_UNKNOWN, NULL, pa);
}
 
static LWLINE* lwline_from_twkb_state(twkb_parse_state *s)
{
        uint32_t npoints;
        POINTARRAY *pa;
 
        LWDEBUG(2,"Entering lwline_from_twkb_state");
 
        if ( s->is_empty )
                return lwline_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        /* Read number of points */
        npoints = twkb_parse_state_uvarint(s);
 
        if ( npoints == 0 )
                return lwline_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        /* Read coordinates */
        pa = ptarray_from_twkb_state(s, npoints);
 
        if( pa == NULL )
                return lwline_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        if( s->check & LW_PARSER_CHECK_MINPOINTS && pa->npoints < 2 )
        {
                lwerror("%s must have at least two points", lwtype_name(s->lwtype));
                return NULL;
        }
 
        return lwline_construct(SRID_UNKNOWN, NULL, pa);
}
 
static LWPOLY* lwpoly_from_twkb_state(twkb_parse_state *s)
{
        uint32_t nrings;
        uint32_t i;
        LWPOLY *poly;
 
        LWDEBUG(2,"Entering lwpoly_from_twkb_state");
 
        if ( s->is_empty )
                return lwpoly_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        /* Read number of rings */
        nrings = twkb_parse_state_uvarint(s);
 
        /* Start w/ empty polygon */
        poly = lwpoly_construct_empty(SRID_UNKNOWN, s->has_z, s->has_m);
 
        LWDEBUGF(4,"Polygon has %d rings", nrings);
 
        /* Empty polygon? */
        if( nrings == 0 )
                return poly;
 
        for( i = 0; i < nrings; i++ )
        {
                /* Ret number of points */
                uint32_t npoints = twkb_parse_state_uvarint(s);
                POINTARRAY *pa = ptarray_from_twkb_state(s, npoints);
 
                /* Skip empty rings */
                if( pa == NULL )
                        continue;
 
                /* Force first and last points to be the same. */
                if( ! ptarray_is_closed_2d(pa) )
                {
                        POINT4D pt;
                        getPoint4d_p(pa, 0, &pt);
                        ptarray_append_point(pa, &pt, LW_FALSE);
                }
 
                /* Check for at least four points. */
                if( s->check & LW_PARSER_CHECK_MINPOINTS && pa->npoints < 4 )
                {
                        LWDEBUGF(2, "%s must have at least four points in each ring", lwtype_name(s->lwtype));
                        lwerror("%s must have at least four points in each ring", lwtype_name(s->lwtype));
                        return NULL;
                }
 
                /* Add ring to polygon */
                if ( lwpoly_add_ring(poly, pa) == LW_FAILURE )
                {
                        LWDEBUG(2, "Unable to add ring to polygon");
                        lwerror("Unable to add ring to polygon");
                }
 
        }
        return poly;
}
 
 
static LWCOLLECTION* lwmultipoint_from_twkb_state(twkb_parse_state *s)
{
        int ngeoms, i;
        LWGEOM *geom = NULL;
        LWCOLLECTION *col = lwcollection_construct_empty(s->lwtype, SRID_UNKNOWN, s->has_z, s->has_m);
 
        LWDEBUG(2,"Entering lwmultipoint_from_twkb_state");
 
        if ( s->is_empty )
                return col;
 
        /* Read number of geometries */
        ngeoms = twkb_parse_state_uvarint(s);
        LWDEBUGF(4,"Number of geometries %d", ngeoms);
 
        /* It has an idlist, we need to skip that */
        if ( s->has_idlist )
        {
                for ( i = 0; i < ngeoms; i++ )
                        twkb_parse_state_varint_skip(s);
        }
 
        for ( i = 0; i < ngeoms; i++ )
        {
                geom = lwpoint_as_lwgeom(lwpoint_from_twkb_state(s));
                if ( lwcollection_add_lwgeom(col, geom) == NULL )
                {
                        lwerror("Unable to add geometry (%p) to collection (%p)", (void *) geom, (void *) col);
                        return NULL;
                }
        }
 
        return col;
}
 
static LWCOLLECTION* lwmultiline_from_twkb_state(twkb_parse_state *s)
{
        int ngeoms, i;
        LWGEOM *geom = NULL;
        LWCOLLECTION *col = lwcollection_construct_empty(s->lwtype, SRID_UNKNOWN, s->has_z, s->has_m);
 
        LWDEBUG(2,"Entering lwmultilinestring_from_twkb_state");
 
        if ( s->is_empty )
                return col;
 
        /* Read number of geometries */
        ngeoms = twkb_parse_state_uvarint(s);
 
        LWDEBUGF(4,"Number of geometries %d",ngeoms);
 
        /* It has an idlist, we need to skip that */
        if ( s->has_idlist )
        {
                for ( i = 0; i < ngeoms; i++ )
                        twkb_parse_state_varint_skip(s);
        }
 
        for ( i = 0; i < ngeoms; i++ )
        {
                geom = lwline_as_lwgeom(lwline_from_twkb_state(s));
                if ( lwcollection_add_lwgeom(col, geom) == NULL )
                {
                        lwerror("Unable to add geometry (%p) to collection (%p)", (void *) geom, (void *) col);
                        return NULL;
                }
        }
 
        return col;
}
 
static LWCOLLECTION* lwmultipoly_from_twkb_state(twkb_parse_state *s)
{
        int ngeoms, i;
        LWGEOM *geom = NULL;
        LWCOLLECTION *col = lwcollection_construct_empty(s->lwtype, SRID_UNKNOWN, s->has_z, s->has_m);
 
        LWDEBUG(2,"Entering lwmultipolygon_from_twkb_state");
 
        if ( s->is_empty )
                return col;
 
        /* Read number of geometries */
        ngeoms = twkb_parse_state_uvarint(s);
        LWDEBUGF(4,"Number of geometries %d",ngeoms);
 
        /* It has an idlist, we need to skip that */
        if ( s->has_idlist )
        {
                for ( i = 0; i < ngeoms; i++ )
                        twkb_parse_state_varint_skip(s);
        }
 
        for ( i = 0; i < ngeoms; i++ )
        {
                geom = lwpoly_as_lwgeom(lwpoly_from_twkb_state(s));
                if ( lwcollection_add_lwgeom(col, geom) == NULL )
                {
                        lwerror("Unable to add geometry (%p) to collection (%p)", (void *) geom, (void *) col);
                        return NULL;
                }
        }
 
        return col;
}
 
 
static LWCOLLECTION* lwcollection_from_twkb_state(twkb_parse_state *s)
{
        int ngeoms, i;
        LWGEOM *geom = NULL;
        LWCOLLECTION *col = lwcollection_construct_empty(s->lwtype, SRID_UNKNOWN, s->has_z, s->has_m);
 
        LWDEBUG(2,"Entering lwcollection_from_twkb_state");
 
        if ( s->is_empty )
                return col;
 
        /* Read number of geometries */
        ngeoms = twkb_parse_state_uvarint(s);
 
        LWDEBUGF(4,"Number of geometries %d",ngeoms);
 
        /* It has an idlist, we need to skip that */
        if ( s->has_idlist )
        {
                for ( i = 0; i < ngeoms; i++ )
                        twkb_parse_state_varint_skip(s);
        }
 
        for ( i = 0; i < ngeoms; i++ )
        {
                geom = lwgeom_from_twkb_state(s);
                if ( lwcollection_add_lwgeom(col, geom) == NULL )
                {
                        lwerror("Unable to add geometry (%p) to collection (%p)", (void *) geom, (void *) col);
                        return NULL;
                }
        }
 
 
        return col;
}
 
 
static void header_from_twkb_state(twkb_parse_state *s)
{
        LWDEBUG(2,"Entering magicbyte_from_twkb_state");
 
        uint8_t extended_dims;
 
        /* Read the first two bytes */
        uint8_t type_precision = byte_from_twkb_state(s);
        uint8_t metadata = byte_from_twkb_state(s);
 
        /* Strip type and precision out of first byte */
        uint8_t type = type_precision & 0x0F;
        int8_t precision = unzigzag8((type_precision & 0xF0) >> 4);
 
        /* Convert TWKB type to internal type */
        s->lwtype = lwtype_from_twkb_type(type);
 
        /* Convert the precision into factor */
        s->factor = pow(10, (double)precision);
 
        /* Strip metadata flags out of second byte */
        s->has_bbox   =  metadata & 0x01;
        s->has_size   = (metadata & 0x02) >> 1;
        s->has_idlist = (metadata & 0x04) >> 2;
        extended_dims = (metadata & 0x08) >> 3;
        s->is_empty   = (metadata & 0x10) >> 4;
 
        /* Flag for higher dims means read a third byte */
        if ( extended_dims )
        {
                int8_t precision_z, precision_m;
 
                extended_dims = byte_from_twkb_state(s);
 
                /* Strip Z/M presence and precision from ext byte */
                s->has_z    = (extended_dims & 0x01);
                s->has_m    = (extended_dims & 0x02) >> 1;
                precision_z = (extended_dims & 0x1C) >> 2;
                precision_m = (extended_dims & 0xE0) >> 5;
 
                /* Convert the precision into factor */
                s->factor_z = pow(10, (double)precision_z);
                s->factor_m = pow(10, (double)precision_m);
        }
        else
        {
                s->has_z = 0;
                s->has_m = 0;
                s->factor_z = 0;
                s->factor_m = 0;
        }
 
        /* Read the size, if there is one */
        if ( s->has_size )
        {
                s->size = twkb_parse_state_uvarint(s);
        }
 
        /* Calculate the number of dimensions */
        s->ndims = 2 + s->has_z + s->has_m;
 
        return;
}
 
 
 
LWGEOM* lwgeom_from_twkb_state(twkb_parse_state *s)
{
        GBOX bbox;
        LWGEOM *geom = NULL;
        uint32_t has_bbox = LW_FALSE;
        int i;
 
        /* Read the first two bytes, and optional */
        /* extended precision info and optional size info */
        header_from_twkb_state(s);
 
        /* Just experienced a geometry header, so now we */
        /* need to reset our coordinate deltas */
        for ( i = 0; i < TWKB_IN_MAXCOORDS; i++ )
        {
                s->coords[i] = 0.0;
        }
 
        /* Read the bounding box, is there is one */
        if ( s->has_bbox )
        {
                /* Initialize */
                has_bbox = s->has_bbox;
                memset(&bbox, 0, sizeof(GBOX));
                bbox.flags = lwflags(s->has_z, s->has_m, 0);
 
                /* X */
                bbox.xmin = twkb_parse_state_double(s, s->factor);
                bbox.xmax = bbox.xmin + twkb_parse_state_double(s, s->factor);
                /* Y */
                bbox.ymin = twkb_parse_state_double(s, s->factor);
                bbox.ymax = bbox.ymin + twkb_parse_state_double(s, s->factor);
                /* Z */
                if ( s->has_z )
                {
                        bbox.zmin = twkb_parse_state_double(s, s->factor_z);
                        bbox.zmax = bbox.zmin + twkb_parse_state_double(s, s->factor_z);
                }
                /* M */
                if ( s->has_m )
                {
                        bbox.mmin = twkb_parse_state_double(s, s->factor_m);
                        bbox.mmax = bbox.mmin + twkb_parse_state_double(s, s->factor_m);
                }
        }
 
        /* Switch to code for the particular type we're dealing with */
        switch( s->lwtype )
        {
                case POINTTYPE:
                        geom = lwpoint_as_lwgeom(lwpoint_from_twkb_state(s));
                        break;
                case LINETYPE:
                        geom = lwline_as_lwgeom(lwline_from_twkb_state(s));
                        break;
                case POLYGONTYPE:
                        geom = lwpoly_as_lwgeom(lwpoly_from_twkb_state(s));
                        break;
                case MULTIPOINTTYPE:
                        geom = lwcollection_as_lwgeom(lwmultipoint_from_twkb_state(s));
                        break;
                case MULTILINETYPE:
                        geom = lwcollection_as_lwgeom(lwmultiline_from_twkb_state(s));
                        break;
                case MULTIPOLYGONTYPE:
                        geom = lwcollection_as_lwgeom(lwmultipoly_from_twkb_state(s));
                        break;
                case COLLECTIONTYPE:
                        geom = lwcollection_as_lwgeom(lwcollection_from_twkb_state(s));
                        break;
                /* Unknown type! */
                default:
                        lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(s->lwtype));
                        break;
        }
 
        if ( has_bbox )
                geom->bbox = gbox_clone(&bbox);
 
        return geom;
}
 
 
LWGEOM* lwgeom_from_twkb(const uint8_t *twkb, size_t twkb_size, char check)
{
        int64_t coords[TWKB_IN_MAXCOORDS] = {0, 0, 0, 0};
        twkb_parse_state s;
 
        LWDEBUG(2,"Entering lwgeom_from_twkb");
        LWDEBUGF(4,"twkb_size: %d",(int) twkb_size);
 
        /* Zero out the state */
        memset(&s, 0, sizeof(twkb_parse_state));
 
        /* Initialize the state appropriately */
        s.twkb = s.pos = twkb;
        s.twkb_end = twkb + twkb_size;
        s.check = check;
        s.coords = coords;
 
        /* Read the rest of the geometry */
        return lwgeom_from_twkb_state(&s);
}