PostGIS  2.3.8dev-r@@SVN_REVISION@@

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

380 {
381  GEOSGeom gout;
382  GEOSGeom geos_bound;
383  GEOSGeom geos_cut_edges, geos_area, collapse_points;
384  GEOSGeometry *vgeoms[3]; /* One for area, one for cut-edges */
385  unsigned int nvgeoms=0;
386 
387  assert (GEOSGeomTypeId(gin) == GEOS_POLYGON ||
388  GEOSGeomTypeId(gin) == GEOS_MULTIPOLYGON);
389 
390  geos_bound = GEOSBoundary(gin);
391  if ( NULL == geos_bound )
392  {
393  return NULL;
394  }
395 
396  LWDEBUGF(3,
397  "Boundaries: %s",
398  lwgeom_to_ewkt(GEOS2LWGEOM(geos_bound, 0)));
399 
400  /* Use noded boundaries as initial "cut" edges */
401 
402 #ifdef LWGEOM_PROFILE_MAKEVALID
403  lwnotice("ST_MakeValid: noding lines");
404 #endif
405 
406 
407  geos_cut_edges = LWGEOM_GEOS_nodeLines(geos_bound);
408  if ( NULL == geos_cut_edges )
409  {
410  GEOSGeom_destroy(geos_bound);
411  lwnotice("LWGEOM_GEOS_nodeLines(): %s", lwgeom_geos_errmsg);
412  return NULL;
413  }
414 
415  /* NOTE: the noding process may drop lines collapsing to points.
416  * We want to retrive any of those */
417  {
418  GEOSGeometry* pi;
419  GEOSGeometry* po;
420 
421 #ifdef LWGEOM_PROFILE_MAKEVALID
422  lwnotice("ST_MakeValid: extracting unique points from bounds");
423 #endif
424 
425  pi = GEOSGeom_extractUniquePoints(geos_bound);
426  if ( NULL == pi )
427  {
428  GEOSGeom_destroy(geos_bound);
429  lwnotice("GEOSGeom_extractUniquePoints(): %s",
431  return NULL;
432  }
433 
434  LWDEBUGF(3,
435  "Boundaries input points %s",
436  lwgeom_to_ewkt(GEOS2LWGEOM(pi, 0)));
437 
438 #ifdef LWGEOM_PROFILE_MAKEVALID
439  lwnotice("ST_MakeValid: extracting unique points from cut_edges");
440 #endif
441 
442  po = GEOSGeom_extractUniquePoints(geos_cut_edges);
443  if ( NULL == po )
444  {
445  GEOSGeom_destroy(geos_bound);
446  GEOSGeom_destroy(pi);
447  lwnotice("GEOSGeom_extractUniquePoints(): %s",
449  return NULL;
450  }
451 
452  LWDEBUGF(3,
453  "Boundaries output points %s",
454  lwgeom_to_ewkt(GEOS2LWGEOM(po, 0)));
455 
456 #ifdef LWGEOM_PROFILE_MAKEVALID
457  lwnotice("ST_MakeValid: find collapse points");
458 #endif
459 
460  collapse_points = GEOSDifference(pi, po);
461  if ( NULL == collapse_points )
462  {
463  GEOSGeom_destroy(geos_bound);
464  GEOSGeom_destroy(pi);
465  GEOSGeom_destroy(po);
466  lwnotice("GEOSDifference(): %s", lwgeom_geos_errmsg);
467  return NULL;
468  }
469 
470  LWDEBUGF(3,
471  "Collapse points: %s",
472  lwgeom_to_ewkt(GEOS2LWGEOM(collapse_points, 0)));
473 
474 #ifdef LWGEOM_PROFILE_MAKEVALID
475  lwnotice("ST_MakeValid: cleanup(1)");
476 #endif
477 
478  GEOSGeom_destroy(pi);
479  GEOSGeom_destroy(po);
480  }
481  GEOSGeom_destroy(geos_bound);
482 
483  LWDEBUGF(3,
484  "Noded Boundaries: %s",
485  lwgeom_to_ewkt(GEOS2LWGEOM(geos_cut_edges, 0)));
486 
487  /* And use an empty geometry as initial "area" */
488  geos_area = GEOSGeom_createEmptyPolygon();
489  if ( ! geos_area )
490  {
491  lwnotice("GEOSGeom_createEmptyPolygon(): %s", lwgeom_geos_errmsg);
492  GEOSGeom_destroy(geos_cut_edges);
493  return NULL;
494  }
495 
496  /*
497  * See if an area can be build with the remaining edges
498  * and if it can, symdifference with the original area.
499  * Iterate this until no more polygons can be created
500  * with left-over edges.
501  */
502  while (GEOSGetNumGeometries(geos_cut_edges))
503  {
504  GEOSGeometry* new_area=0;
505  GEOSGeometry* new_area_bound=0;
506  GEOSGeometry* symdif=0;
507  GEOSGeometry* new_cut_edges=0;
508 
509 #ifdef LWGEOM_PROFILE_MAKEVALID
510  lwnotice("ST_MakeValid: building area from %d edges", GEOSGetNumGeometries(geos_cut_edges));
511 #endif
512 
513  /*
514  * ASSUMPTION: cut_edges should already be fully noded
515  */
516 
517  new_area = LWGEOM_GEOS_buildArea(geos_cut_edges);
518  if ( ! new_area ) /* must be an exception */
519  {
520  GEOSGeom_destroy(geos_cut_edges);
521  GEOSGeom_destroy(geos_area);
522  lwnotice("LWGEOM_GEOS_buildArea() threw an error: %s",
524  return NULL;
525  }
526 
527  if ( GEOSisEmpty(new_area) )
528  {
529  /* no more rings can be build with thes edges */
530  GEOSGeom_destroy(new_area);
531  break;
532  }
533 
534  /*
535  * We succeeded in building a ring !
536  */
537 
538 #ifdef LWGEOM_PROFILE_MAKEVALID
539  lwnotice("ST_MakeValid: ring built with %d cut edges, saving boundaries", GEOSGetNumGeometries(geos_cut_edges));
540 #endif
541 
542  /*
543  * Save the new ring boundaries first (to compute
544  * further cut edges later)
545  */
546  new_area_bound = GEOSBoundary(new_area);
547  if ( ! new_area_bound )
548  {
549  /* We did check for empty area already so
550  * this must be some other error */
551  lwnotice("GEOSBoundary('%s') threw an error: %s",
552  lwgeom_to_ewkt(GEOS2LWGEOM(new_area, 0)),
554  GEOSGeom_destroy(new_area);
555  GEOSGeom_destroy(geos_area);
556  return NULL;
557  }
558 
559 #ifdef LWGEOM_PROFILE_MAKEVALID
560  lwnotice("ST_MakeValid: running SymDifference with new area");
561 #endif
562 
563  /*
564  * Now symdif new and old area
565  */
566  symdif = GEOSSymDifference(geos_area, new_area);
567  if ( ! symdif ) /* must be an exception */
568  {
569  GEOSGeom_destroy(geos_cut_edges);
570  GEOSGeom_destroy(new_area);
571  GEOSGeom_destroy(new_area_bound);
572  GEOSGeom_destroy(geos_area);
573  lwnotice("GEOSSymDifference() threw an error: %s",
575  return NULL;
576  }
577 
578  GEOSGeom_destroy(geos_area);
579  GEOSGeom_destroy(new_area);
580  geos_area = symdif;
581  symdif = 0;
582 
583  /*
584  * Now let's re-set geos_cut_edges with what's left
585  * from the original boundary.
586  * ASSUMPTION: only the previous cut-edges can be
587  * left, so we don't need to reconsider
588  * the whole original boundaries
589  *
590  * NOTE: this is an expensive operation.
591  *
592  */
593 
594 #ifdef LWGEOM_PROFILE_MAKEVALID
595  lwnotice("ST_MakeValid: computing new cut_edges (GEOSDifference)");
596 #endif
597 
598  new_cut_edges = GEOSDifference(geos_cut_edges, new_area_bound);
599  GEOSGeom_destroy(new_area_bound);
600  if ( ! new_cut_edges ) /* an exception ? */
601  {
602  /* cleanup and throw */
603  GEOSGeom_destroy(geos_cut_edges);
604  GEOSGeom_destroy(geos_area);
605  /* TODO: Shouldn't this be an lwerror ? */
606  lwnotice("GEOSDifference() threw an error: %s",
608  return NULL;
609  }
610  GEOSGeom_destroy(geos_cut_edges);
611  geos_cut_edges = new_cut_edges;
612  }
613 
614 #ifdef LWGEOM_PROFILE_MAKEVALID
615  lwnotice("ST_MakeValid: final checks");
616 #endif
617 
618  if ( ! GEOSisEmpty(geos_area) )
619  {
620  vgeoms[nvgeoms++] = geos_area;
621  }
622  else
623  {
624  GEOSGeom_destroy(geos_area);
625  }
626 
627  if ( ! GEOSisEmpty(geos_cut_edges) )
628  {
629  vgeoms[nvgeoms++] = geos_cut_edges;
630  }
631  else
632  {
633  GEOSGeom_destroy(geos_cut_edges);
634  }
635 
636  if ( ! GEOSisEmpty(collapse_points) )
637  {
638  vgeoms[nvgeoms++] = collapse_points;
639  }
640  else
641  {
642  GEOSGeom_destroy(collapse_points);
643  }
644 
645  if ( 1 == nvgeoms )
646  {
647  /* Return cut edges */
648  gout = vgeoms[0];
649  }
650  else
651  {
652  /* Collect areas and lines (if any line) */
653  gout = GEOSGeom_createCollection(GEOS_GEOMETRYCOLLECTION, vgeoms, nvgeoms);
654  if ( ! gout ) /* an exception again */
655  {
656  /* cleanup and throw */
657  /* TODO: Shouldn't this be an lwerror ? */
658  lwnotice("GEOSGeom_createCollection() threw an error: %s",
660  /* TODO: cleanup! */
661  return NULL;
662  }
663  }
664 
665  return gout;
666 
667 }
void lwnotice(const char *fmt,...)
Write a notice out to the notice handler.
Definition: lwutil.c:89
char * lwgeom_to_ewkt(const LWGEOM *lwgeom)
Return an alloced string.
Definition: lwgeom.c:482
char lwgeom_geos_errmsg[LWGEOM_GEOS_ERRMSG_MAXSIZE]
LWGEOM * GEOS2LWGEOM(const GEOSGeometry *geom, char want3d)
#define LWDEBUGF(level, msg,...)
Definition: lwgeom_log.h:88
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: