LWGEOM_GEOS_makeValidPolygon()

static GEOSGeometry* LWGEOM_GEOS_makeValidPolygon ( const GEOSGeometry *  gin )

static

Definition at line 379 of file liblwgeom/lwgeom_geos_clean.c.

References GEOS2LWGEOM(), LWDEBUGF, LWGEOM_GEOS_buildArea(), lwgeom_geos_errmsg, LWGEOM_GEOS_nodeLines(), lwgeom_to_ewkt(), and lwnotice().

Referenced by LWGEOM_GEOS_makeValid().

{
        GEOSGeom gout;
        GEOSGeom geos_bound;
        GEOSGeom geos_cut_edges, geos_area, collapse_points;
        GEOSGeometry *vgeoms[3]; /* One for area, one for cut-edges */
        unsigned int nvgeoms=0;

        assert (GEOSGeomTypeId(gin) == GEOS_POLYGON ||
                GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);

        geos_bound = GEOSBoundary(gin);
        if ( NULL == geos_bound )
        {
                return NULL;
        }

        LWDEBUGF(3,
                       "Boundaries: %s",
                       lwgeom_to_ewkt(GEOS2LWGEOM(geos_bound, 0)));

        /* Use noded boundaries as initial "cut" edges */

#ifdef LWGEOM_PROFILE_MAKEVALID
  lwnotice("ST_MakeValid: noding lines");
#endif


        geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
        if ( NULL == geos_cut_edges )
        {
                GEOSGeom_destroy(geos_bound);
                lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
                return NULL;
        }

        /* NOTE: the noding process may drop lines collapsing to points.
         *       We want to retrive any of those */
        {
                GEOSGeometry* pi;
                GEOSGeometry* po;

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: extracting unique points from bounds");
#endif

                pi = GEOSGeom_extractUniquePoints(geos_bound);
                if ( NULL == pi )
                {
                        GEOSGeom_destroy(geos_bound);
                        lwnotice("GEOSGeom_extractUniquePoints(): %s",
                                 lwgeom_geos_errmsg);
                        return NULL;
                }

                LWDEBUGF(3,
                               "Boundaries input points %s",
                               lwgeom_to_ewkt(GEOS2LWGEOM(pi, 0)));

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: extracting unique points from cut_edges");
#endif

                po = GEOSGeom_extractUniquePoints(geos_cut_edges);
                if ( NULL == po )
                {
                        GEOSGeom_destroy(geos_bound);
                        GEOSGeom_destroy(pi);
                        lwnotice("GEOSGeom_extractUniquePoints(): %s",
                                 lwgeom_geos_errmsg);
                        return NULL;
                }

                LWDEBUGF(3,
                               "Boundaries output points %s",
                               lwgeom_to_ewkt(GEOS2LWGEOM(po, 0)));

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: find collapse points");
#endif

                collapse_points = GEOSDifference(pi, po);
                if ( NULL == collapse_points )
                {
                        GEOSGeom_destroy(geos_bound);
                        GEOSGeom_destroy(pi);
                        GEOSGeom_destroy(po);
                        lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
                        return NULL;
                }

                LWDEBUGF(3,
                               "Collapse points: %s",
                               lwgeom_to_ewkt(GEOS2LWGEOM(collapse_points, 0)));

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: cleanup(1)");
#endif

                GEOSGeom_destroy(pi);
                GEOSGeom_destroy(po);
        }
        GEOSGeom_destroy(geos_bound);

        LWDEBUGF(3,
                       "Noded Boundaries: %s",
                       lwgeom_to_ewkt(GEOS2LWGEOM(geos_cut_edges, 0)));

        /* And use an empty geometry as initial "area" */
        geos_area = GEOSGeom_createEmptyPolygon();
        if ( ! geos_area )
        {
                lwnotice("GEOSGeom_createEmptyPolygon(): %s", lwgeom_geos_errmsg);
                GEOSGeom_destroy(geos_cut_edges);
                return NULL;
        }

        /*
         * See if an area can be build with the remaining edges
         * and if it can, symdifference with the original area.
         * Iterate this until no more polygons can be created
         * with left-over edges.
         */
        while (GEOSGetNumGeometries(geos_cut_edges))
        {
                GEOSGeometry* new_area=0;
                GEOSGeometry* new_area_bound=0;
                GEOSGeometry* symdif=0;
                GEOSGeometry* new_cut_edges=0;

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges));
#endif

                /*
                 * ASSUMPTION: cut_edges should already be fully noded
                 */

                new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
                if ( ! new_area )   /* must be an exception */
                {
                        GEOSGeom_destroy(geos_cut_edges);
                        GEOSGeom_destroy(geos_area);
                        lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s",
                                 lwgeom_geos_errmsg);
                        return NULL;
                }

                if ( GEOSisEmpty(new_area) )
                {
                        /* no more rings can be build with thes edges */
                        GEOSGeom_destroy(new_area);
                        break;
                }

                /*
                 * We succeeded in building a ring !
                 */

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: ring built with %d cut edges, saving boundaries", GEOSGetNumGeometries(geos_cut_edges));
#endif

                /*
                 * Save the new ring boundaries first (to compute
                 * further cut edges later)
                 */
                new_area_bound = GEOSBoundary(new_area);
                if ( ! new_area_bound )
                {
                        /* We did check for empty area already so
                         * this must be some other error */
                        lwnotice("GEOSBoundary('%s') threw an error: %s",
                                 lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0)),
                                 lwgeom_geos_errmsg);
                        GEOSGeom_destroy(new_area);
                        GEOSGeom_destroy(geos_area);
                        return NULL;
                }

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: running SymDifference with new area");
#endif

                /*
                 * Now symdif new and old area
                 */
                symdif = GEOSSymDifference(geos_area, new_area);
                if ( ! symdif )   /* must be an exception */
                {
                        GEOSGeom_destroy(geos_cut_edges);
                        GEOSGeom_destroy(new_area);
                        GEOSGeom_destroy(new_area_bound);
                        GEOSGeom_destroy(geos_area);
                        lwnotice("GEOSSymDifference() threw an error: %s",
                                 lwgeom_geos_errmsg);
                        return NULL;
                }

                GEOSGeom_destroy(geos_area);
                GEOSGeom_destroy(new_area);
                geos_area = symdif;
                symdif = 0;

                /*
                 * Now let's re-set geos_cut_edges with what's left
                 * from the original boundary.
                 * ASSUMPTION: only the previous cut-edges can be
                 *             left, so we don't need to reconsider
                 *             the whole original boundaries
                 *
                 * NOTE: this is an expensive operation.
                 *
                 */

