lwgeom_subdivide_recursive()

static void lwgeom_subdivide_recursive ( const LWGEOM *  geom, uint8_t  dimension, uint32_t  maxvertices, uint32_t  depth, LWCOLLECTION *  col  )

static

Definition at line 2261 of file lwgeom.c.

{
        const uint32_t maxdepth = 50;
 
        if (!geom)
                return;
 
        const GBOX *box_in = lwgeom_get_bbox(geom);
        if (!box_in)
                return;
 
        GBOX clip;
        gbox_duplicate(box_in, &clip);
        double width = clip.xmax - clip.xmin;
        double height = clip.ymax - clip.ymin;
 
        if ( geom->type == POLYHEDRALSURFACETYPE || geom->type == TINTYPE )
                lwerror("%s: unsupported geometry type '%s'", __func__, lwtype_name(geom->type));
 
        if ( width == 0.0 && height == 0.0 )
        {
                if ( geom->type == POINTTYPE && dimension == 0)
                        lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
                return;
        }
 
        if (width == 0.0)
        {
                clip.xmax += FP_TOLERANCE;
                clip.xmin -= FP_TOLERANCE;
                width = 2 * FP_TOLERANCE;
        }
        if (height == 0.0)
        {
                clip.ymax += FP_TOLERANCE;
                clip.ymin -= FP_TOLERANCE;
                height = 2 * FP_TOLERANCE;
        }
 
        /* Always just recurse into collections */
        if ( lwgeom_is_collection(geom) && geom->type != MULTIPOINTTYPE )
        {
                LWCOLLECTION *incol = (LWCOLLECTION*)geom;
                /* Don't increment depth yet, since we aren't actually
                 * subdividing geometries yet */
                for (uint32_t i = 0; i < incol->ngeoms; i++ )
                        lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col);
                return;
        }
 
        if (lwgeom_dimension(geom) < dimension)
        {
                /* We've hit a lower dimension object produced by clipping at
                 * a shallower recursion level. Ignore it. */
                return;
        }
 
        /* But don't go too far. 2^50 ~= 10^15, that's enough subdivision */
        /* Just add what's left */
        if ( depth > maxdepth )
        {
                lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
                return;
        }
 
        uint32_t nvertices = lwgeom_count_vertices(geom);
 
        /* Skip empties entirely */
        if (nvertices == 0)
                return;
 
        /* If it is under the vertex tolerance, just add it, we're done */
        if (nvertices <= maxvertices)
        {
                lwcollection_add_lwgeom(col, lwgeom_clone_deep(geom));
                return;
        }
 
        uint8_t split_ordinate = (width > height) ? 0 : 1;
        double center = (split_ordinate == 0) ? (clip.xmin + clip.xmax) / 2 : (clip.ymin + clip.ymax) / 2;
        double pivot = DBL_MAX;
        if (geom->type == POLYGONTYPE)
        {
                uint32_t ring_to_trim = 0;
                double ring_area = 0;
                double pivot_eps = DBL_MAX;
                double pt_eps = DBL_MAX;
                POINTARRAY *pa;
                LWPOLY *lwpoly = (LWPOLY *)geom;
 
                /* if there are more points in holes than in outer ring */
                if (nvertices >= 2 * lwpoly->rings[0]->npoints)
                {
                        /* trim holes starting from biggest */
                        for (uint32_t i = 1; i < lwpoly->nrings; i++)
                        {
                                double current_ring_area = fabs(ptarray_signed_area(lwpoly->rings[i]));
                                if (current_ring_area >= ring_area)
                                {
                                        ring_area = current_ring_area;
                                        ring_to_trim = i;
                                }
                        }
                }
 
                pa = lwpoly->rings[ring_to_trim];
 
                /* find most central point in chosen ring */
                for (uint32_t i = 0; i < pa->npoints; i++)
                {
                        double pt;
                        if (split_ordinate == 0)
                                pt = getPoint2d_cp(pa, i)->x;
                        else
                                pt = getPoint2d_cp(pa, i)->y;
                        pt_eps = fabs(pt - center);
                        if (pivot_eps > pt_eps)
                        {
                                pivot = pt;
                                pivot_eps = pt_eps;
                        }
                }
        }
        GBOX subbox1, subbox2;
        gbox_duplicate(&clip, &subbox1);
        gbox_duplicate(&clip, &subbox2);
 
        if (pivot == DBL_MAX)
                pivot = center;
 
        if (split_ordinate == 0)
        {
                if (FP_NEQUALS(subbox1.xmax, pivot) && FP_NEQUALS(subbox1.xmin, pivot))
                        subbox1.xmax = subbox2.xmin = pivot;
                else
                        subbox1.xmax = subbox2.xmin = center;
        }
        else
        {
                if (FP_NEQUALS(subbox1.ymax, pivot) && FP_NEQUALS(subbox1.ymin, pivot))
                        subbox1.ymax = subbox2.ymin = pivot;
                else
                        subbox1.ymax = subbox2.ymin = center;
        }
 
        ++depth;
 
        {
                LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
                    geom->srid, subbox1.xmin, subbox1.ymin, subbox1.xmax, subbox1.ymax);
                LWGEOM *clipped = lwgeom_intersection(geom, subbox);
                lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
                lwgeom_free(subbox);
                if (clipped && !lwgeom_is_empty(clipped))
                {
                        lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
                        lwgeom_free(clipped);
                }
        }
        {
                LWGEOM *subbox = (LWGEOM *)lwpoly_construct_envelope(
                    geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
                LWGEOM *clipped = lwgeom_intersection(geom, subbox);
                lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
                lwgeom_free(subbox);
                if (clipped && !lwgeom_is_empty(clipped))
                {
                        lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
                        lwgeom_free(clipped);
                }
        }
}

References FP_NEQUALS, FP_TOLERANCE, gbox_duplicate(), LWCOLLECTION::geoms, getPoint2d_cp(), LW_TRUE, lwcollection_add_lwgeom(), lwerror(), lwgeom_clone_deep(), lwgeom_count_vertices(), lwgeom_dimension(), lwgeom_free(), lwgeom_get_bbox(), lwgeom_intersection(), lwgeom_is_collection(), lwgeom_is_empty(), lwgeom_simplify_in_place(), lwgeom_subdivide_recursive(), lwpoly_construct_envelope(), lwtype_name(), MULTIPOINTTYPE, LWCOLLECTION::ngeoms, POINTARRAY::npoints, LWPOLY::nrings, pivot(), POINTTYPE, POLYGONTYPE, POLYHEDRALSURFACETYPE, ptarray_signed_area(), LWPOLY::rings, LWGEOM::srid, TINTYPE, LWGEOM::type, POINT2D::x, GBOX::xmax, GBOX::xmin, POINT2D::y, GBOX::ymax, and GBOX::ymin.

Referenced by lwgeom_subdivide(), and lwgeom_subdivide_recursive().

Here is the call graph for this function:

Here is the caller graph for this function: