lwout_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) 2013 Nicklas Avén
 *
 **********************************************************************/
 
 
#include "lwout_twkb.h"
 
/*
* GeometryType, and dimensions
*/
static uint8_t lwgeom_twkb_type(const LWGEOM *geom)
{
        uint8_t twkb_type = 0;
 
        LWDEBUG(2, "Entered lwgeom_twkb_type");
 
        switch ( geom->type )
        {
                case POINTTYPE:
                        twkb_type = WKB_POINT_TYPE;
                        break;
                case LINETYPE:
                        twkb_type = WKB_LINESTRING_TYPE;
                        break;
                case TRIANGLETYPE:
                case POLYGONTYPE:
                        twkb_type = WKB_POLYGON_TYPE;
                        break;
                case MULTIPOINTTYPE:
                        twkb_type = WKB_MULTIPOINT_TYPE;
                        break;
                case MULTILINETYPE:
                        twkb_type = WKB_MULTILINESTRING_TYPE;
                        break;
                case MULTIPOLYGONTYPE:
                        twkb_type = WKB_MULTIPOLYGON_TYPE;
                        break;
                case TINTYPE:
                case COLLECTIONTYPE:
                        twkb_type = WKB_GEOMETRYCOLLECTION_TYPE;
                        break;
                default:
                        lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(geom->type));
        }
        return twkb_type;
}
 
 
static size_t sizeof_bbox(TWKB_STATE *ts, int ndims)
{
        int i;
        uint8_t buf[16];
        size_t size = 0;
        LWDEBUGF(2, "Entered %s", __func__);
        for ( i = 0; i < ndims; i++ )
        {
                size += varint_s64_encode_buf(ts->bbox_min[i], buf);
                size += varint_s64_encode_buf((ts->bbox_max[i] - ts->bbox_min[i]), buf);
        }
        return size;
}
static void write_bbox(TWKB_STATE *ts, int ndims)
{
        int i;
        LWDEBUGF(2, "Entered %s", __func__);
        for ( i = 0; i < ndims; i++ )
        {
                bytebuffer_append_varint(ts->header_buf, ts->bbox_min[i]);
                bytebuffer_append_varint(ts->header_buf, (ts->bbox_max[i] - ts->bbox_min[i]));
        }
}
 
 
static int ptarray_to_twkb_buf(const POINTARRAY *pa, TWKB_GLOBALS *globals, TWKB_STATE *ts, int register_npoints, uint32_t minpoints)
{
        uint32_t ndims = FLAGS_NDIMS(pa->flags);
        uint32_t i, j;
        bytebuffer_t b;
        bytebuffer_t *b_p;
        int64_t nextdelta[MAX_N_DIMS];
        int npoints = 0;
        size_t npoints_offset = 0;
        uint32_t max_points_left = pa->npoints;
 
        LWDEBUGF(2, "Entered %s", __func__);
 
        /* Dispense with the empty case right away */
        if ( pa->npoints == 0 && register_npoints )
        {
                LWDEBUGF(4, "Register npoints:%d", pa->npoints);
                bytebuffer_append_uvarint(ts->geom_buf, pa->npoints);
                return 0;
        }
 
        /* If npoints is more than 127 it is unpredictable how many bytes npoints will need */
        /* Then we have to store the deltas in a temp buffer to later add them after npoints */
        /* If noints is below 128 we know 1 byte will be needed */
        /* Then we can make room for that 1 byte at once and write to */
        /* ordinary buffer */
        if( pa->npoints > 127 )
        {
                /* Independent buffer to hold the coordinates, so we can put the npoints */
                /* into the stream once we know how many points we actually have */
                bytebuffer_init_with_size(&b, (size_t)3 * ndims * pa->npoints);
                b_p = &b;
        }
        else
        {
                /* We give an alias to our ordinary buffer */
                b_p = ts->geom_buf;
                if ( register_npoints )
                {
                        /* We do not store a pointer to the place where we want the npoints value */
                        /* Instead we store how far from the beginning of the buffer we want the value */
                        /* That is because we otherwise will get in trouble if the buffer is reallocated */
                        npoints_offset = b_p->writecursor - b_p->buf_start;
 
                        /* We just move the cursor 1 step to make room for npoints byte */
                        /* We use the function append_byte even if we have no value yet, */
                        /* since that gives us the check for big enough buffer and moves the cursor */
                        bytebuffer_append_byte(b_p, 0);
                }
        }
 
        for ( i = 0; i < pa->npoints; i++ )
        {
                double *dbl_ptr = (double*)getPoint_internal(pa, i);
                int64_t diff = 0;
 
                /* Write this coordinate to the buffer as a varint */
                for ( j = 0; j < ndims; j++ )
                {
                        /* To get the relative coordinate we don't get the distance */
                        /* from the last point but instead the distance from our */
                        /* last accumulated point. This is important to not build up an */
                        /* accumulated error when rounding the coordinates */
                        nextdelta[j] = (int64_t) llround(globals->factor[j] * dbl_ptr[j]) - ts->accum_rels[j];
                        LWDEBUGF(4, "deltavalue: %lld, ", nextdelta[j]);
                        diff += llabs(nextdelta[j]);
                }
 
                /* Skipping the first point is not allowed */
                /* If the sum(abs()) of all the deltas was zero, */
                /* then this was a duplicate point, so we can ignore it */
                if ( i > 0 && diff == 0 &&  max_points_left > minpoints )
                {
                        max_points_left--;
                        continue;
                }
 
                /* We really added a point, so... */
                npoints++;
 
                /* Write this vertex to the temporary buffer as varints */
                for ( j = 0; j < ndims; j++ )
                {
                        ts->accum_rels[j] += nextdelta[j];
                        bytebuffer_append_varint(b_p, nextdelta[j]);
                }
 
                /* See if this coordinate expands the bounding box */
                if( globals->variant & TWKB_BBOX )
                {
                        for ( j = 0; j < ndims; j++ )
                        {
                                if( ts->accum_rels[j] > ts->bbox_max[j] )
                                        ts->bbox_max[j] = ts->accum_rels[j];
 
                                if( ts->accum_rels[j] < ts->bbox_min[j] )
                                        ts->bbox_min[j] = ts->accum_rels[j];
                        }
                }
 
        }
 
        if ( pa->npoints > 127 )
        {
                /* Now write the temporary results into the main buffer */
                /* First the npoints */
                if ( register_npoints )
                        bytebuffer_append_uvarint(ts->geom_buf, npoints);
                /* Now the coordinates */
                bytebuffer_append_bytebuffer(ts->geom_buf, b_p);
 
                /* Clear our temporary buffer */
                bytebuffer_destroy_buffer(&b);
        }
        else
        {
                /* If we didn't use a temp buffer, we just write that npoints value */
                /* to where it belongs*/
                if ( register_npoints )
                        varint_u64_encode_buf(npoints, b_p->buf_start + npoints_offset);
        }
 
        return 0;
}
 
/******************************************************************
* POINTS
*******************************************************************/
 
static int lwpoint_to_twkb_buf(const LWPOINT *pt, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        LWDEBUGF(2, "Entered %s", __func__);
 
        /* Set the coordinates (don't write npoints) */
        ptarray_to_twkb_buf(pt->point, globals, ts, 0, 1);
        return 0;
}
 
/******************************************************************
* LINESTRINGS
*******************************************************************/
 
static int lwline_to_twkb_buf(const LWLINE *line, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        LWDEBUGF(2, "Entered %s", __func__);
 
        /* Set the coordinates (do write npoints) */
        ptarray_to_twkb_buf(line->points, globals, ts, 1, 2);
        return 0;
}
 
static int
lwtriangle_to_twkb_buf(const LWTRIANGLE *tri, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        LWDEBUGF(2, "Entered %s", __func__);
        bytebuffer_append_uvarint(ts->geom_buf, (uint64_t)1);
 
        /* Set the coordinates (do write npoints) */
        ptarray_to_twkb_buf(tri->points, globals, ts, 1, 2);
        return 0;
}
 
/******************************************************************
* POLYGONS
*******************************************************************/
 
static int lwpoly_to_twkb_buf(const LWPOLY *poly, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        uint32_t i;
 
        /* Set the number of rings */
        bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) poly->nrings);
 
        for ( i = 0; i < poly->nrings; i++ )
        {
                /* Set the coordinates (do write npoints) */
                ptarray_to_twkb_buf(poly->rings[i], globals, ts, 1, 4);
        }
 
        return 0;
}
 
 
 
/******************************************************************
* MULTI-GEOMETRYS (MultiPoint, MultiLinestring, MultiPolygon)
*******************************************************************/
 
