lwgeom_wrapx.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 2016 Sandro Santilli <strk@kbt.io>
 *
 **********************************************************************/

#include "../postgis_config.h"
/*#define POSTGIS_DEBUG_LEVEL 4*/
#include "lwgeom_geos.h"
#include "liblwgeom_internal.h"

#include <string.h>
#include <assert.h>

LWGEOM* lwgeom_wrapx(const LWGEOM* lwgeom_in, double cutx, double amount);
static LWCOLLECTION* lwcollection_wrapx(const LWCOLLECTION* lwcoll_in, double cutx, double amount);

static LWGEOM*
lwgeom_split_wrapx(const LWGEOM* geom_in, double cutx, double amount)
{
        LWGEOM *blade, *split;
        POINTARRAY *bladepa;
        POINT4D pt;
        const GBOX *box_in;
        AFFINE affine = {
                1, 0, 0,
                0, 1, 0,
                0, 0, 1,
                amount, 0, 0,
        };

        /* Extract box */
        /* TODO: check if the bbox should be force-recomputed */
        box_in = lwgeom_get_bbox(geom_in);
        if ( ! box_in ) {
                /* must be empty */
                return lwgeom_clone_deep(geom_in);
        }

        LWDEBUGF(2, "BOX X range is %g..%g, cutx:%g, amount:%g", box_in->xmin, box_in->xmax, cutx, amount);

        /* Check if geometry is fully on the side needing shift */
        if ( ( amount < 0 && box_in->xmin >= cutx ) || ( amount > 0 && box_in->xmax <= cutx ) )
        {
                split = lwgeom_clone_deep(geom_in);
                lwgeom_affine(split, &affine);
                LWDEBUGG(2, split, "returning the translated geometry");
                return split;
        }

        /* Check if geometry is fully on the side needing no shift */
        if ( ( amount < 0 && box_in->xmax <= cutx ) || ( amount > 0 && box_in->xmin >= cutx ) )
        {
                split = lwgeom_clone_deep(geom_in);
                LWDEBUGG(2, split, "returning the cloned geometry");
                return split;
        }

        /* We need splitting here */

        /* construct blade */
        bladepa = ptarray_construct(0, 0, 2);
        pt.x = cutx;
        pt.y = box_in->ymin - 1;
        ptarray_set_point4d(bladepa, 0, &pt);
        pt.y = box_in->ymax + 1;
        ptarray_set_point4d(bladepa, 1, &pt);
        blade = lwline_as_lwgeom(lwline_construct(geom_in->srid, NULL, bladepa));

        LWDEBUG(2, "splitting the geometry");

        /* split by blade */
        split = lwgeom_split(geom_in, blade);
        lwgeom_free(blade);
        if ( ! split ) {
                lwerror("%s:%d - lwgeom_split_wrapx:  %s", __FILE__, __LINE__, lwgeom_geos_errmsg);
                return NULL;
        }
        LWDEBUGG(2, split, "split geometry");


        /* iterate over components, translate if needed */
        const LWCOLLECTION *col = lwgeom_as_lwcollection(split);
        if ( ! col ) {
                /* not split, this is unexpected */
                lwnotice("WARNING: unexpected lack of split in lwgeom_split_wrapx");
                return lwgeom_clone_deep(geom_in);
        }
        LWCOLLECTION *col_out = lwcollection_wrapx(col, cutx, amount);
        lwgeom_free(split);

        /* unary-union the result (homogenize too ?) */
        LWGEOM* out = lwgeom_unaryunion(lwcollection_as_lwgeom(col_out));
        LWDEBUGF(2, "col_out:%p, unaryunion_out:%p", col_out, out);
        LWDEBUGG(2, out, "unary-unioned");

        lwcollection_free(col_out);

        return out;
}

static LWCOLLECTION*
lwcollection_wrapx(const LWCOLLECTION* lwcoll_in, double cutx, double amount)
{
        LWGEOM** wrap_geoms;
        LWCOLLECTION* out;
        int i;
        int outtype = lwcoll_in->type;

        wrap_geoms = lwalloc(lwcoll_in->ngeoms * sizeof(LWGEOM*));
        if ( ! wrap_geoms )
        {
                lwerror("Out of virtual memory");
                return NULL;
        }

        for (i=0; i<lwcoll_in->ngeoms; ++i)
        {
                LWDEBUGF(3, "Wrapping collection element %d", i);
                wrap_geoms[i] = lwgeom_wrapx(lwcoll_in->geoms[i], cutx, amount);
                /* an exception should prevent this from ever returning NULL */
                if ( ! wrap_geoms[i] ) {
                        lwnotice("Error wrapping geometry, cleaning up");
                        while ((--i)>=0) {
                                lwnotice("cleaning geometry %d (%p)", i, wrap_geoms[i]);
                                lwgeom_free(wrap_geoms[i]);
                        }
                        lwfree(wrap_geoms);
                        lwnotice("cleanup complete");
                        return NULL;
                }
          if ( outtype != COLLECTIONTYPE ) {
                        if ( MULTITYPE[wrap_geoms[i]->type] != outtype )
                        {
                                outtype = COLLECTIONTYPE;
                        }
                }
        }

        /* Now wrap_geoms has wrap_geoms_size geometries */
        out = lwcollection_construct(outtype, lwcoll_in->srid, NULL,
                                     lwcoll_in->ngeoms, wrap_geoms);

        return out;
}

/* exported */
LWGEOM*
lwgeom_wrapx(const LWGEOM* lwgeom_in, double cutx, double amount)
{
        /* Nothing to wrap in an empty geom */
        if ( lwgeom_is_empty(lwgeom_in) )
        {
                LWDEBUG(2, "geom is empty, cloning");
                return lwgeom_clone_deep(lwgeom_in);
        }

        /* Nothing to wrap if shift amount is zero */
        if ( amount == 0 )
        {
                LWDEBUG(2, "amount is zero, cloning");
                return lwgeom_clone_deep(lwgeom_in);
        }

        switch (lwgeom_in->type)
        {
        case LINETYPE:
        case POLYGONTYPE:
                LWDEBUG(2, "split-wrapping line or polygon");
                return lwgeom_split_wrapx(lwgeom_in, cutx, amount);

        case POINTTYPE:
        {
                const LWPOINT *pt = lwgeom_as_lwpoint(lwgeom_clone_deep(lwgeom_in));
                POINT4D pt4d;
                getPoint4d_p(pt->point, 0, &pt4d);

                LWDEBUGF(2, "POINT X is %g, cutx:%g, amount:%g", pt4d.x, cutx, amount);

                if ( ( amount < 0 && pt4d.x > cutx ) || ( amount > 0 && pt4d.x < cutx ) )
                {
                        pt4d.x += amount;
                        ptarray_set_point4d(pt->point, 0, &pt4d);
                }
                return lwpoint_as_lwgeom(pt);
        }

        case MULTIPOINTTYPE:
        case MULTIPOLYGONTYPE:
        case MULTILINETYPE:
        case COLLECTIONTYPE:
                LWDEBUG(2, "collection-wrapping multi");
                return lwcollection_as_lwgeom(
                                                lwcollection_wrapx((const LWCOLLECTION*)lwgeom_in, cutx, amount)
                                         );

        default:
                lwerror("Wrapping of %s geometries is unsupported",
                        lwtype_name(lwgeom_in->type));
                return NULL;
        }

}