geography_centroid_from_mpoly()

LWPOINT * geography_centroid_from_mpoly	(	const LWMPOLY *	mpoly,
		bool	use_spheroid,
		SPHEROID *	s
	)

Split polygons into triangles and use centroid of the triangle with the triangle area as weight to calculate the centroid of a (multi)polygon.

Definition at line 309 of file geography_centroid.c.

References geography_centroid_from_wpoints(), LWMPOLY::geoms, getPoint2d_cp(), LW_TRUE, lwgeom_area_sphere(), lwgeom_area_spheroid(), lwgeom_free(), lwgeom_set_geodetic(), lwpoint_free(), lwpoint_get_x(), lwpoint_get_y(), lwpoly_add_ring(), lwpoly_as_lwgeom(), lwpoly_construct_empty(), POINT3DM::m, LWMPOLY::ngeoms, POINTARRAY::npoints, LWPOLY::nrings, ptarray_construct_empty(), ptarray_insert_point(), LWPOLY::rings, LWMPOLY::srid, POINT3DM::x, POINT4D::x, POINT3DM::y, and POINT4D::y.

Referenced by geography_centroid().

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

Here is the call graph for this function:

Here is the caller graph for this function: