geography_centroid.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) 2017 Danny Götte <danny.goette@fem.tu-ilmenau.de>
 *
 **********************************************************************/

#include "postgres.h"

#include "../postgis_config.h"

#include <math.h>

#include "liblwgeom.h"         /* For standard geometry types. */
#include "lwgeom_pg.h"       /* For pg macros. */
#include "lwgeom_transform.h" /* For SRID functions */

Datum geography_centroid(PG_FUNCTION_ARGS);

/* internal functions */
LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size);
LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s);
LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s);
LWPOINT* cart_to_lwpoint(const double_t x_sum, const double_t y_sum, const double_t z_sum, const double_t weight_sum, const uint32_t srid);
POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat);

PG_FUNCTION_INFO_V1(geography_centroid);
Datum geography_centroid(PG_FUNCTION_ARGS)
{
        LWGEOM *lwgeom = NULL;
        LWGEOM *lwgeom_out = NULL;
    LWPOINT *lwpoint_out = NULL;
        GSERIALIZED *g = NULL;
        GSERIALIZED *g_out = NULL;
    uint32_t srid;
    bool use_spheroid = true;
    SPHEROID s;
        uint32_t type;

        /* Get our geometry object loaded into memory. */
        g = PG_GETARG_GSERIALIZED_P(0);
        lwgeom = lwgeom_from_gserialized(g);

        if (g == NULL)
        {
                PG_RETURN_NULL();
        }

        srid = lwgeom_get_srid(lwgeom);

        /* on empty input, return empty output */
        if (gserialized_is_empty(g))
        {
                LWCOLLECTION* empty = lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
                lwgeom_out = lwcollection_as_lwgeom(empty);
                lwgeom_set_geodetic(lwgeom_out, true);
                g_out = gserialized_from_lwgeom(lwgeom_out, 0);
                PG_RETURN_POINTER(g_out);
        }

    /* Initialize spheroid */
    spheroid_init_from_srid(fcinfo, srid, &s);

    /* Set to sphere if requested */
    use_spheroid = PG_GETARG_BOOL(1);
        if ( ! use_spheroid )
                s.a = s.b = s.radius;

        type = gserialized_get_type(g);

        switch (type)
        {

        case POINTTYPE:
    {
        /* centroid of a point is itself */
        PG_RETURN_POINTER(g);
        break;
    }

        case MULTIPOINTTYPE:
    {
        LWMPOINT* mpoints = lwgeom_as_lwmpoint(lwgeom);

        /* average between all points */
        uint32_t size = mpoints->ngeoms;
        POINT3DM* points = palloc(size*sizeof(POINT3DM));

                uint32_t i;
                for (i = 0; i < size; i++) {
            points[i].x = lwpoint_get_x(mpoints->geoms[i]);
            points[i].y = lwpoint_get_y(mpoints->geoms[i]);
            points[i].m = 1;
        }

                lwpoint_out = geography_centroid_from_wpoints(srid, points, size);
        pfree(points);
        break;
    }

        case LINETYPE:
    {
        LWLINE* line = lwgeom_as_lwline(lwgeom);

        /* reuse mline function */
        LWMLINE* mline = lwmline_construct_empty(srid, 0, 0);
        lwmline_add_lwline(mline, line);

        lwpoint_out = geography_centroid_from_mline(mline, &s);
        lwmline_free(mline);
                break;
    }

        case MULTILINETYPE:
    {
        LWMLINE* mline = lwgeom_as_lwmline(lwgeom);
        lwpoint_out = geography_centroid_from_mline(mline, &s);
                break;
    }

        case POLYGONTYPE:
    {
        LWPOLY* poly = lwgeom_as_lwpoly(lwgeom);

        /* reuse mpoly function */
        LWMPOLY* mpoly = lwmpoly_construct_empty(srid, 0, 0);
        lwmpoly_add_lwpoly(mpoly, poly);

        lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
        lwmpoly_free(mpoly);
        break;
    }

        case MULTIPOLYGONTYPE:
    {
        LWMPOLY* mpoly = lwgeom_as_lwmpoly(lwgeom);
        lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
        break;
    }
        default:
                elog(ERROR, "ST_Centroid(geography) unhandled geography type");
                PG_RETURN_NULL();
        }

        PG_FREE_IF_COPY(g, 0);

    lwgeom_out = lwpoint_as_lwgeom(lwpoint_out);
    lwgeom_set_geodetic(lwgeom_out, true);
    g_out = gserialized_from_lwgeom(lwgeom_out, 0);

        PG_RETURN_POINTER(g_out);
}


LWPOINT* geography_centroid_from_wpoints(const uint32_t srid, const POINT3DM* points, const uint32_t size)
{
    double_t x_sum = 0;
    double_t y_sum = 0;
    double_t z_sum = 0;
    double_t weight_sum = 0;

    double_t weight = 1;
    POINT3D* point;

        uint32_t i;
        for (i = 0; i < size; i++ )
    {
                point = lonlat_to_cart(points[i].x, points[i].y);
        weight = points[i].m;

        x_sum += point->x * weight;
        y_sum += point->y * weight;
        z_sum += point->z * weight;

        weight_sum += weight;

                lwfree(point);
    }

    return cart_to_lwpoint(x_sum, y_sum, z_sum, weight_sum, srid);
}

POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat)
{
    double_t lat, lon;
    double_t sin_lat;

    POINT3D* point = lwalloc(sizeof(POINT3D));;

    // prepare coordinate for trigonometric functions from [-90, 90] -> [0, pi]
    lat = (raw_lat + 90) / 180 * M_PI;

    // prepare coordinate for trigonometric functions from [-180, 180] -> [-pi, pi]
    lon = raw_lon / 180 * M_PI;

    /* calculate value only once */
    sin_lat = sinl(lat);

    /* convert to 3D cartesian coordinates */
    point->x = sin_lat * cosl(lon);
    point->y = sin_lat * sinl(lon);
    point->z = cosl(lat);

    return point;
}

LWPOINT* cart_to_lwpoint(const double_t x_sum, const double_t y_sum, const double_t z_sum, const double_t weight_sum, const uint32_t srid)
{
    double_t x = x_sum / weight_sum;
    double_t y = y_sum / weight_sum;
    double_t z = z_sum / weight_sum;

    /* x-y-z vector length */
    double_t r = sqrtl(powl(x, 2) + powl(y, 2) + powl(z, 2));

    double_t lon = atan2l(y, x) * 180 / M_PI;
    double_t lat = acosl(z / r) * 180 / M_PI - 90;

        return lwpoint_make2d(srid, lon, lat);
}

LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s)
{
    double_t tolerance = 0.0;
    uint32_t size = 0;
        uint32_t i, k;

    /* get total number of points */
    for (i = 0; i < mline->ngeoms; i++) {
        size += (mline->geoms[i]->points->npoints - 1) * 2;
    }

    POINT3DM* points = palloc(size*sizeof(POINT3DM));
    uint32_t j = 0;

    for (i = 0; i < mline->ngeoms; i++) {
        LWLINE* line = mline->geoms[i];

        /* add both points of line segment as weighted point */
        for (k = 0; k < line->points->npoints - 1; k++) {
            const POINT2D* p1 = getPoint2d_cp(line->points, k);
            const POINT2D* p2 = getPoint2d_cp(line->points, k+1);

            /* use line-segment length as weight */
            LWPOINT* lwp1 = lwpoint_make2d(mline->srid, p1->x, p1->y);
            LWPOINT* lwp2 = lwpoint_make2d(mline->srid, p2->x, p2->y);
            LWGEOM* lwgeom1 = lwpoint_as_lwgeom(lwp1);
            LWGEOM* lwgeom2 = lwpoint_as_lwgeom(lwp2);
            lwgeom_set_geodetic(lwgeom1, LW_TRUE);
            lwgeom_set_geodetic(lwgeom2, LW_TRUE);

            /* use point distance as weight */
            double_t weight = lwgeom_distance_spheroid(lwgeom1, lwgeom2, s, tolerance);

            points[j].x = p1->x;
            points[j].y = p1->y;
            points[j].m = weight;
            j++;

            points[j].x = p2->x;
            points[j].y = p2->y;
            points[j].m = weight;
            j++;

            lwgeom_free(lwgeom1);
            lwgeom_free(lwgeom2);
        }
    }

    LWPOINT* result = geography_centroid_from_wpoints(mline->srid, points, size);
    pfree(points);
    return result;
}


LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s)
{
    uint32_t size = 0;
        uint32_t i, ir, ip;
    for (ip = 0; ip < mpoly->ngeoms; ip++) {
                for (ir = 0; ir < mpoly->geoms[ip]->nrings; ir++) {
                size += mpoly->geoms[ip]->rings[ir]->npoints - 1;
                }
    }

    POINT3DM* points = palloc(size*sizeof(POINT3DM));
    uint32_t j = 0;

    /* use first point as reference to create triangles */
    const POINT4D* reference_point = (const POINT4D*) getPoint2d_cp(mpoly->geoms[0]->rings[0], 0);

    for (ip = 0; ip < mpoly->ngeoms; ip++) {
        LWPOLY* poly = mpoly->geoms[ip];

        for (ir = 0; ir < poly->nrings; ir++) {
            POINTARRAY* ring = poly->rings[ir];

            /* split into triangles (two points + reference point) */
            for (i = 0; i < ring->npoints - 1; i++) {
                LWPOLY* poly_tri;
                LWGEOM* geom_tri;
                LWPOINT* tri_centroid;
                POINT3DM triangle[3];
                double_t weight;
                const POINT4D* p1 = (const POINT4D*) getPoint2d_cp(ring, i);
                const POINT4D* p2 = (const POINT4D*) getPoint2d_cp(ring, i+1);
                POINTARRAY* pa = ptarray_construct_empty(0, 0, 4);

                ptarray_insert_point(pa, p1, 0);
                ptarray_insert_point(pa, p2, 1);
                ptarray_insert_point(pa, reference_point, 2);
                ptarray_insert_point(pa, p1, 3);

                poly_tri = lwpoly_construct_empty(mpoly->srid, 0, 0);
                lwpoly_add_ring(poly_tri, pa);

                geom_tri = lwpoly_as_lwgeom(poly_tri);
                lwgeom_set_geodetic(geom_tri, LW_TRUE);

                /* Calculate the weight of the triangle. If counter clockwise,
                 * the weight is negative (e.g. for holes in polygons)
                 */

                if ( use_spheroid )
                        weight = lwgeom_area_spheroid(geom_tri, s);
                else
                        weight = lwgeom_area_sphere(geom_tri, s);


                triangle[0].x = p1->x;
                triangle[0].y = p1->y;
                triangle[0].m = 1;

                triangle[1].x = p2->x;
                triangle[1].y = p2->y;
                triangle[1].m = 1;

                triangle[2].x = reference_point->x;
                triangle[2].y = reference_point->y;
                triangle[2].m = 1;

                /* get center of triangle */
                tri_centroid = geography_centroid_from_wpoints(mpoly->srid, triangle, 3);

                points[j].x = lwpoint_get_x(tri_centroid);
                points[j].y = lwpoint_get_y(tri_centroid);
                points[j].m = weight;
                j++;

                lwpoint_free(tri_centroid);
                lwgeom_free(geom_tri);
            }
        }
    }
    LWPOINT* result = geography_centroid_from_wpoints(mpoly->srid, points, size);
    pfree(points);
    return result;
}