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

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

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