#ifdef LWGEOM_PROFILE_MAKEVALID
    lwnotice("ST_MakeValid: computing new cut_edges (GEOSDifference)");
#endif

                new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
                GEOSGeom_destroy(new_area_bound);
                if ( ! new_cut_edges )   /* an exception ? */
                {
                        /* cleanup and throw */
                        GEOSGeom_destroy(geos_cut_edges);
                        GEOSGeom_destroy(geos_area);
                        /* TODO: Shouldn't this be an lwerror ? */
                        lwnotice("GEOSDifference() threw an error: %s",
                                 lwgeom_geos_errmsg);
                        return NULL;
                }
                GEOSGeom_destroy(geos_cut_edges);
                geos_cut_edges = new_cut_edges;
        }

#ifdef LWGEOM_PROFILE_MAKEVALID
  lwnotice("ST_MakeValid: final checks");
#endif

        if ( ! GEOSisEmpty(geos_area) )
        {
                vgeoms[nvgeoms++] = geos_area;
        }
        else
        {
                GEOSGeom_destroy(geos_area);
        }

        if ( ! GEOSisEmpty(geos_cut_edges) )
        {
                vgeoms[nvgeoms++] = geos_cut_edges;
        }
        else
        {
                GEOSGeom_destroy(geos_cut_edges);
        }

        if ( ! GEOSisEmpty(collapse_points) )
        {
                vgeoms[nvgeoms++] = collapse_points;
        }
        else
        {
                GEOSGeom_destroy(collapse_points);
        }

        if ( 1 == nvgeoms )
        {
                /* Return cut edges */
                gout = vgeoms[0];
        }
        else
        {
                /* Collect areas and lines (if any line) */
                gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
                if ( ! gout )   /* an exception again */
                {
                        /* cleanup and throw */
                        /* TODO: Shouldn't this be an lwerror ? */
                        lwnotice("GEOSGeom_createCollection() threw an error: %s",
                                 lwgeom_geos_errmsg);
                        /* TODO: cleanup! */
                        return NULL;
                }
        }

        return gout;

}

Here is the call graph for this function:

Here is the caller graph for this function: