lwgeom_subdivide_recursive()

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

static

Definition at line 2265 of file lwgeom.c.

{
        const uint32_t maxdepth = 50;
        GBOX clip, subbox1, subbox2;
        uint32_t nvertices = 0;
        uint32_t i, n = 0;
        uint32_t split_ordinate;
        double width;
        double height;
        double pivot = DBL_MAX;
        double center = DBL_MAX;
        LWPOLY *lwpoly = NULL;
        LWGEOM *clipped;
        const GBOX *box_in;
        if (!geom) return 0;
        box_in = lwgeom_get_bbox(geom);
        if (!box_in) return 0;
        gbox_duplicate(box_in, &clip);
        width = clip.xmax - clip.xmin;
        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 1;
                }
                else
                        return 0;
        }
 
        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;
                int n = 0;
                /* Don't increment depth yet, since we aren't actually
                 * subdividing geometries yet */
                for ( i = 0; i < incol->ngeoms; i++ )
                        n += lwgeom_subdivide_recursive(incol->geoms[i], dimension, maxvertices, depth, col);
                return n;
        }
 
        if (lwgeom_dimension(geom) < dimension)
        {
                /* We've hit a lower dimension object produced by clipping at
                 * a shallower recursion level. Ignore it. */
                return 0;
        }
 
        /* 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 1;
        }
 
        nvertices = lwgeom_count_vertices(geom);
 
        /* Skip empties entirely */
        if (nvertices == 0) return 0;
 
        /* 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 1;
        }
 
        split_ordinate = (width > height) ? 0 : 1;
        if (split_ordinate == 0)
                center = (clip.xmin + clip.xmax) / 2;
        else
                center = (clip.ymin + clip.ymax) / 2;
 
        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;
                pivot = DBL_MAX;
                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 (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 (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_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);
        clipped = lwgeom_intersection(geom, subbox);
        lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
        lwgeom_free(subbox);
        if (clipped && !lwgeom_is_empty(clipped))
        {
                n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
                lwgeom_free(clipped);
        }
 
        subbox = (LWGEOM*) lwpoly_construct_envelope(geom->srid, subbox2.xmin, subbox2.ymin, subbox2.xmax, subbox2.ymax);
        clipped = lwgeom_intersection(geom, subbox);
        lwgeom_simplify_in_place(clipped, 0.0, LW_TRUE);
        lwgeom_free(subbox);
        if (clipped && !lwgeom_is_empty(clipped))
        {
                n += lwgeom_subdivide_recursive(clipped, dimension, maxvertices, depth, col);
                lwgeom_free(clipped);
        }
 
        return n;
}

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(), 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().

Here is the call graph for this function:

Here is the caller graph for this function: