PostGIS  2.5.1dev-r@@SVN_REVISION@@

◆ LWGEOM_GEOS_makeValidPolygon()

static GEOSGeometry* LWGEOM_GEOS_makeValidPolygon ( const GEOSGeometry *  gin)
static

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

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

Referenced by LWGEOM_GEOS_makeValid().

341 {
342  GEOSGeom gout;
343  GEOSGeom geos_bound;
344  GEOSGeom geos_cut_edges, geos_area, collapse_points;
345  GEOSGeometry* vgeoms[3]; /* One for area, one for cut-edges */
346  unsigned int nvgeoms = 0;
347 
348  assert(GEOSGeomTypeId(gin) == GEOS_POLYGON || GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);
349 
350  geos_bound = GEOSBoundary(gin);
351  if (!geos_bound) return NULL;
352 
353  /* Use noded boundaries as initial "cut" edges */
354 
355  geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
356  if (!geos_cut_edges)
357  {
358  GEOSGeom_destroy(geos_bound);
359  lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
360  return NULL;
361  }
362 
363  /* NOTE: the noding process may drop lines collapsing to points.
364  * We want to retrieve any of those */
365  {
366  GEOSGeometry* pi;
367  GEOSGeometry* po;
368 
369  pi = GEOSGeom_extractUniquePoints(geos_bound);
370  if (!pi)
371  {
372  GEOSGeom_destroy(geos_bound);
373  lwnotice("GEOSGeom_extractUniquePoints(): %s", lwgeom_geos_errmsg);
374  return NULL;
375  }
376 
377  po = GEOSGeom_extractUniquePoints(geos_cut_edges);
378  if (!po)
379  {
380  GEOSGeom_destroy(geos_bound);
381  GEOSGeom_destroy(pi);
382  lwnotice("GEOSGeom_extractUniquePoints(): %s", lwgeom_geos_errmsg);
383  return NULL;
384  }
385 
386  collapse_points = GEOSDifference(pi, po);
387  if (!collapse_points)
388  {
389  GEOSGeom_destroy(geos_bound);
390  GEOSGeom_destroy(pi);
391  GEOSGeom_destroy(po);
392  lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
393  return NULL;
394  }
395 
396  GEOSGeom_destroy(pi);
397  GEOSGeom_destroy(po);
398  }
399  GEOSGeom_destroy(geos_bound);
400 
401  /* And use an empty geometry as initial "area" */
402  geos_area = GEOSGeom_createEmptyPolygon();
403  if (!geos_area)
404  {
405  lwnotice("GEOSGeom_createEmptyPolygon(): %s", lwgeom_geos_errmsg);
406  GEOSGeom_destroy(geos_cut_edges);
407  return NULL;
408  }
409 
410  /*
411  * See if an area can be build with the remaining edges
412  * and if it can, symdifference with the original area.
413  * Iterate this until no more polygons can be created
414  * with left-over edges.
415  */
416  while (GEOSGetNumGeometries(geos_cut_edges))
417  {
418  GEOSGeometry* new_area = 0;
419  GEOSGeometry* new_area_bound = 0;
420  GEOSGeometry* symdif = 0;
421  GEOSGeometry* new_cut_edges = 0;
422 
423 #ifdef LWGEOM_PROFILE_MAKEVALID
424  lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges));
425 #endif
426 
427  /*
428  * ASSUMPTION: cut_edges should already be fully noded
429  */
430 
431  new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
432  if (!new_area) /* must be an exception */
433  {
434  GEOSGeom_destroy(geos_cut_edges);
435  GEOSGeom_destroy(geos_area);
436  lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s", lwgeom_geos_errmsg);
437  return NULL;
438  }
439 
440  if (GEOSisEmpty(new_area))
441  {
442  /* no more rings can be build with thes edges */
443  GEOSGeom_destroy(new_area);
444  break;
445  }
446 
447  /*
448  * We succeeded in building a ring!
449  * Save the new ring boundaries first (to compute
450  * further cut edges later)
451  */
452  new_area_bound = GEOSBoundary(new_area);
453  if (!new_area_bound)
454  {
455  /* We did check for empty area already so this must be some other error */
456  lwnotice("GEOSBoundary('%s') threw an error: %s",
457  lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0)),
459  GEOSGeom_destroy(new_area);
460  GEOSGeom_destroy(geos_area);
461  return NULL;
462  }
463 
464  /*
465  * Now symdif new and old area
466  */
467  symdif = GEOSSymDifference(geos_area, new_area);
468  if (!symdif) /* must be an exception */
469  {
470  GEOSGeom_destroy(geos_cut_edges);
471  GEOSGeom_destroy(new_area);
472  GEOSGeom_destroy(new_area_bound);
473  GEOSGeom_destroy(geos_area);
474  lwnotice("GEOSSymDifference() threw an error: %s", lwgeom_geos_errmsg);
475  return NULL;
476  }
477 
478  GEOSGeom_destroy(geos_area);
479  GEOSGeom_destroy(new_area);
480  geos_area = symdif;
481  symdif = 0;
482 
483  /*
484  * Now let's re-set geos_cut_edges with what's left
485  * from the original boundary.
486  * ASSUMPTION: only the previous cut-edges can be
487  * left, so we don't need to reconsider
488  * the whole original boundaries
489  *
490  * NOTE: this is an expensive operation.
491  *
492  */
493 
494  new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
495  GEOSGeom_destroy(new_area_bound);
496  if (!new_cut_edges) /* an exception ? */
497  {
498  /* cleanup and throw */
499  GEOSGeom_destroy(geos_cut_edges);
500  GEOSGeom_destroy(geos_area);
501  lwerror("GEOSDifference() threw an error: %s", lwgeom_geos_errmsg);
502  return NULL;
503  }
504  GEOSGeom_destroy(geos_cut_edges);
505  geos_cut_edges = new_cut_edges;
506  }
507 
508  if (!GEOSisEmpty(geos_area))
509  vgeoms[nvgeoms++] = geos_area;
510  else
511  GEOSGeom_destroy(geos_area);
512 
513  if (!GEOSisEmpty(geos_cut_edges))
514  vgeoms[nvgeoms++] = geos_cut_edges;
515  else
516  GEOSGeom_destroy(geos_cut_edges);
517 
518  if (!GEOSisEmpty(collapse_points))
519  vgeoms[nvgeoms++] = collapse_points;
520  else
521  GEOSGeom_destroy(collapse_points);
522 
523  if (1 == nvgeoms)
524  {
525  /* Return cut edges */
526  gout = vgeoms[0];
527  }
528  else
529  {
530  /* Collect areas and lines (if any line) */
531  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
532  if (!gout) /* an exception again */
533  {
534  /* cleanup and throw */
535  lwerror("GEOSGeom_createCollection() threw an error: %s", lwgeom_geos_errmsg);
536  /* TODO: cleanup! */
537  return NULL;
538  }
539  }
540 
541  return gout;
542 }
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, uint8_t want3d)
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:177
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:556
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
static GEOSGeometry * LWGEOM_GEOS_nodeLines(const GEOSGeometry *lines)
GEOSGeometry * LWGEOM_GEOS_buildArea(const GEOSGeometry *geom_in)
Here is the call graph for this function:
Here is the caller graph for this function: