Changeset 14727

Show
Ignore:
Timestamp:
06/19/08 17:44:44 (6 months ago)
Author:
rouault
Message:

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

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/gdal/ogr/ogrgeometryfactory.cpp

    r14672 r14727  
    780780 */ 
    781781 
    782 enum 
    783 { 
    784     CONTAINS, 
    785     IS_CONTAINED_BY, 
    786     NOT_RELATED 
     782typedef struct _sPolyExtended sPolyExtended; 
     783 
     784struct _sPolyExtended 
     785{ 
     786    OGRPolygon*     poPolygon; 
     787    OGREnvelope     sEnvelope; 
     788    OGRLinearRing*  poExteriorRing; 
     789    OGRPoint        poAPoint; 
     790    int             nInitialIndex; 
     791    int             bIsTopLevel; 
     792    OGRPolygon*     poEnclosingPolygon; 
     793    double          dfArea; 
    787794}; 
     795 
     796static int OGRGeometryFactoryCompareArea(const void* p1, const void* p2) 
     797{ 
     798    const sPolyExtended* psPoly1 = (const sPolyExtended*) p1; 
     799    const sPolyExtended* psPoly2 = (const sPolyExtended*) p2; 
     800    if (psPoly2->dfArea < psPoly1->dfArea) 
     801        return -1; 
     802    else if (psPoly2->dfArea > psPoly1->dfArea) 
     803        return 1; 
     804    else 
     805        return 0; 
     806} 
     807 
     808static int OGRGeometryFactoryCompareByIndex(const void* p1, const void* p2) 
     809{ 
     810    const sPolyExtended* psPoly1 = (const sPolyExtended*) p1; 
     811    const sPolyExtended* psPoly2 = (const sPolyExtended*) p2; 
     812    if (psPoly1->nInitialIndex < psPoly2->nInitialIndex) 
     813        return -1; 
     814    else if (psPoly1->nInitialIndex > psPoly2->nInitialIndex) 
     815        return 1; 
     816    else 
     817        return 0; 
     818} 
     819 
    788820 
    789821OGRGeometry* OGRGeometryFactory::organizePolygons( OGRGeometry **papoPolygons, 
    …… 
    830862 
    831863/* -------------------------------------------------------------------- */ 
    832 /*      Setup per polygon relation, envelope and area information.      */ 
    833 /* -------------------------------------------------------------------- */ 
    834     OGREnvelope* envelopes = new OGREnvelope[nPolygonCount]; 
    835     int** relations = new int*[nPolygonCount]; 
    836     double* areas = new double[nPolygonCount]; 
     864/*      Setup per polygon envelope and area information.                */ 
     865/* -------------------------------------------------------------------- */ 
     866    sPolyExtended* asPolyEx = (sPolyExtended*)VSICalloc(sizeof(sPolyExtended), nPolygonCount); 
     867 
    837868    int go_on = TRUE; 
    838869    int bMixedUpGeometries = FALSE; 
    …… 
    841872    for(i=0;i<nPolygonCount;i++) 
    842873    { 
    843         relations[i] = new int[nPolygonCount]; 
    844         papoPolygons[i]->getEnvelope(&envelopes[i]); 
    845  
    846         //fprintf(stderr, "[%d] %d points\n", i, ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints()); 
     874        asPolyEx[i].nInitialIndex = i; 
     875        asPolyEx[i].poPolygon = (OGRPolygon*)papoPolygons[i]; 
     876        papoPolygons[i]->getEnvelope(&asPolyEx[i].sEnvelope); 
    847877 
    848878        if( wkbFlatten(papoPolygons[i]->getGeometryType()) == wkbPolygon 
    …… 
    850880            && ((OGRPolygon *)papoPolygons[i])->getExteriorRing()->getNumPoints() >= 4) 
    851881        { 
    852             areas[i] = ((OGRPolygon *)papoPolygons[i])->get_Area(); 
     882            asPolyEx[i].dfArea = asPolyEx[i].poPolygon->get_Area(); 
     883            asPolyEx[i].poExteriorRing = asPolyEx[i].poPolygon->getExteriorRing(); 
     884            asPolyEx[i].poExteriorRing->getPoint(0, &asPolyEx[i].poAPoint); 
    853885        } 
    854886        else 
    …… 
    866898                bNonPolygon = TRUE; 
    867899        } 
    868              
    869     } 
    870          
    871     /* This a several step algorithm : 
    872        1) Compute in a matrix how polygons relate to each other 
    873        (this is the moment for detecting pathological intersections and exiting) 
    874        2) For each polygon, find the smallest enclosing polygon 
    875        3) For each polygon, compute its inclusion depth (0 means toplevel) 
    876        4) For each polygon of odd depth (= inner ring), add it to its outer ring 
    877          
     900    } 
     901 
     902 
     903    /* This a several steps algorithm : 
     904       1) Sort polygons by descending areas 
     905       2) For each polygon of rank i, find its smallest enclosing polygon 
     906          among the polygons of rank [i-1 ... 0]. If there are no such polygon, 
     907          this is a toplevel polygon. Otherwise, depending on if the enclosing 
     908          polygon is toplevel or not, we can decide if we are toplevel or not 
     909       3) Re-sort the polygons to retrieve their inital order (nicer for some applications) 
     910       4) For each non toplevel polygon (= inner ring), add it to its outer ring 
     911       5) Add the toplevel polygons to the multipolygon 
     912 
    878913       Complexity : O(nPolygonCount^2) 
    879914    */ 
    880          
     915 
    881916    /* Compute how each polygon relate to the other ones 
    882917       To save a bit of computation we always begin the computation by a test  
    …… 
    895930    */ 
    896931 
     932    if (!bMixedUpGeometries) 
     933    { 
     934        /* STEP 1 : Sort polygons by descending area */ 
     935        qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareArea); 
     936    } 
     937    papoPolygons = NULL; /* just to use to avoid it afterwards */ 
     938 
    897939/* -------------------------------------------------------------------- */ 
    898940/*      Compute relationships, if things seem well structured.          */ 
    899941/* -------------------------------------------------------------------- */ 
    900     for(i=0; !bMixedUpGeometries && go_on && i<nPolygonCount; i++) 
    901     { 
    902         for(j=i+1; go_on && j<nPolygonCount;j++) 
    903         { 
    904             OGRPoint pointI, pointJ; 
    905             const OGRLinearRing* exteriorRingI = ((OGRPolygon *)papoPolygons[i])->getExteriorRing(); 
    906             const OGRLinearRing* exteriorRingJ = ((OGRPolygon *)papoPolygons[j])->getExteriorRing(); 
    907             exteriorRingI->getPoint(0, &pointI); 
    908             exteriorRingJ->getPoint(0, &pointJ); 
    909  
    910             if (areas[i] > areas[j] && envelopes[i].Contains(envelopes[j]) && 
    911                 ((bUseFastVersion && exteriorRingI->isPointInRing(&pointJ)) || 
    912                  (!bUseFastVersion && papoPolygons[i]->Contains(papoPolygons[j])))) 
     942 
     943    /* The first (largest) polygon is necessarily top-level */ 
     944    asPolyEx[0].bIsTopLevel = TRUE; 
     945    asPolyEx[0].poEnclosingPolygon = NULL; 
     946 
     947    int nCountTopLevel = 1; 
     948 
     949    /* STEP 2 */ 
     950    for(i=1; !bMixedUpGeometries && go_on && i<nPolygonCount; i++) 
     951    { 
     952        for(j=i-1; go_on && j>=0;j--) 
     953        { 
     954            if (asPolyEx[j].sEnvelope.Contains(asPolyEx[i].sEnvelope) && 
     955                ((bUseFastVersion && asPolyEx[j].poExteriorRing->isPointInRing(&asPolyEx[i].poAPoint)) || 
     956                 (!bUseFastVersion && asPolyEx[j].poPolygon->Contains(asPolyEx[i].poPolygon)))) 
    913957            { 
    914                 relations[i][j] = CONTAINS; 
    915                 relations[j][i] = IS_CONTAINED_BY; 
    916             } 
    917             else if (areas[j] > areas[i] && envelopes[j].Contains(envelopes[i]) && 
    918                      ((bUseFastVersion && exteriorRingJ->isPointInRing(&pointI)) || 
    919                       (!bUseFastVersion && papoPolygons[j]->Contains(papoPolygons[i])))) 
    920             { 
    921                 relations[j][i] = CONTAINS; 
    922                 relations[i][j] = IS_CONTAINED_BY; 
     958                if (asPolyEx[j].bIsTopLevel) 
     959                { 
     960                    /* We are a lake */ 
     961                    asPolyEx[i].bIsTopLevel = FALSE; 
     962                    asPolyEx[i].poEnclosingPolygon = asPolyEx[j].poPolygon; 
     963                } 
     964                else 
     965                { 
     966                    /* We are included in a something not toplevel (a lake), */ 
     967                    /* so in OGCSF we are considered as toplevel too */ 
     968                    nCountTopLevel ++; 
     969                    asPolyEx[i].bIsTopLevel = TRUE; 
     970                    asPolyEx[i].poEnclosingPolygon = NULL; 
     971                } 
     972                break; 
    923973            } 
    924974            /* We use Overlaps instead of Intersects to be more  
    925975               tolerant about touching polygons */  
    926             else if ( bUseFastVersion || !envelopes[i].Intersects(envelopes[j]) 
    927                      || !papoPolygons[i]->Overlaps(papoPolygons[j]) ) 
     976            else if ( bUseFastVersion || !asPolyEx[i].sEnvelope.Intersects(asPolyEx[j].sEnvelope) 
     977                     || !asPolyEx[i].poPolygon->Overlaps(asPolyEx[j].poPolygon) ) 
    928978            { 
    929                 relations[i][j] = NOT_RELATED; 
    930                 relations[j][i] = NOT_RELATED; 
     979 
    931980            } 
    932981            else 
    …… 
    940989                char* wkt1; 
    941990                char* wkt2; 
    942                 papoPolygons[i]->exportToWkt(&wkt1); 
    943                 papoPolygons[j]->exportToWkt(&wkt2); 
     991                asPolyEx[i].poPolygon->exportToWkt(&wkt1); 
     992                asPolyEx[j].poPolygon->exportToWkt(&wkt2); 
    944993                CPLDebug( "OGR",  
    945994                          "Bad intersection for polygons %d and %d\n" 
    …… 
    9521001            } 
    9531002        } 
     1003 
     1004        if (j < 0) 
     1005        { 
     1006            /* We come here because we are not included in anything */ 
     1007            /* We are toplevel */ 
     1008            nCountTopLevel ++; 
     1009            asPolyEx[i].bIsTopLevel = TRUE; 
     1010            asPolyEx[i].poEnclosingPolygon = NULL; 
     1011        } 
    9541012    } 
    9551013 
    …… 
    9701028        { 
    9711029            ((OGRGeometryCollection*)geom)-> 
    972                 addGeometryDirectly( papoPolygons[i] ); 
     1030                addGeometryDirectly( asPolyEx[i].poPolygon ); 
    9731031        } 
    9741032    } 
    …… 
    9801038    else 
    9811039    { 
    982         int* directContainerIndex = new int[nPolygonCount]; 
    983  
    984         /* Find the smallest enclosing polygon of each polygon */ 
     1040        /* STEP 3: Resort in initial order */ 
     1041        qsort(asPolyEx, nPolygonCount, sizeof(sPolyExtended), OGRGeometryFactoryCompareByIndex); 
     1042 
     1043        /* STEP 4: Add holes as rings of their enclosing polygon */ 
    9851044        for(i=0;i<nPolygonCount;i++) 
    9861045        { 
    987             int jSmallestContainer = -1; 
    988             double areaSmallestContainer = 0; 
    989             for(j=0;j<nPolygonCount;j++) 
     1046            if (asPolyEx[i].bIsTopLevel == FALSE) 
    9901047            { 
    991                 if (i != j) 
     1048                asPolyEx[i].poEnclosingPolygon->addRing( 
     1049                    asPolyEx[i].poPolygon->getExteriorRing()); 
     1050                delete asPolyEx[i].poPolygon; 
     1051            } 
     1052            else if (nCountTopLevel == 1) 
     1053            { 
     1054                geom = asPolyEx[i].poPolygon; 
     1055            } 
     1056        } 
     1057 
     1058        /* STEP 5: Add toplevel polygons */ 
     1059        if (nCountTopLevel > 1) 
     1060        { 
     1061            for(i=0;i<nPolygonCount;i++) 
     1062            { 
     1063                if (asPolyEx[i].bIsTopLevel) 
    9921064                { 
    993                     if (relations[i][j] == IS_CONTAINED_BY) 
     1065                    if (geom == NULL) 
    9941066                    { 
    995                         if (jSmallestContainer < 0 || areas[j] < areaSmallestContainer) 
    996                         { 
    997                             jSmallestContainer = j; 
    998                             areaSmallestContainer = areas[j]; 
    999                         } 
     1067                        geom = new OGRMultiPolygon(); 
     1068                        ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon); 
     1069                    } 
     1070                    else 
     1071                    { 
     1072                        ((OGRMultiPolygon*)geom)->addGeometryDirectly(asPolyEx[i].poPolygon); 
    10001073                    } 
    10011074                } 
    10021075            } 
    1003             directContainerIndex[i] = jSmallestContainer; 
    1004         } 
    1005  
    1006         /* Compute the inclusion depth of each polygon */ 
    1007         int* containedDepth = new int [nPolygonCount]; 
    1008         for(i=0;i<nPolygonCount;i++) 
    1009         { 
    1010             int depth = 0; 
    1011             int j = directContainerIndex[i]; 
    1012             while (j >= 0) 
    1013             { 
    1014                 j = directContainerIndex[j]; 
    1015                 depth++; 
    1016             } 
    1017             containedDepth[i] = depth; 
    1018 //          fprintf(stderr, "%d is of depth %d\n", i, depth); 
    1019         } 
    1020  
    1021         int nbTopLevelPolygons = 0; 
    1022         OGRPolygon** tempPolygons = new OGRPolygon*[nPolygonCount];  
    1023  
    1024         /* Create a copy of toplevel polygons */ 
    1025         for(i=0;i<nPolygonCount;i++) 
    1026         { 
    1027             if ((containedDepth[i] % 2) == 0) 
    1028             { 
    1029                 nbTopLevelPolygons ++; 
    1030                 tempPolygons[i] = (OGRPolygon*)papoPolygons[i]; 
    1031                 papoPolygons[i] = NULL; 
    1032                 if (nbTopLevelPolygons == 1) 
    1033                     geom = tempPolygons[i]; 
    1034             } 
    1035         } 
    1036  
    1037         /* Add interior rings to toplevel polygons */ 
    1038         for(i=0;i<nPolygonCount;i++) 
    1039         { 
    1040             if ((containedDepth[i] % 2) == 1) 
    1041             { 
    1042                 tempPolygons[directContainerIndex[i]]->addRing( 
    1043                     ((OGRPolygon *)papoPolygons[i])->getExteriorRing()); 
    1044                 delete papoPolygons[i]; 
    1045             } 
    1046         } 
    1047  
    1048         if (nbTopLevelPolygons > 1) 
    1049         { 
    1050             geom = new OGRMultiPolygon(); 
    1051  
    1052             /* Add toplevel polygons to the multipolygon */ 
    1053             for(i=0;i<nPolygonCount;i++) 
    1054             { 
    1055                 if ((containedDepth[i] % 2) == 0) 
    1056                 { 
    1057                     ((OGRMultiPolygon*)geom)->addGeometryDirectly( 
    1058                         tempPolygons[i]); 
    1059                     tempPolygons[i] = NULL; 
    1060                 } 
    1061             } 
    1062         } 
    1063  
    1064         delete[] tempPolygons; 
    1065         delete[] directContainerIndex; 
    1066         delete[] containedDepth; 
    1067     } 
    1068  
    1069     for(i=0;i<nPolygonCount;i++) 
    1070     { 
    1071         delete[] relations[i]; 
    1072     } 
    1073     delete[] relations; 
    1074     delete[] areas; 
    1075     delete[] envelopes; 
     1076        } 
     1077    } 
     1078 
     1079    CPLFree(asPolyEx); 
    10761080 
    10771081    return geom;