Changeset 13599

Show
Ignore:
Timestamp:
01/26/08 10:33:38 (6 months ago)
Author:
rouault
Message:

Cleanup to remove unused old classification code for multipolygons (#2174)

Files:

Legend:

Unmodified
Added
Removed
Modified
Copied
Moved

  • trunk/gdal/ogr/ogrsf_frmts/shape/shape2ogr.cpp

    r13551 r13599  
    4545} 
    4646 
    47 /* A helper class which represents 2D bounding box */ 
    48 class RingExtent 
    49 { 
    50 public: 
    51     RingExtent(); 
    52     /* Calculates bounding box for given set of points */ 
    53     void calculate( int ringParts, double *ringX, double *ringY ); 
    54     /* Does this extent contains rhs*/ 
    55     bool contains( const RingExtent &rhs ); 
    56     /* Does this extent contains given point*/ 
    57     bool contains( double x, double y ); 
    58 private: 
    59     bool empty; /* extent is empty */ 
    60     double xMin, yMin; /* lower left vertex of the bounding box */ 
    61     double xMax, yMax; /* upper right vertex of the bounding box */ 
    62 }; 
    63  
    64 RingExtent::RingExtent(): 
    65     empty( true ), 
    66     xMin( 0 ), //nan ? 
    67     yMin( 0 ),  
    68     xMax( 0 ),  
    69     yMax( 0 ) 
    70 { 
    71 } 
    72  
    73 void RingExtent::calculate( int ringParts, double *ringX, double *ringY ) 
    74 { 
    75     empty = ringParts <= 0; 
    76     if ( !empty ) 
    77     { 
    78         xMin = xMax = *ringX; 
    79         yMin = yMax = *ringY; 
    80         while( --ringParts > 0 ) 
    81         { 
    82             if ( xMin > *ringX ) 
    83                 xMin = *ringX; 
    84             else if ( xMax < *ringX ) 
    85                 xMax = *ringX; 
    86             if ( yMin > *ringY ) 
    87                yMin = *ringY; 
    88             else if ( yMax < *ringY ) 
    89                 yMax = *ringY; 
    90             ++ringX;  
    91             ++ringY;  
    92         } 
    93     } 
    94 } 
    95  
    96 inline bool RingExtent::contains( double x, double y ) 
    97 { 
    98     if ( empty ) 
    99         return false; 
    100  
    101     return xMin <= x && x <= xMax && 
    102         yMin <= y && y <= yMax; 
    103 } 
    104  
    105 inline bool RingExtent::contains( const RingExtent &extent ) 
    106 { 
    107     return contains( extent.xMin, extent.yMin ) && 
    108         contains( extent.xMax, extent.yMax ); 
    109 } 
    110  
    111  
    11247/************************************************************************/ 
    11348/*                        RingStartEnd                                  */ 
    …… 
    13065        *start = psShape->panPartStart[ring]; 
    13166    } 
    132 } 
    133      
    134 /************************************************************************/ 
    135 /*                          RingDirection()                             */ 
    136 /*                                                                      */ 
    137 /*      return: 1 CW (clockwise)                                        */ 
    138 /*             -1 CCW (counterclockwise)                                */ 
    139 /*              0 error                                                 */ 
    140 /************************************************************************/ 
    141  
    142 int RingDirection ( SHPObject *Shape, int ring )  
    143 { 
    144     int    i, start, end, v, next; 
    145     double *x, *y, dx0, dy0, dx1, dy1, area2;  
    146  
    147     /* Note: first and last coordinate MUST be identical according to shapefile specification */ 
    148     if ( ring < 0 || ring >= Shape->nParts ) return ( 0 ); 
    149  
    150     x = Shape->padfX; 
    151     y = Shape->padfY; 
    152  
    153     RingStartEnd ( Shape, ring, &start, &end ); 
    154  
    155     if ( start == end ) return ( 0 ); // empty ring!?!? 
    156  
    157     /* Find the lowest rightmost vertex */ 
    158     v = start;   
    159     for ( i = start + 1; i < end; i++ ) 
    160     { 
    161         /* => v < end */ 
    162         if ( y[i] < y[v] || ( y[i] == y[v] && x[i] > x[v] ) ) 
    163         { 
    164             v = i; 
    165         } 
    166     } 
    167      
    168     /* Vertices may be duplicate, we have to go to nearest different in each direction */ 
    169     /* preceding */ 
    170     next = v - 1; 
    171     while ( 1 ) 
    172     { 
    173         if ( next < start )  
    174         { 
    175             next = end - 1;  
    176         } 
    177  
    178         if( !epsilonEqual(x[next], x[v], EPSILON)  
    179             || !epsilonEqual(y[next], y[v], EPSILON) ) 
    180         { 
    181             break; 
    182         } 
    183  
    184         if ( next == v ) /* So we cannot get into endless loop */ 
    185         { 
    186             break; 
    187         } 
    188  
    189         next--; 
    190     } 
    191              
    192     dx0 = x[next] - x[v]; 
    193     dy0 = y[next] - y[v]; 
    194      
    195      
    196     /* following */ 
    197     next = v + 1; 
    198     while ( 1 ) 
    199     { 
    200         if ( next >= end )  
    201         { 
    202             next = start;  
    203         } 
    204  
    205         if ( !epsilonEqual(x[next], x[v], EPSILON)  
    206              || !epsilonEqual(y[next], y[v], EPSILON) ) 
    207         { 
    208             break; 
    209         } 
    210  
    211         if ( next == v ) /* So we cannot get into endless loop */ 
    212         { 
    213             break; 
    214         } 
    215  
    216         next++; 
    217     } 
    218  
    219     dx1 = x[next] - x[v]; 
    220     dy1 = y[next] - y[v]; 
    221  
    222     area2 = dx1 * dy0 - dx0 * dy1; 
    223      
    224     if ( area2 > 0 )      /* CCW */ 
    225         return -1; 
    226     else if ( area2 < 0 )  /* CW */ 
    227         return 1; 
    228      
    229     return 0; /* error */ 
    230 } 
    231  
    232 /************************************************************************/ 
    233 /*                          PointInRing()                               */ 
    234 /*                                                                      */ 
    235 /*      return: 1 point is inside the ring                              */ 
    236 /*              0 point is outside the ring                             */ 
    237 /*                                                                      */ 
    238 /*              for point exactly on the boundary it returns 0 or 1     */ 
    239 /************************************************************************/ 
    240  
    241 int PointInRing ( SHPObject *Shape, int ring, double x, double y )  
    242 { 
    243     int    i, start, end, c = 0; 
    244     double *xp, *yp;  
    245  
    246     if ( ring < 0 || ring >= Shape->nParts ) return ( 0 ); 
    247  
    248     xp = Shape->padfX; 
    249     yp = Shape->padfY; 
    250  
    251     RingStartEnd ( Shape, ring, &start, &end ); 
    252  
    253     /* This code was originaly written by Randolph Franklin, here it is slightly modified. */ 
    254     for (i = start; i < end; i++) { 
    255         if ( ( ( ( yp[i] <= y ) && ( y < yp[i+1] ) ) ||  
    256                ( ( yp[i+1] <= y ) && ( y < yp[i] ) ) ) && 
    257              ( x < (xp[i+1] - xp[i]) * (y - yp[i]) / (yp[i+1] - yp[i]) + xp[i] ) 
    258              )  
    259         { 
    260             c = !c; 
    261         } 
    262     } 
    263     return c; 
    264 } 
    265  
    266 /************************************************************************/ 
    267 /*                          RingInRing()                                */ 
    268 /*                                                                      */ 
    269 /* Checks point by point using PointInRing if oRing contains iRing      */ 
    270 /*                                                                      */ 
    271 /*      return: 1 oRing contains iRing                                  */ 
    272 /*              0 iRing is not contained by oRing                       */ 
    273 /*                                                                      */ 
    274 /*              for point exactly on the boundary it returns 0 or 1     */ 
    275 /************************************************************************/ 
    276 int RingInRing( SHPObject *Shape, int oRing, int iRing ) 
    277 { 
    278     int iRingStart, iRingEnd; 
    279     RingStartEnd ( Shape, iRing, &iRingStart, &iRingEnd ); 
    280     while( iRingStart < iRingEnd ) 
    281     { 
    282         if ( !PointInRing( Shape, oRing, Shape->padfX[iRingStart], Shape->padfY[iRingStart] ) ) 
    283             return 0; 
    284         ++iRingStart; 
    285     } 
    286     return 1; 
    28767} 
    28868     
    …… 
    461241        } 
    462242 
    463         // OGRGeometryFactory::organizePolygons does not require GOES anymore  
    464         else if( TRUE /* OGRGeometryFactory::haveGEOS()*/ ) 
     243        else 
    465244        { 
    466245            OGRPolygon** tabPolygons = new OGRPolygon*[psShape->nParts]; 
    …… 
    485264            delete[] tabPolygons; 
    486265        } 
    487         else  
    488         { 
    489             /* Multipart polygon. */ 
    490  
    491             int nOuter = 0;         /* Number of outer rings */  
    492             int *direction = NULL;  /* ring direction (1 CW,outer, -1 CCW, inner) */  
    493             int *outer = NULL;      /* list of outer rings */  
    494                  
    495             /* Outer ring index for inner rings, -1 for outer rings */  
    496             int *outside = NULL;    
    497  
    498             direction = (int *) CPLMalloc( psShape->nParts * sizeof(int) ); 
    499             outer = (int *) CPLMalloc( psShape->nParts * sizeof(int) ); 
    500             outside = (int *) CPLMalloc( psShape->nParts * sizeof(int) ); 
    501              
    502             /* First distinguish outer from inner rings */ 
    503             for( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    504             { 
    505                 direction[iRing] = RingDirection ( psShape, iRing); 
    506                 if ( direction[iRing] == 1 ) 
    507                 { 
    508                     outer[nOuter] = iRing; 
    509                     nOuter++; 
    510                 } 
    511                 outside[iRing] = -1; 
    512             } 
    513  
    514             if ( nOuter == 1 ) 
    515             {  
    516                 /* One outer ring? 
    517                  * we do not need to check anything as we assume that shapefile is valid 
    518                  * ie. if there is one outer ring then all inner rings are inside. 
    519                  */ 
    520                 for( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    521                 { 
    522                     if ( direction[iRing] == -1 ) 
    523                     { 
    524                         outside[iRing] = outer[ 0 ]; 
    525                     } 
    526                 } 
    527             } 
    528             else 
    529             { 
    530                 /* Calculate ring extents */ 
    531                 RingExtent *extents = new RingExtent[ psShape->nParts ]; 
    532  
    533                 for( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    534                 { 
    535                     int start = 0; 
    536                     int end = 0; 
    537                     RingStartEnd ( psShape, iRing, &start, &end ); 
    538                     extents[ iRing ].calculate( end-start, psShape->padfX + start, psShape->padfY + start ); 
    539                 } 
    540  
    541                 /* Outer ring index */ 
    542                 int oRing; 
    543  
    544                 /* This loops tries to assign inner to outer rings. 
    545                  * The fundamental assumption is that if the extent of an inner ring is contained 
    546                  * in the extent of exactly one outer ring then this inner ring is inside the outer ring. 
    547                  * Otherwise that shapefile is broken. 
    548                  * 
    549                  * XXX - mloskot - This loop does incorrect classification, see Bug 1217 
    550                  */ 
    551                 for( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    552                 { 
    553                     if ( direction[ iRing ] != 1 ) 
    554                     { 
    555                         /* Inner ring */ 
    556  
    557                         /* Find the outer ring for iRing */ 
    558                         for( oRing = 0; oRing < nOuter; oRing++ ) 
    559                         { 
    560                             if ( extents[ outer[ oRing ] ].contains( extents[ iRing ] ) ) 
    561                             { 
    562                                 /* The outer ring contains the inner ring, we have to execute 
    563                                  * some additional tests. 
    564                                  */  
    565                                 if ( outside[ iRing ] == -1 ) 
    566                                 { 
    567                                     /* this is the first outer ring, assume it containes this inner iRing */ 
    568                                     outside[ iRing ] = outer[ oRing ]; 
    569                                 } 
    570                                 else 
    571                                 {  
    572                                     /* This is another outer ring, which contains iRing, have to distinguish 
    573                                      * between this and previous outer rings using RingInPoint. 
    574                                      * Here is a place for some optimization as sometimes we may check 
    575                                      * the same ring many times. 
    576                                      */ 
    577                                     if ( !RingInRing( psShape, outside[ iRing ], iRing ) ) 
    578                                     { 
    579                                         outside[ iRing ] = outer[ oRing ]; 
    580                                     } 
    581                                 }                             
    582                             } 
    583                         } 
    584                     } 
    585                 } 
    586  
    587                 /* Destroy ring extents object. */ 
    588                 delete[] extents; 
    589             } 
    590  
    591             /* The inner rings are preliminary sorted, but: 
    592              * - Some of them are not assigned to any outer ring. Probably broken shapefile 
    593              * - Some inner rings are assigned to outer rings using extents only. 
    594              * Additional geometry tests are ommited here due to perfomance penalties. 
    595              * Promote any unassigned inner rings to be outside rings. 
    596              */ 
    597             for ( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    598             { 
    599                 /* Outer rings */ 
    600                 if( direction[iRing] != 1 && outside[iRing] == -1 ) 
    601                 { 
    602                     direction[iRing] = 1; /* this isn't exactly true! */ 
    603                     outer[nOuter++] = iRing; 
    604                 } 
    605             } 
    606  
    607             if ( nOuter == 1 ) 
    608             { 
    609                 /* One outer ring and more inner rings */ 
    610                 OGRPolygon    *poOGRPoly = NULL; 
    611                 OGRLinearRing *poRing = NULL; 
    612  
    613                 /* Outer */ 
    614                 poOGR = poOGRPoly = new OGRPolygon(); 
    615                 poRing = CreateLinearRing ( psShape, outer[0] ); 
    616                 poOGRPoly->addRingDirectly( poRing ); 
    617  
    618                 /* Inner */ 
    619                 for( iRing = 0; iRing < psShape->nParts; iRing++ ) 
    620                 { 
    621                     if ( direction[iRing] == -1 ) 
    622                     { 
    623                         poRing = CreateLinearRing ( psShape, iRing ); 
    624                         poOGRPoly->addRingDirectly( poRing ); 
    625                     } 
    626                 } 
    627             } 
    628             else 
    629             {  
    630                 OGRGeometryCollection *poOGRGeCo = NULL; 
    631  
    632                 poOGR = poOGRGeCo = new OGRMultiPolygon(); 
    633  
    634                 for ( iRing = 0; iRing < nOuter; iRing++ ) 
    635                 { 
    636                     /* Outer rings */ 
    637                     int oRing;  
    638  
    639                     OGRPolygon    *poOGRPoly = NULL; 
    640                     OGRLinearRing *poRing = NULL; 
    641  
    642                     oRing = outer[iRing]; 
    643  
    644                     /* Outer */ 
    645                     poOGRPoly = new OGRPolygon(); 
    646                     poRing = CreateLinearRing ( psShape, oRing ); 
    647                     poOGRPoly->addRingDirectly( poRing ); 
    648  
    649                     /* Inner */ 
    650                     for( int iRing2 = 0; iRing2 < psShape->nParts; iRing2++ ) 
    651                     { 
    652                         if ( outside[iRing2] == oRing ) 
    653                         { 
    654                             poRing = CreateLinearRing ( psShape, iRing2 ); 
    655                             poOGRPoly->addRingDirectly( poRing ); 
    656                         } 
    657                     } 
    658                      
    659                     poOGRGeCo->addGeometryDirectly (poOGRPoly); 
    660                 } 
    661             } 
    662  
    663             CPLFree( direction ); 
    664             CPLFree( outer ); 
    665             CPLFree( outside ); 
    666         } /* End of multipart polygon processing. */ 
    667266 
    668267        if( psShape->nSHPType == SHPT_POLYGON )