Changeset 13553

Show
Ignore:
Timestamp:
01/20/08 07:07:33 (4 months ago)
Author:
rouault
Message:

Fix performance problems in OGRGeometryFactory::organizePolygons (bug #1217)

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved
  • branches/1.5/gdal/ogr/ogrgeometryfactory.cpp

    r13316 r13553  
    720720 * Analyse a set of rings (passed as simple polygons), and based on a  
    721721 * geometric analysis convert them into a polygon with inner rings,  
    722  * or a MultiPolygon if dealing with more than one polygon.  The  
    723  * contains analysis is done with GEOS if available, otherwise using a 
    724  * less robust approach based on ring envelopes.  
    725  * 
    726  * All the input geometries must be OGRPolygons with only an exterior 
    727  * ring and no interior rings.  
     722 * or a MultiPolygon if dealing with more than one polygon. 
     723 * 
     724 * All the input geometries must be OGRPolygons with only a valid exterior 
     725 * ring (at least 4 points) and no interior rings.  
    728726 * 
    729727 * The passed in geometries become the responsibility of the method, but the 
    730728 * papoPolygons "pointer array" remains owned by the caller. 
     729 * 
     730 * For faster computation, a polygon is considered to be inside 
     731 * another one if a single point of its external ring is included into the other one. 
     732 * (unless 'OGR_DEBUG_ORGANIZE_POLYGONS' configuration option is set to TRUE. 
     733 * In that case, a slower algorithm that tests exact topological relationships  
     734 * is used if GEOS is available.) 
    731735 *  
    732736 * @param papoPolygons array of geometry pointers - should all be OGRPolygons. 
     
    750754                                                   int *pbIsValidGeometry ) 
    751755{ 
     756    int bUseFastVersion; 
    752757    int i, j; 
    753758    OGRGeometry* geom = NULL; 
     
    767772    } 
    768773 
    769 /* -------------------------------------------------------------------- */ 
    770 /*      A wee bit of a warning.                                         */ 
    771 /* -------------------------------------------------------------------- */ 
    772     static int firstTime = 1; 
    773     if (!haveGEOS() && firstTime) 
    774     { 
    775         CPLDebug("OGR", 
    776                  "GDAL should be built with GEOS support enabled in order the " 
    777                  "SHAPE driver to provide reliable results on complex polygons."); 
    778         firstTime = 0; 
     774    if (CSLTestBoolean(CPLGetConfigOption("OGR_DEBUG_ORGANIZE_POLYGONS", "NO"))) 
     775    { 
     776        /* -------------------------------------------------------------------- */ 
     777        /*      A wee bit of a warning.                                         */ 
     778        /* -------------------------------------------------------------------- */ 
     779        static int firstTime = 1; 
     780        if (!haveGEOS() && firstTime) 
     781        { 
     782            CPLDebug("OGR", 
     783                    "In OGR_DEBUG_ORGANIZE_POLYGONS mode, GDAL should be built with GEOS support enabled in order " 
     784                    "OGRGeometryFactory::organizePolygons to provide reliable results on complex polygons."); 
     785            firstTime = 0; 
     786        } 
     787        bUseFastVersion = !haveGEOS(); 
     788    } 
     789    else 
     790    { 
     791        bUseFastVersion = TRUE; 
    779792    } 
    780793 
     
    794807        papoPolygons[i]->getEnvelope(&envelopes[i]); 
    795808 
     809        //fprintf(stderr, "[%d] %d points\n", i, ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints()); 
     810 
    796811        if( wkbFlatten(papoPolygons[i]->getGeometryType()) == wkbPolygon 
    797             && ((OGRPolygon *) papoPolygons[i])->getNumInteriorRings() == 0 ) 
     812            && ((OGRPolygon *) papoPolygons[i])->getNumInteriorRings() == 0 
     813            && ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints() >= 4) 
    798814        { 
    799815            areas[i] = ((OGRPolygon *)papoPolygons[i])->get_Area(); 
     
    806822                    CE_Warning, CPLE_AppDefined,  
    807823                    "organizePolygons() received an unexpected geometry.\n" 
    808                     "Either a polygon with interior rings, or a non-Polygon\n" 
    809                     "geometry.  Return arguments as a collection." ); 
     824                    "Either a polygon with interior rings, or a polygon with less than 4 points,\n" 
     825                    "or a non-Polygon geometry.  Return arguments as a collection." ); 
    810826                bMixedUpGeometries = TRUE; 
    811827            } 
     
    849865        for(j=i+1; go_on && j<nPolygonCount;j++) 
    850866        { 
    851             if (areas[i] > areas[j] 
    852                 && envelopes[i].Contains(envelopes[j]) 
    853                 && (!haveGEOS() ||papoPolygons[i]->Contains(papoPolygons[j]))) 
     867            OGRPoint pointI, pointJ; 
     868            const OGRLinearRing* exteriorRingI = ((OGRPolygon *)papoPolygons[i])->getExteriorRing(); 
     869            const OGRLinearRing* exteriorRingJ = ((OGRPolygon *)papoPolygons[j])->getExteriorRing(); 
     870            exteriorRingI->getPoint(0, &pointI); 
     871            exteriorRingJ->getPoint(0, &pointJ); 
     872 
     873            if (areas[i] > areas[j] && envelopes[i].Contains(envelopes[j]) && 
     874                ((bUseFastVersion && exteriorRingI->isPointInRing(&pointJ)) || 
     875                 (!bUseFastVersion && papoPolygons[i]->Contains(papoPolygons[j])))) 
    854876            { 
    855877                relations[i][j] = CONTAINS; 
    856878                relations[j][i] = IS_CONTAINED_BY; 
    857879            } 
    858             else if (areas[j] > areas[i] 
    859                      && envelopes[j].Contains(envelopes[i]) 
    860                      && (!haveGEOS()  
    861                          || papoPolygons[j]->Contains(papoPolygons[i])) ) 
     880            else if (areas[j] > areas[i] && envelopes[j].Contains(envelopes[i]) && 
     881                     ((bUseFastVersion && exteriorRingJ->isPointInRing(&pointI)) || 
     882                      (!bUseFastVersion && papoPolygons[j]->Contains(papoPolygons[i])))) 
    862883            { 
    863884                relations[j][i] = CONTAINS; 
    864885                relations[i][j] = IS_CONTAINED_BY; 
    865886            } 
    866  
    867887            /* We use Overlaps instead of Intersects to be more  
    868888               tolerant about touching polygons */  
    869             else if ( !envelopes[i].Intersects(envelopes[j]) 
    870                      || (haveGEOS() && !papoPolygons[i]->Overlaps(papoPolygons[j])) ) 
     889            else if ( bUseFastVersion || !envelopes[i].Intersects(envelopes[j]) 
     890                     || !papoPolygons[i]->Overlaps(papoPolygons[j]) ) 
    871891            { 
    872892                relations[i][j] = NOT_RELATED; 
     
    880900                   polygons */ 
    881901                go_on = FALSE; 
    882 #ifdef notdef 
     902#ifdef DEBUG 
     903                char* wkt1; 
     904                char* wkt2; 
     905                papoPolygons[i]->exportToWkt(&wkt1); 
     906                papoPolygons[j]->exportToWkt(&wkt2); 
    883907                CPLDebug( "OGR",  
    884908                          "Bad intersection for polygons %d and %d\n" 
     
    886910                          "geom %d: %s",  
    887911                          i, j, i, wkt1, j, wkt2 ); 
    888                 char* wkt1; 
    889                 char* wkt2; 
    890                 papoPolygons[i]->exportToWkt(&wkt1); 
    891                 papoPolygons[j]->exportToWkt(&wkt2); 
    892                 fprintf(stderr, "geom %d : %s\n", i, wkt1); 
    893                 fprintf(stderr, "geom %d : %s\n", j, wkt2); 
    894912                CPLFree(wkt1); 
    895913                CPLFree(wkt2); 
  • branches/1.5/gdal/ogr/ogrsf_frmts/bna/ogrbnalayer.cpp

    r13327 r13553  
    676676            if (!isValidGeometry) 
    677677            { 
    678 #ifdef HAVE_GEOS 
    679678                CPLError(CE_Warning, CPLE_AppDefined,  
    680679                        "Geometry of polygon of fid %d starting at line %d cannot be translated to Simple Geometry. " 
     
    682681                        fid, 
    683682                        offsetAndLineFeaturesTable[fid].line + 1); 
    684 #endif 
    685683            } 
    686684        } 
  • branches/1.5/gdal/ogr/ogrsf_frmts/shape/shape2ogr.cpp

    r13317 r13553  
    460460            poOGRPoly->addRingDirectly( poRing ); 
    461461        } 
    462         else if( OGRGeometryFactory::haveGEOS() ) 
     462 
     463        // OGRGeometryFactory::organizePolygons does not require GOES anymore  
     464        else if( TRUE /* OGRGeometryFactory::haveGEOS()*/ ) 
    463465        { 
    464466            OGRPolygon** tabPolygons = new OGRPolygon*[psShape->nParts];