Changeset 14734

Show
Ignore:
Timestamp:
06/20/08 15:40:35 (5 months ago)
Author:
rouault
Message:

Improve memory efficiency (and possibly speed) of OGRGeometryFactory::organizePolygons (#2428)

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • branches/1.5/gdal/ogr/ogrgeometryfactory.cpp

    r13553 r14734  
    743743 */ 
    744744 
    745 enum 
    746 { 
    747     CONTAINS, 
    748     IS_CONTAINED_BY, 
    749     NOT_RELATED 
     745typedef struct _sPolyExtended sPolyExtended; 
     746 
     747struct _sPolyExtended 
     748{ 
     749    OGRPolygon*     poPolygon; 
     750    OGREnvelope     sEnvelope; 
     751    OGRLinearRing*  poExteriorRing; 
     752    OGRPoint        poAPoint; 
     753    int             nInitialIndex; 
     754    int             bIsTopLevel; 
     755    OGRPolygon*     poEnclosingPolygon; 
     756    double          dfArea; 
    750757}; 
     758 
     759static int OGRGeometryFactoryCompareArea(const void* p1, const void* p2) 
     760{ 
     761    const sPolyExtended* psPoly1 = (const sPolyExtended*) p1; 
     762    const sPolyExtended* psPoly2 = (const sPolyExtended*) p2; 
     763    if (psPoly2->dfArea < psPoly1->dfArea) 
     764        return -1; 
     765    else if (psPoly2->dfArea > psPoly1->dfArea) 
     766        return 1; 
     767    else 
     768        return 0; 
     769} 
     770 
     771static int OGRGeometryFactoryCompareByIndex(const void* p1, const void* p2) 
     772{ 
     773    const sPolyExtended* psPoly1 = (const sPolyExtended*) p1; 
     774    const sPolyExtended* psPoly2 = (const sPolyExtended*) p2; 
     775    if (psPoly1->nInitialIndex < psPoly2->nInitialIndex) 
     776        return -1; 
     777    else if (psPoly1->nInitialIndex > psPoly2->nInitialIndex) 
     778        return 1; 
     779    else 
     780        return 0; 
     781} 
     782 
    751783 
    752784OGRGeometry* OGRGeometryFactory::organizePolygons( OGRGeometry **papoPolygons, 
    …… 
    793825 
    794826/* -------------------------------------------------------------------- */ 
    795 /*      Setup per polygon relation, envelope and area information.      */ 
    796 /* -------------------------------------------------------------------- */ 
    797     OGREnvelope* envelopes = new OGREnvelope[nPolygonCount]; 
    798     int** relations = new int*[nPolygonCount]; 
    799     double* areas = new double[nPolygonCount]; 
     827/*      Setup per polygon envelope and area information.                */ 
     828/* -------------------------------------------------------------------- */ 
     829    sPolyExtended* asPolyEx = new sPolyExtended[nPolygonCount]; 
     830 
    800831    int go_on = TRUE; 
    801832    int bMixedUpGeometries = FALSE; 
    …… 
    804835    for(i=0;i<nPolygonCount;i++) 
    805836    { 
    806         relations[i] = new int[nPolygonCount]; 
    807         papoPolygons[i]->getEnvelope(&envelopes[i]); 
    808  
    809         //fprintf(stderr, "[%d] %d points\n", i, ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints()); 
     837        asPolyEx[i].nInitialIndex = i; 
     838        asPolyEx[i].poPolygon = (OGRPolygon*)papoPolygons[i]; 
     839        papoPolygons[i]->getEnvelope(&asPolyEx[i].sEnvelope); 
    810840 
    811841        if( wkbFlatten(papoPolygons[i]->getGeometryType()) == wkbPolygon 
    …… 
    813843            && ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints() >= 4) 
    814844        { 
    815             areas[i] = ((OGRPolygon *)papoPolygons[i])->get_Area(); 
     845            asPolyEx[i].dfArea = asPolyEx[i].poPolygon->get_Area(); 
     846            asPolyEx[i].poExteriorRing = asPolyEx[i].poPolygon->getExteriorRing(); 
     847            asPolyEx[i].poExteriorRing->getPoint(0, &asPolyEx[i].poAPoint); 
    816848        } 
    817849        else 
    …… 
    829861                bNonPolygon = TRUE; 
    830862        } 
    831              
    832     } 
    833          
    834     /* This a several step algorithm : 
    835        1) Compute in a matrix how polygons relate to each other 
    836        (this is the moment for detecting pathological intersections and exiting) 
    837        2) For each polygon, find the smallest enclosing polygon 
    838        3) For each polygon, compute its inclusion depth (0 means toplevel) 
    839        4) For each polygon of odd depth (= inner ring), add it to its outer ring 
    840          
     863    } 
     864 
     865 
     866    /* This a several steps algorithm : 
     867       1) Sort polygons by descending areas 
     868       2) For each polygon of rank i, find its smallest enclosing polygon 
     869          among the polygons of rank [i-1 ... 0]. If there are no such polygon, 
     870          this is a toplevel polygon. Otherwise, depending on if the enclosing 
     871          polygon is toplevel or not, we can decide if we are toplevel or not 
     872       3) Re-sort the polygons to retrieve their inital order (nicer for some applications) 
     873       4) For each non toplevel polygon (= inner ring), add it to its outer ring 
     874       5) Add the toplevel polygons to the multipolygon 
     875 
    841876       Complexity : O(nPolygonCount^2) 
    842877    */ 
    843          
     878 
    844879    /* Compute how each polygon relate to the other ones 
    845880       To save a bit of computation we always begin the computation by a test  
    …… 
    858893    */ 
    859894 
     895    if (!bMixedUpGeometries) 
     896    { 
     897        /* STEP 1 : Sort polygons by descending area */ 
     898        qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareArea); 
     899    } 
     900    papoPolygons = NULL; /* just to use to avoid it afterwards */ 
     901 
    860902/* -------------------------------------------------------------------- */ 
    861903/*      Compute relationships, if things seem well structured.          */ 
    862904/* -------------------------------------------------------------------- */ 
    863     for(i=0; !bMixedUpGeometries && go_on && i<nPolygonCount; i++) 
    864     { 
    865         for(j=i+1; go_on && j<nPolygonCount;j++) 
    866         { 
    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])))) 
     905 
     906    /* The first (largest) polygon is necessarily top-level */ 
     907    asPolyEx[0].bIsTopLevel = TRUE; 
     908    asPolyEx[0].poEnclosingPolygon = NULL; 
     909 
     910    int nCountTopLevel = 1; 
     911 
     912    /* STEP 2 */ 
     913    for(i=1; !bMixedUpGeometries && go_on && i<nPolygonCount; i++) 
     914    { 
     915        for(j=i-1; go_on && j>=0;j--) 
     916        { 
     917            if (asPolyEx[j].sEnvelope.Contains(asPolyEx[i].sEnvelope) && 
     918                ((bUseFastVersion && asPolyEx[j].poExteriorRing->isPointInRing(&asPolyEx[i].poAPoint)) || 
     919                 (!bUseFastVersion && asPolyEx[j].poPolygon->Contains(asPolyEx[i].poPolygon)))) 
    876920            { 
    877                 relations[i][j] = CONTAINS; 
    878                 relations[j][i] = IS_CONTAINED_BY; 
    879             } 
    880             else if (areas[j] > areas[i] && envelopes[j].Contains(envelopes[i]) && 
    881                      ((bUseFastVersion && exteriorRingJ->isPointInRing(&pointI)) || 
    882                       (!bUseFastVersion && papoPolygons[j]->Contains(papoPolygons[i])))) 
    883             { 
    884                 relations[j][i] = CONTAINS; 
    885                 relations[i][j] = IS_CONTAINED_BY; 
     921                if (asPolyEx[j].bIsTopLevel) 
     922                { 
     923                    /* We are a lake */ 
     924                    asPolyEx[i].bIsTopLevel = FALSE; 
     925                    asPolyEx[i].poEnclosingPolygon = asPolyEx[j].poPolygon; 
     926                } 
     927                else 
     928                { 
     929                    /* We are included in a something not toplevel (a lake), */ 
     930                    /* so in OGCSF we are considered as toplevel too */ 
     931                    nCountTopLevel ++; 
     932                    asPolyEx[i].bIsTopLevel = TRUE; 
     933                    asPolyEx[i].poEnclosingPolygon = NULL; 
     934                } 
     935                break; 
    886936            } 
    887937            /* We use Overlaps instead of Intersects to be more  
    888938               tolerant about touching polygons */  
    889             else if ( bUseFastVersion || !envelopes[i].Intersects(envelopes[j]) 
    890                      || !papoPolygons[i]->Overlaps(papoPolygons[j]) ) 
     939            else if ( bUseFastVersion || !asPolyEx[i].sEnvelope.Intersects(asPolyEx[j].sEnvelope) 
     940                     || !asPolyEx[i].poPolygon->Overlaps(asPolyEx[j].poPolygon) ) 
    891941            { 
    892                 relations[i][j] = NOT_RELATED; 
    893                 relations[j][i] = NOT_RELATED; 
     942 
    894943            } 
    895944            else 
    …… 
    903952                char* wkt1; 
    904953                char* wkt2; 
    905                 papoPolygons[i]->exportToWkt(&wkt1); 
    906                 papoPolygons[j]->exportToWkt(&wkt2); 
     954                asPolyEx[i].poPolygon->exportToWkt(&wkt1); 
     955                asPolyEx[j].poPolygon->exportToWkt(&wkt2); 
    907956                CPLDebug( "OGR",  
    908957                          "Bad intersection for polygons %d and %d\n" 
    …… 
    915964            } 
    916965        } 
     966 
     967        if (j < 0) 
     968        { 
     969            /* We come here because we are not included in anything */ 
     970            /* We are toplevel */ 
     971            nCountTopLevel ++; 
     972            asPolyEx[i].bIsTopLevel = TRUE; 
     973            asPolyEx[i].poEnclosingPolygon = NULL; 
     974        } 
    917975    } 
    918976 
    …… 
    933991        { 
    934992            ((OGRGeometryCollection*)geom)-> 
    935                 addGeometryDirectly( papoPolygons[i] ); 
     993                addGeometryDirectly( asPolyEx[i].poPolygon ); 
    936994        } 
    937995    } 
    …… 
    9431001    else 
    9441002    { 
    945         int* directContainerIndex = new int[nPolygonCount]; 
    946  
    947         /* Find the smallest enclosing polygon of each polygon */ 
     1003        /* STEP 3: Resort in initial order */ 
     1004        qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareByIndex); 
     1005 
     1006        /* STEP 4: Add holes as rings of their enclosing polygon */ 
    9481007        for(i=0;i<nPolygonCount;i++) 
    9491008        { 
    950             int jSmallestContainer = -1; 
    951             double areaSmallestContainer = 0; 
    952             for(j=0;j<nPolygonCount;j++) 
     1009            if (asPolyEx[i].bIsTopLevel == FALSE) 
    9531010            { 
    954                 if (i != j) 
     1011                asPolyEx[i].poEnclosingPolygon->addRing( 
     1012                    asPolyEx[i].poPolygon->getExteriorRing()); 
     1013                delete asPolyEx[i].poPolygon; 
     1014            } 
     1015            else if (nCountTopLevel == 1) 
     1016            { 
     1017                geom = asPolyEx[i].poPolygon; 
     1018            } 
     1019        } 
     1020 
     1021        /* STEP 5: Add toplevel polygons */ 
     1022        if (nCountTopLevel > 1) 
     1023        { 
     1024            for(i=0;i<nPolygonCount;i++) 
     1025            { 
     1026                if (asPolyEx[i].bIsTopLevel) 
    9551027                { 
    956                     if (relations[i][j] == IS_CONTAINED_BY) 
     1028                    if (geom == NULL) 
    9571029                    { 
    958                         if (jSmallestContainer < 0 || areas[j] < areaSmallestContainer) 
    959                         { 
    960                             jSmallestContainer = j; 
    961                             areaSmallestContainer = areas[j]; 
    962                         } 
     1030                        geom = new OGRMultiPolygon(); 
     1031                        ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon); 
     1032                    } 
     1033                    else 
     1034                    { 
     1035                        ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon); 
    9631036                    } 
    9641037                } 
    9651038            } 
    966             directContainerIndex[i] = jSmallestContainer; 
    967         } 
    968  
    969         /* Compute the inclusion depth of each polygon */ 
    970         int* containedDepth = new int [nPolygonCount]; 
    971         for(i=0;i<nPolygonCount;i++) 
    972         { 
    973             int depth = 0; 
    974             int j = directContainerIndex[i]; 
    975             while (j >= 0) 
    976             { 
    977                 j = directContainerIndex[j]; 
    978                 depth++; 
    979             } 
    980             containedDepth[i] = depth; 
    981 //          fprintf(stderr, "%d is of depth %d\n", i, depth); 
    982         } 
    983  
    984         int nbTopLevelPolygons = 0; 
    985         OGRPolygon** tempPolygons = new OGRPolygon*[nPolygonCount];  
    986  
    987         /* Create a copy of toplevel polygons */ 
    988         for(i=0;i<nPolygonCount;i++) 
    989         { 
    990             if ((containedDepth[i] % 2) == 0) 
    991             { 
    992                 nbTopLevelPolygons ++; 
    993                 tempPolygons[i] = (OGRPolygon*)papoPolygons[i]; 
    994                 papoPolygons[i] = NULL; 
    995                 if (nbTopLevelPolygons == 1) 
    996                     geom = tempPolygons[i]; 
    997             } 
    998         } 
    999  
    1000         /* Add interior rings to toplevel polygons */ 
    1001         for(i=0;i<nPolygonCount;i++) 
    1002         { 
    1003             if ((containedDepth[i] % 2) == 1) 
    1004             { 
    1005                 tempPolygons[directContainerIndex[i]]->addRing( 
    1006                     ((OGRPolygon *)papoPolygons[i])->getExteriorRing()); 
    1007                 delete papoPolygons[i]; 
    1008             } 
    1009         } 
    1010  
    1011         if (nbTopLevelPolygons > 1) 
    1012         { 
    1013             geom = new OGRMultiPolygon(); 
    1014  
    1015             /* Add toplevel polygons to the multipolygon */ 
    1016             for(i=0;i<nPolygonCount;i++) 
    1017             { 
    1018                 if ((containedDepth[i] % 2) == 0) 
    1019                 { 
    1020                     ((OGRMultiPolygon*)geom)->addGeometryDirectly( 
    1021                         tempPolygons[i]); 
    1022                     tempPolygons[i] = NULL; 
    1023                 } 
    1024             } 
    1025         } 
    1026  
    1027         delete[] tempPolygons; 
    1028         delete[] directContainerIndex; 
    1029         delete[] containedDepth; 
    1030     } 
    1031  
    1032     for(i=0;i<nPolygonCount;i++) 
    1033     { 
    1034         delete[] relations[i]; 
    1035     } 
    1036     delete[] relations; 
    1037     delete[] areas; 
    1038     delete[] envelopes; 
     1039        } 
     1040    } 
     1041 
     1042    delete[] asPolyEx; 
    10391043 
    10401044    return geom;