PostGIS  2.2.7dev-r@@SVN_REVISION@@
static GEOSGeometry* LWGEOM_GEOS_makeValidPolygon ( const GEOSGeometry *  gin)
static

Definition at line 390 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().

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