static int lwmulti_to_twkb_buf(const LWCOLLECTION *col, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        uint32_t i;
        int nempty = 0;
 
        LWDEBUGF(2, "Entered %s", __func__);
        LWDEBUGF(4, "Number of geometries in multi is %d", col->ngeoms);
 
        /* Deal with special case for MULTIPOINT: skip any empty points */
        if ( col->type == MULTIPOINTTYPE )
        {
                for ( i = 0; i < col->ngeoms; i++ )
                        if ( lwgeom_is_empty(col->geoms[i]) )
                                nempty++;
        }
 
        /* Set the number of geometries */
        bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) (col->ngeoms - nempty));
 
        /* We've been handed an idlist, so write it in */
        if ( ts->idlist )
        {
                for ( i = 0; i < col->ngeoms; i++ )
                {
                        /* Skip empty points in multipoints, we can't represent them */
                        if ( col->type == MULTIPOINTTYPE && lwgeom_is_empty(col->geoms[i]) )
                                continue;
 
                        bytebuffer_append_varint(ts->geom_buf, ts->idlist[i]);
                }
 
                /* Empty it out to nobody else uses it now */
                ts->idlist = NULL;
        }
 
        for ( i = 0; i < col->ngeoms; i++ )
        {
                /* Skip empty points in multipoints, we can't represent them */
                if ( col->type == MULTIPOINTTYPE && lwgeom_is_empty(col->geoms[i]) )
                        continue;
 
                lwgeom_to_twkb_buf(col->geoms[i], globals, ts);
        }
        return 0;
}
 
/******************************************************************
* GEOMETRYCOLLECTIONS
*******************************************************************/
 
static int lwcollection_to_twkb_buf(const LWCOLLECTION *col, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        uint32_t i;
 
        LWDEBUGF(2, "Entered %s", __func__);
        LWDEBUGF(4, "Number of geometries in collection is %d", col->ngeoms);
 
        /* Set the number of geometries */
        bytebuffer_append_uvarint(ts->geom_buf, (uint64_t) col->ngeoms);
 
        /* We've been handed an idlist, so write it in */
        if ( ts->idlist )
        {
                for ( i = 0; i < col->ngeoms; i++ )
                        bytebuffer_append_varint(ts->geom_buf, ts->idlist[i]);
 
                /* Empty it out to nobody else uses it now */
                ts->idlist = NULL;
        }
 
        /* Write in the sub-geometries */
        for ( i = 0; i < col->ngeoms; i++ )
        {
                lwgeom_write_to_buffer(col->geoms[i], globals, ts);
        }
        return 0;
}
 
 
/******************************************************************
* Handle whole TWKB
*******************************************************************/
 
static int lwgeom_to_twkb_buf(const LWGEOM *geom, TWKB_GLOBALS *globals, TWKB_STATE *ts)
{
        LWDEBUGF(2, "Entered %s", __func__);
 
        switch ( geom->type )
        {
                case POINTTYPE:
                {
                        LWDEBUGF(4,"Type found is Point, %d", geom->type);
                        return lwpoint_to_twkb_buf((LWPOINT*) geom, globals, ts);
                }
                case LINETYPE:
                {
                        LWDEBUGF(4,"Type found is Linestring, %d", geom->type);
                        return lwline_to_twkb_buf((LWLINE*) geom, globals, ts);
                }
                case TRIANGLETYPE:
                {
                        LWDEBUGF(4, "Type found is Triangle, %d", geom->type);
                        return lwtriangle_to_twkb_buf((LWTRIANGLE *)geom, globals, ts);
                }
                /* Polygon has 'nrings' and 'rings' elements */
                case POLYGONTYPE:
                {
                        LWDEBUGF(4,"Type found is Polygon, %d", geom->type);
                        return lwpoly_to_twkb_buf((LWPOLY*)geom, globals, ts);
                }
 
                /* All these Collection types have 'ngeoms' and 'geoms' elements */
                case MULTIPOINTTYPE:
                case MULTILINETYPE:
                case MULTIPOLYGONTYPE:
                {
                        LWDEBUGF(4,"Type found is Multi, %d", geom->type);
                        return lwmulti_to_twkb_buf((LWCOLLECTION*)geom, globals, ts);
                }
                case COLLECTIONTYPE:
                case TINTYPE:
                {
                        LWDEBUGF(4,"Type found is collection, %d", geom->type);
                        return lwcollection_to_twkb_buf((LWCOLLECTION*) geom, globals, ts);
                }
                /* Unknown type! */
                default:
                        lwerror("%s: Unsupported geometry type: %s", __func__, lwtype_name(geom->type));
        }
 
        return 0;
}
 
 
static int lwgeom_write_to_buffer(const LWGEOM *geom, TWKB_GLOBALS *globals, TWKB_STATE *parent_state)
{
        int i, is_empty, has_z = 0, has_m = 0, ndims;
        size_t bbox_size = 0, optional_precision_byte = 0;
        uint8_t flag = 0, type_prec = 0;
        bytebuffer_t header_bytebuffer, geom_bytebuffer;
 
        TWKB_STATE child_state;
        memset(&child_state, 0, sizeof(TWKB_STATE));
        child_state.header_buf = &header_bytebuffer;
        child_state.geom_buf = &geom_bytebuffer;
        child_state.idlist = parent_state->idlist;
 
        bytebuffer_init_with_size(child_state.header_buf, 16);
        bytebuffer_init_with_size(child_state.geom_buf, 64);
 
        /* Read dimensionality from input */
        ndims = lwgeom_ndims(geom);
        is_empty = lwgeom_is_empty(geom);
        if ( ndims > 2 )
        {
                has_z = lwgeom_has_z(geom);
                has_m = lwgeom_has_m(geom);
        }
 
        /* Do we need extended precision? If we have a Z or M we do. */
        optional_precision_byte = (has_z || has_m);
 
        /* Both X and Y dimension use the same precision */
        globals->factor[0] = pow(10, globals->prec_xy);
        globals->factor[1] = globals->factor[0];
 
        /* Z and M dimensions have their own precisions */
        if ( has_z )
                globals->factor[2] = pow(10, globals->prec_z);
        if ( has_m )
                globals->factor[2 + has_z] = pow(10, globals->prec_m);
 
        /* Reset stats */
        for ( i = 0; i < MAX_N_DIMS; i++ )
        {
                /* Reset bbox calculation */
                child_state.bbox_max[i] = INT64_MIN;
                child_state.bbox_min[i] = INT64_MAX;
                /* Reset accumulated delta values to get absolute values on next point */
                child_state.accum_rels[i] = 0;
        }
 
        /* TYPE/PRECISION BYTE */
        if ( abs(globals->prec_xy) > 7 )
                lwerror("%s: X/Z precision cannot be greater than 7 or less than -7", __func__);
 
        /* Read the TWKB type number from the geometry */
        TYPE_PREC_SET_TYPE(type_prec, lwgeom_twkb_type(geom));
        /* Zig-zag the precision value before encoding it since it is a signed value */
        TYPE_PREC_SET_PREC(type_prec, zigzag8(globals->prec_xy));
        /* Write the type and precision byte */
        bytebuffer_append_byte(child_state.header_buf, type_prec);
 
        /* METADATA BYTE */
        /* Set first bit if we are going to store bboxes */
        FIRST_BYTE_SET_BBOXES(flag, (globals->variant & TWKB_BBOX) && ! is_empty);
        /* Set second bit if we are going to store resulting size */
        FIRST_BYTE_SET_SIZES(flag, globals->variant & TWKB_SIZE);
        /* There will be no ID-list (for now) */
        FIRST_BYTE_SET_IDLIST(flag, parent_state->idlist && ! is_empty);
        /* Are there higher dimensions */
        FIRST_BYTE_SET_EXTENDED(flag, optional_precision_byte);
        /* Empty? */
        FIRST_BYTE_SET_EMPTY(flag, is_empty);
        /* Write the header byte */
        bytebuffer_append_byte(child_state.header_buf, flag);
 
        /* EXTENDED PRECISION BYTE (OPTIONAL) */
        /* If needed, write the extended dim byte */
        if( optional_precision_byte )
        {
                uint8_t flag = 0;
 
                if ( has_z && ( globals->prec_z > 7 || globals->prec_z < 0 ) )
                        lwerror("%s: Z precision cannot be negative or greater than 7", __func__);
 
                if ( has_m && ( globals->prec_m > 7 || globals->prec_m < 0 ) )
                        lwerror("%s: M precision cannot be negative or greater than 7", __func__);
 
                HIGHER_DIM_SET_HASZ(flag, has_z);
                HIGHER_DIM_SET_HASM(flag, has_m);
                HIGHER_DIM_SET_PRECZ(flag, globals->prec_z);
                HIGHER_DIM_SET_PRECM(flag, globals->prec_m);
                bytebuffer_append_byte(child_state.header_buf, flag);
        }
 
        /* It the geometry is empty, we're almost done */
        if ( is_empty )
        {
                /* If this output is sized, write the size of */
                /* all following content, which is zero because */
                /* there is none */
                if ( globals->variant & TWKB_SIZE )
                        bytebuffer_append_byte(child_state.header_buf, 0);
 
                bytebuffer_append_bytebuffer(parent_state->geom_buf, child_state.header_buf);
                bytebuffer_destroy_buffer(child_state.header_buf);
                bytebuffer_destroy_buffer(child_state.geom_buf);
                return 0;
        }
 
        /* Write the TWKB into the output buffer */
        lwgeom_to_twkb_buf(geom, globals, &child_state);
 
        /*If we have a header_buf, we know that this function is called inside a collection*/
        /*and then we have to merge the bboxes of the included geometries*/
        /*and put the result to the parent (the collection)*/
        if( (globals->variant & TWKB_BBOX) && parent_state->header_buf )
        {
                LWDEBUG(4,"Merge bboxes");
                for ( i = 0; i < ndims; i++ )
                {
                        if(child_state.bbox_min[i]<parent_state->bbox_min[i])
                                parent_state->bbox_min[i] = child_state.bbox_min[i];
                        if(child_state.bbox_max[i]>parent_state->bbox_max[i])
                                parent_state->bbox_max[i] = child_state.bbox_max[i];
                }
        }
 
        /* Did we have a box? If so, how big? */
        bbox_size = 0;
        if( globals->variant & TWKB_BBOX )
        {
                LWDEBUG(4,"We want boxes and will calculate required size");
                bbox_size = sizeof_bbox(&child_state, ndims);
        }
 
        /* Write the size if wanted */
        if( globals->variant & TWKB_SIZE )
        {
                /* Here we have to add what we know will be written to header */
                /* buffer after size value is written */
                size_t size_to_register = bytebuffer_getlength(child_state.geom_buf);
                size_to_register += bbox_size;
                bytebuffer_append_uvarint(child_state.header_buf, size_to_register);
        }
 
        if( globals->variant & TWKB_BBOX )
                write_bbox(&child_state, ndims);
 
        bytebuffer_append_bytebuffer(parent_state->geom_buf,child_state.header_buf);
        bytebuffer_append_bytebuffer(parent_state->geom_buf,child_state.geom_buf);
 
        bytebuffer_destroy_buffer(child_state.header_buf);
        bytebuffer_destroy_buffer(child_state.geom_buf);
        return 0;
}
 
 
lwvarlena_t *
lwgeom_to_twkb_with_idlist(const LWGEOM *geom,
                           int64_t *idlist,
                           uint8_t variant,
                           int8_t precision_xy,
                           int8_t precision_z,
                           int8_t precision_m)
{
        LWDEBUGF(2, "Entered %s", __func__);
        LWDEBUGF(2, "variant value %x", variant);
 
        TWKB_GLOBALS tg;
        TWKB_STATE ts;
        bytebuffer_t geom_bytebuffer;
 
        memset(&ts, 0, sizeof(TWKB_STATE));
        memset(&tg, 0, sizeof(TWKB_GLOBALS));
 
        tg.variant = variant;
        tg.prec_xy = precision_xy;
        tg.prec_z = precision_z;
        tg.prec_m = precision_m;
 
        if ( idlist && ! lwgeom_is_collection(geom) )
        {
                lwerror("Only collections can support ID lists");
                return NULL;
        }
 
        if ( ! geom )
        {
                LWDEBUG(4,"Cannot convert NULL into TWKB.");
                lwerror("Cannot convert NULL into TWKB");
                return NULL;
        }
 
        ts.idlist = idlist;
        ts.header_buf = NULL;
        ts.geom_buf = &geom_bytebuffer;
        bytebuffer_init_with_size(ts.geom_buf, 512);
        lwgeom_write_to_buffer(geom, &tg, &ts);
 
        lwvarlena_t *v = bytebuffer_get_buffer_varlena(ts.geom_buf);
        bytebuffer_destroy_buffer(ts.geom_buf);
        return v;
}
 
lwvarlena_t *
lwgeom_to_twkb(const LWGEOM *geom, uint8_t variant, int8_t precision_xy, int8_t precision_z, int8_t precision_m)
{
        return lwgeom_to_twkb_with_idlist(geom, NULL, variant, precision_xy, precision_z, precision_m);
}