LCOV - code coverage report
Current view: directory - alg - polygonize.cpp (source / functions) Found Hit Coverage
Test: gdal_filtered.info Lines: 232 206 88.8 %
Date: 2012-12-26 Functions: 10 9 90.0 %

       1                 : /******************************************************************************
       2                 :  * $Id: polygonize.cpp 22501 2011-06-04 21:28:47Z rouault $
       3                 :  * Project:  GDAL
       4                 :  * Purpose:  Raster to Polygon Converter
       5                 :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6                 :  *
       7                 :  ******************************************************************************
       8                 :  * Copyright (c) 2008, Frank Warmerdam
       9                 :  *
      10                 :  * Permission is hereby granted, free of charge, to any person obtaining a
      11                 :  * copy of this software and associated documentation files (the "Software"),
      12                 :  * to deal in the Software without restriction, including without limitation
      13                 :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      14                 :  * and/or sell copies of the Software, and to permit persons to whom the
      15                 :  * Software is furnished to do so, subject to the following conditions:
      16                 :  *
      17                 :  * The above copyright notice and this permission notice shall be included
      18                 :  * in all copies or substantial portions of the Software.
      19                 :  *
      20                 :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      21                 :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      22                 :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      23                 :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      24                 :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      25                 :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      26                 :  * DEALINGS IN THE SOFTWARE.
      27                 :  ****************************************************************************/
      28                 : 
      29                 : #include "gdal_alg_priv.h"
      30                 : #include "cpl_conv.h"
      31                 : #include "cpl_string.h"
      32                 : #include <vector>
      33                 : 
      34                 : CPL_CVSID("$Id: polygonize.cpp 22501 2011-06-04 21:28:47Z rouault $");
      35                 : 
      36                 : #define GP_NODATA_MARKER -51502112
      37                 : 
      38                 : #ifdef OGR_ENABLED
      39                 : 
      40                 : /************************************************************************/
      41                 : /* ==================================================================== */
      42                 : /*                               RPolygon                               */
      43                 : /*                  */
      44                 : /*      This is a helper class to hold polygons while they are being    */
      45                 : /*      formed in memory, and to provide services to coalesce a much    */
      46                 : /*      of edge sections into complete rings.                           */
      47                 : /* ==================================================================== */
      48                 : /************************************************************************/
      49                 : 
      50             209 : class RPolygon {
      51                 : public:
      52             209 :     RPolygon( int nValue ) { nPolyValue = nValue; nLastLineUpdated = -1; }
      53                 : 
      54                 :     int              nPolyValue;
      55                 :     int              nLastLineUpdated;
      56                 : 
      57                 :     std::vector< std::vector<int> > aanXY;
      58                 : 
      59                 :     void             AddSegment( int x1, int y1, int x2, int y2 );
      60                 :     void             Dump();
      61                 :     void             Coalesce();
      62                 :     void             Merge( int iBaseString, int iSrcString, int iDirection );
      63                 : };
      64                 : 
      65                 : /************************************************************************/
      66                 : /*                                Dump()                                */
      67                 : /************************************************************************/
      68               0 : void RPolygon::Dump()
      69                 : {
      70                 :     size_t iString;
      71                 : 
      72                 :     printf( "RPolygon: Value=%d, LastLineUpdated=%d\n",
      73               0 :             nPolyValue, nLastLineUpdated );
      74                 :     
      75               0 :     for( iString = 0; iString < aanXY.size(); iString++ )
      76                 :     {
      77               0 :         std::vector<int> &anString = aanXY[iString];
      78                 :         size_t iVert;
      79                 :      
      80               0 :         printf( "  String %d:\n", (int) iString );
      81               0 :         for( iVert = 0; iVert < anString.size(); iVert += 2 )
      82                 :         {
      83               0 :             printf( "    (%d,%d)\n", anString[iVert], anString[iVert+1] );
      84                 :         }
      85                 :     }
      86               0 : }
      87                 : 
      88                 : /************************************************************************/
      89                 : /*                              Coalesce()                              */
      90                 : /************************************************************************/
      91                 : 
      92             201 : void RPolygon::Coalesce()
      93                 : 
      94                 : {
      95                 :     size_t iBaseString;
      96                 : 
      97                 : /* -------------------------------------------------------------------- */
      98                 : /*      Iterate over loops starting from the first, trying to merge     */
      99                 : /*      other segments into them.                                       */
     100                 : /* -------------------------------------------------------------------- */
     101             411 :     for( iBaseString = 0; iBaseString < aanXY.size(); iBaseString++ )
     102                 :     {
     103             210 :         std::vector<int> &anBase = aanXY[iBaseString];
     104             210 :         int bMergeHappened = TRUE;
     105                 : 
     106                 : /* -------------------------------------------------------------------- */
     107                 : /*      Keep trying to merge the following strings into our target      */
     108                 : /*      "base" string till we have tried them all once without any      */
     109                 : /*      mergers.                                                        */
     110                 : /* -------------------------------------------------------------------- */
     111             890 :         while( bMergeHappened )
     112                 :         {
     113                 :             size_t iString;
     114                 : 
     115             470 :             bMergeHappened = FALSE;
     116                 : 
     117                 : /* -------------------------------------------------------------------- */
     118                 : /*      Loop over the following strings, trying to find one we can      */
     119                 : /*      merge onto the end of our base string.                          */
     120                 : /* -------------------------------------------------------------------- */
     121            1077 :             for( iString = iBaseString+1; 
     122                 :                  iString < aanXY.size(); 
     123                 :                  iString++ )
     124                 :             {
     125             607 :                 std::vector<int> &anString = aanXY[iString];
     126                 : 
     127             607 :                 if( anBase[anBase.size()-2] == anString[0]
     128                 :                     && anBase[anBase.size()-1] == anString[1] )
     129                 :                 {
     130              16 :                     Merge( iBaseString, iString, 1 );
     131              16 :                     bMergeHappened = TRUE;
     132                 :                 }
     133             591 :                 else if( anBase[anBase.size()-2] == anString[anString.size()-2]
     134                 :                          && anBase[anBase.size()-1] == anString[anString.size()-1] )
     135                 :                 {
     136             246 :                     Merge( iBaseString, iString, -1 );
     137             246 :                     bMergeHappened = TRUE;
     138                 :                 }
     139                 :             }
     140                 :         }
     141                 : 
     142                 :         /* At this point our loop *should* be closed! */
     143                 : 
     144                 :         CPLAssert( anBase[0] == anBase[anBase.size()-2]
     145             210 :                    && anBase[1] == anBase[anBase.size()-1] );
     146                 :     }
     147                 : 
     148             201 : }
     149                 : 
     150                 : /************************************************************************/
     151                 : /*                               Merge()                                */
     152                 : /************************************************************************/
     153                 : 
     154             262 : void RPolygon::Merge( int iBaseString, int iSrcString, int iDirection )
     155                 : 
     156                 : {
     157             262 :     std::vector<int> &anBase = aanXY[iBaseString];
     158             262 :     std::vector<int> &anString = aanXY[iSrcString];
     159                 :     int iStart, iEnd, i;
     160                 : 
     161             262 :     if( iDirection == 1 )
     162                 :     {
     163              16 :         iStart = 1;
     164              16 :         iEnd = anString.size() / 2; 
     165                 :     }
     166                 :     else
     167                 :     {
     168             246 :         iStart = anString.size() / 2 - 2;
     169             246 :         iEnd = -1; 
     170                 :     }
     171                 :     
     172             840 :     for( i = iStart; i != iEnd; i += iDirection )
     173                 :     {
     174             578 :         anBase.push_back( anString[i*2+0] );
     175             578 :         anBase.push_back( anString[i*2+1] );
     176                 :     }
     177                 :     
     178             262 :     if( iSrcString < ((int) aanXY.size())-1 )
     179              36 :         aanXY[iSrcString] = aanXY[aanXY.size()-1];
     180                 : 
     181             262 :     size_t nSize = aanXY.size(); 
     182             262 :     aanXY.resize(nSize-1);
     183             262 : }
     184                 : 
     185                 : /************************************************************************/
     186                 : /*                             AddSegment()                             */
     187                 : /************************************************************************/
     188                 : 
     189            1634 : void RPolygon::AddSegment( int x1, int y1, int x2, int y2 )
     190                 : 
     191                 : {
     192            1634 :     nLastLineUpdated = MAX(y1, y2);
     193                 : 
     194                 : /* -------------------------------------------------------------------- */
     195                 : /*      Is there an existing string ending with this?                   */
     196                 : /* -------------------------------------------------------------------- */
     197                 :     size_t iString;
     198                 : 
     199            4434 :     for( iString = 0; iString < aanXY.size(); iString++ )
     200                 :     {
     201            3942 :         std::vector<int> &anString = aanXY[iString];
     202            3942 :         size_t nSSize = anString.size();
     203                 :         
     204            3942 :         if( anString[nSSize-2] == x1 
     205                 :             && anString[nSSize-1] == y1 )
     206                 :         {
     207                 :             int nTemp;
     208                 : 
     209            1094 :             nTemp = x2;
     210            1094 :             x2 = x1;
     211            1094 :             x1 = nTemp;
     212                 : 
     213            1094 :             nTemp = y2;
     214            1094 :             y2 = y1;
     215            1094 :             y1 = nTemp;
     216                 :         }
     217                 : 
     218            3942 :         if( anString[nSSize-2] == x2 
     219                 :             && anString[nSSize-1] == y2 )
     220                 :         {
     221                 :             // We are going to add a segment, but should we just extend 
     222                 :             // an existing segment already going in the right direction?
     223                 : 
     224            1142 :             int nLastLen = MAX(ABS(anString[nSSize-4]-anString[nSSize-2]),
     225                 :                                ABS(anString[nSSize-3]-anString[nSSize-1]));
     226                 :             
     227            1142 :             if( nSSize >= 4 
     228                 :                 && (anString[nSSize-4] - anString[nSSize-2]
     229                 :                     == (anString[nSSize-2] - x1)*nLastLen)
     230                 :                 && (anString[nSSize-3] - anString[nSSize-1]
     231                 :                     == (anString[nSSize-1] - y1)*nLastLen) )
     232                 :             {
     233             588 :                 anString.pop_back();
     234             588 :                 anString.pop_back();
     235                 :             }
     236                 : 
     237            1142 :             anString.push_back( x1 );
     238            1142 :             anString.push_back( y1 );
     239            1142 :             return;
     240                 :         }
     241                 :     }
     242                 : 
     243                 : /* -------------------------------------------------------------------- */
     244                 : /*      Create a new string.                                            */
     245                 : /* -------------------------------------------------------------------- */
     246             492 :     size_t nSize = aanXY.size();
     247             492 :     aanXY.resize(nSize + 1);
     248             492 :     std::vector<int> &anString = aanXY[nSize];
     249                 : 
     250             492 :     anString.push_back( x1 );
     251             492 :     anString.push_back( y1 );
     252             492 :     anString.push_back( x2 );
     253             492 :     anString.push_back( y2 );
     254                 : 
     255             492 :     return;
     256                 : }
     257                 : 
     258                 : /************************************************************************/
     259                 : /* ==================================================================== */
     260                 : /*     End of RPolygon                                                  */
     261                 : /* ==================================================================== */
     262                 : /************************************************************************/
     263                 : 
     264                 : /************************************************************************/
     265                 : /*                              AddEdges()                              */
     266                 : /*                                                                      */
     267                 : /*      Examine one pixel and compare to its neighbour above            */
     268                 : /*      (previous) and right.  If they are different polygon ids        */
     269                 : /*      then add the pixel edge to this polygon and the one on the      */
     270                 : /*      other side of the edge.                                         */
     271                 : /************************************************************************/
     272                 : 
     273            1831 : static void AddEdges( GInt32 *panThisLineId, GInt32 *panLastLineId, 
     274                 :                       GInt32 *panPolyIdMap, GInt32 *panPolyValue,
     275                 :                       RPolygon **papoPoly, int iX, int iY )
     276                 : 
     277                 : {
     278            1831 :     int nThisId = panThisLineId[iX];
     279            1831 :     int nRightId = panThisLineId[iX+1];
     280            1831 :     int nPreviousId = panLastLineId[iX];
     281            1831 :     int iXReal = iX - 1;
     282                 : 
     283            1831 :     if( nThisId != -1 )
     284            1687 :         nThisId = panPolyIdMap[nThisId];
     285            1831 :     if( nRightId != -1 )
     286            1687 :         nRightId = panPolyIdMap[nRightId];
     287            1831 :     if( nPreviousId != -1 )
     288            1687 :         nPreviousId = panPolyIdMap[nPreviousId];
     289                 : 
     290            1831 :     if( nThisId != nPreviousId )
     291                 :     {
     292             461 :         if( nThisId != -1 )
     293                 :         {
     294             400 :             if( papoPoly[nThisId] == NULL )
     295               0 :                 papoPoly[nThisId] = new RPolygon( panPolyValue[nThisId] );
     296                 : 
     297             400 :             papoPoly[nThisId]->AddSegment( iXReal, iY, iXReal+1, iY );
     298                 :         }
     299             461 :         if( nPreviousId != -1 )
     300                 :         {
     301             400 :             if( papoPoly[nPreviousId] == NULL )
     302               0 :                 papoPoly[nPreviousId] = new RPolygon(panPolyValue[nPreviousId]);
     303                 : 
     304             400 :             papoPoly[nPreviousId]->AddSegment( iXReal, iY, iXReal+1, iY );
     305                 :         }
     306                 :     }
     307                 : 
     308            1831 :     if( nThisId != nRightId )
     309                 :     {
     310             494 :         if( nThisId != -1 )
     311                 :         {
     312             417 :             if( papoPoly[nThisId] == NULL )
     313               0 :                 papoPoly[nThisId] = new RPolygon(panPolyValue[nThisId]);
     314                 : 
     315             417 :             papoPoly[nThisId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1 );
     316                 :         }
     317                 : 
     318             494 :         if( nRightId != -1 )
     319                 :         {
     320             417 :             if( papoPoly[nRightId] == NULL )
     321             209 :                 papoPoly[nRightId] = new RPolygon(panPolyValue[nRightId]);
     322                 : 
     323             417 :             papoPoly[nRightId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1 );
     324                 :         }
     325                 :     }
     326            1831 : }
     327                 : 
     328                 : /************************************************************************/
     329                 : /*                         EmitPolygonToLayer()                         */
     330                 : /************************************************************************/
     331                 : 
     332                 : static CPLErr
     333             201 : EmitPolygonToLayer( OGRLayerH hOutLayer, int iPixValField,
     334                 :                     RPolygon *poRPoly, double *padfGeoTransform )
     335                 : 
     336                 : {
     337                 :     OGRFeatureH hFeat;
     338                 :     OGRGeometryH hPolygon;
     339                 : 
     340                 : /* -------------------------------------------------------------------- */
     341                 : /*      Turn bits of lines into coherent rings.                         */
     342                 : /* -------------------------------------------------------------------- */
     343             201 :     poRPoly->Coalesce();
     344                 : 
     345                 : /* -------------------------------------------------------------------- */
     346                 : /*      Create the polygon geometry.                                    */
     347                 : /* -------------------------------------------------------------------- */
     348                 :     size_t iString;
     349                 : 
     350             201 :     hPolygon = OGR_G_CreateGeometry( wkbPolygon );
     351                 :     
     352             411 :     for( iString = 0; iString < poRPoly->aanXY.size(); iString++ )
     353                 :     {
     354             210 :         std::vector<int> &anString = poRPoly->aanXY[iString];
     355             210 :         OGRGeometryH hRing = OGR_G_CreateGeometry( wkbLinearRing );
     356                 : 
     357                 :         int iVert;
     358                 : 
     359                 :         // we go last to first to ensure the linestring is allocated to 
     360                 :         // the proper size on the first try.
     361            1426 :         for( iVert = anString.size()/2 - 1; iVert >= 0; iVert-- )
     362                 :         {
     363                 :             double dfX, dfY;
     364                 :             int    nPixelX, nPixelY;
     365                 :             
     366            1216 :             nPixelX = anString[iVert*2];
     367            1216 :             nPixelY = anString[iVert*2+1];
     368                 : 
     369            1216 :             dfX = padfGeoTransform[0] 
     370            1216 :                 + nPixelX * padfGeoTransform[1]
     371            2432 :                 + nPixelY * padfGeoTransform[2];
     372            1216 :             dfY = padfGeoTransform[3] 
     373            1216 :                 + nPixelX * padfGeoTransform[4]
     374            2432 :                 + nPixelY * padfGeoTransform[5];
     375                 : 
     376            1216 :             OGR_G_SetPoint_2D( hRing, iVert, dfX, dfY );
     377                 :         }
     378                 : 
     379             210 :         OGR_G_AddGeometryDirectly( hPolygon, hRing );
     380                 :     }
     381                 :     
     382                 : /* -------------------------------------------------------------------- */
     383                 : /*      Create the feature object.                                      */
     384                 : /* -------------------------------------------------------------------- */
     385             201 :     hFeat = OGR_F_Create( OGR_L_GetLayerDefn( hOutLayer ) );
     386                 : 
     387             201 :     OGR_F_SetGeometryDirectly( hFeat, hPolygon );
     388                 : 
     389             201 :     if( iPixValField >= 0 )
     390             201 :         OGR_F_SetFieldInteger( hFeat, iPixValField, poRPoly->nPolyValue );
     391                 : 
     392                 : /* -------------------------------------------------------------------- */
     393                 : /*      Write the to the layer.                                         */
     394                 : /* -------------------------------------------------------------------- */
     395             201 :     CPLErr eErr = CE_None;
     396                 : 
     397             201 :     if( OGR_L_CreateFeature( hOutLayer, hFeat ) != OGRERR_NONE )
     398               0 :         eErr = CE_Failure;
     399                 : 
     400             201 :     OGR_F_Destroy( hFeat );
     401                 : 
     402             201 :     return eErr;
     403                 : }
     404                 : 
     405                 : /************************************************************************/
     406                 : /*                          GPMaskImageData()                           */
     407                 : /*                                                                      */
     408                 : /*      Mask out image pixels to a special nodata value if the mask     */
     409                 : /*      band is zero.                                                   */
     410                 : /************************************************************************/
     411                 : 
     412                 : static CPLErr 
     413              28 : GPMaskImageData( GDALRasterBandH hMaskBand, GByte* pabyMaskLine, int iY, int nXSize, 
     414                 :                  GInt32 *panImageLine )
     415                 : 
     416                 : {
     417                 :     CPLErr eErr;
     418                 : 
     419                 :     eErr = GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1, 
     420              28 :                          pabyMaskLine, nXSize, 1, GDT_Byte, 0, 0 );
     421              28 :     if( eErr == CE_None )
     422                 :     {
     423                 :         int i;
     424             168 :         for( i = 0; i < nXSize; i++ )
     425                 :         {
     426             140 :             if( pabyMaskLine[i] == 0 )
     427              32 :                 panImageLine[i] = GP_NODATA_MARKER;
     428                 :         }
     429                 :     }
     430                 : 
     431              28 :     return eErr;
     432                 : }
     433                 : #endif // OGR_ENABLED
     434                 : 
     435                 : /************************************************************************/
     436                 : /*                           GDALPolygonize()                           */
     437                 : /************************************************************************/
     438                 : 
     439                 : /**
     440                 :  * Create polygon coverage from raster data.
     441                 :  *
     442                 :  * This function creates vector polygons for all connected regions of pixels in
     443                 :  * the raster sharing a common pixel value.  Optionally each polygon may be
     444                 :  * labelled with the pixel value in an attribute.  Optionally a mask band
     445                 :  * can be provided to determine which pixels are eligible for processing.
     446                 :  *
     447                 :  * Note that currently the source pixel band values are read into a
     448                 :  * signed 32bit integer buffer (Int32), so floating point or complex 
     449                 :  * bands will be implicitly truncated before processing. If you want to use a
     450                 :  * version using 32bit float buffers, see GDALFPolygonize() at fpolygonize.cpp.
     451                 :  *
     452                 :  * Polygon features will be created on the output layer, with polygon 
     453                 :  * geometries representing the polygons.  The polygon geometries will be
     454                 :  * in the georeferenced coordinate system of the image (based on the
     455                 :  * geotransform of the source dataset).  It is acceptable for the output
     456                 :  * layer to already have features.  Note that GDALPolygonize() does not
     457                 :  * set the coordinate system on the output layer.  Application code should
     458                 :  * do this when the layer is created, presumably matching the raster 
     459                 :  * coordinate system. 
     460                 :  *
     461                 :  * The algorithm used attempts to minimize memory use so that very large
     462                 :  * rasters can be processed.  However, if the raster has many polygons 
     463                 :  * or very large/complex polygons, the memory use for holding polygon 
     464                 :  * enumerations and active polygon geometries may grow to be quite large. 
     465                 :  *
     466                 :  * The algorithm will generally produce very dense polygon geometries, with
     467                 :  * edges that follow exactly on pixel boundaries for all non-interior pixels.
     468                 :  * For non-thematic raster data (such as satellite images) the result will
     469                 :  * essentially be one small polygon per pixel, and memory and output layer
     470                 :  * sizes will be substantial.  The algorithm is primarily intended for 
     471                 :  * relatively simple thematic imagery, masks, and classification results. 
     472                 :  * 
     473                 :  * @param hSrcBand the source raster band to be processed.
     474                 :  * @param hMaskBand an optional mask band.  All pixels in the mask band with a 
     475                 :  * value other than zero will be considered suitable for collection as 
     476                 :  * polygons.  
     477                 :  * @param hOutLayer the vector feature layer to which the polygons should
     478                 :  * be written. 
     479                 :  * @param iPixValField the attribute field index indicating the feature
     480                 :  * attribute into which the pixel value of the polygon should be written.
     481                 :  * @param papszOptions a name/value list of additional options
     482                 :  * <dl>
     483                 :  * <dt>"8CONNECTED":</dt> May be set to "8" to use 8 connectedness.
     484                 :  * Otherwise 4 connectedness will be applied to the algorithm
     485                 :  * </dl>
     486                 :  * @param pfnProgress callback for reporting algorithm progress matching the
     487                 :  * GDALProgressFunc() semantics.  May be NULL.
     488                 :  * @param pProgressArg callback argument passed to pfnProgress.
     489                 :  * 
     490                 :  * @return CE_None on success or CE_Failure on a failure.
     491                 :  */
     492                 : 
     493                 : CPLErr CPL_STDCALL
     494               6 : GDALPolygonize( GDALRasterBandH hSrcBand, 
     495                 :                 GDALRasterBandH hMaskBand,
     496                 :                 OGRLayerH hOutLayer, int iPixValField, 
     497                 :                 char **papszOptions,
     498                 :                 GDALProgressFunc pfnProgress, 
     499                 :                 void * pProgressArg )
     500                 : 
     501                 : {
     502                 : #ifndef OGR_ENABLED
     503                 :     CPLError(CE_Failure, CPLE_NotSupported, "GDALPolygonize() unimplemented in a non OGR build");
     504                 :     return CE_Failure;
     505                 : #else
     506               6 :     VALIDATE_POINTER1( hSrcBand, "GDALPolygonize", CE_Failure );
     507               6 :     VALIDATE_POINTER1( hOutLayer, "GDALPolygonize", CE_Failure );
     508                 : 
     509               6 :     if( pfnProgress == NULL )
     510               5 :         pfnProgress = GDALDummyProgress;
     511                 : 
     512               6 :     int nConnectedness = CSLFetchNameValue( papszOptions, "8CONNECTED" ) ? 8 : 4;
     513                 : 
     514                 : /* -------------------------------------------------------------------- */
     515                 : /*      Confirm our output layer will support feature creation.         */
     516                 : /* -------------------------------------------------------------------- */
     517               6 :     if( !OGR_L_TestCapability( hOutLayer, OLCSequentialWrite ) )
     518                 :     {
     519                 :         CPLError( CE_Failure, CPLE_AppDefined,
     520                 :                   "Output feature layer does not appear to support creation\n"
     521               0 :                   "of features in GDALPolygonize()." );
     522               0 :         return CE_Failure;
     523                 :     }
     524                 : 
     525                 : /* -------------------------------------------------------------------- */
     526                 : /*      Allocate working buffers.                                       */
     527                 : /* -------------------------------------------------------------------- */
     528               6 :     CPLErr eErr = CE_None;
     529               6 :     int nXSize = GDALGetRasterBandXSize( hSrcBand );
     530               6 :     int nYSize = GDALGetRasterBandYSize( hSrcBand );
     531               6 :     GInt32 *panLastLineVal = (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);
     532               6 :     GInt32 *panThisLineVal = (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);
     533               6 :     GInt32 *panLastLineId =  (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);
     534               6 :     GInt32 *panThisLineId =  (GInt32 *) VSIMalloc2(sizeof(GInt32),nXSize + 2);
     535               6 :     GByte *pabyMaskLine = (hMaskBand != NULL) ? (GByte *) VSIMalloc(nXSize) : NULL;
     536               6 :     if (panLastLineVal == NULL || panThisLineVal == NULL ||
     537                 :         panLastLineId == NULL || panThisLineId == NULL ||
     538                 :         (hMaskBand != NULL && pabyMaskLine == NULL))
     539                 :     {
     540                 :         CPLError(CE_Failure, CPLE_OutOfMemory,
     541               0 :                  "Could not allocate enough memory for temporary buffers");
     542               0 :         CPLFree( panThisLineId );
     543               0 :         CPLFree( panLastLineId );
     544               0 :         CPLFree( panThisLineVal );
     545               0 :         CPLFree( panLastLineVal );
     546               0 :         CPLFree( pabyMaskLine );
     547               0 :         return CE_Failure;
     548                 :     }
     549                 : 
     550                 : /* -------------------------------------------------------------------- */
     551                 : /*      Get the geotransform, if there is one, so we can convert the    */
     552                 : /*      vectors into georeferenced coordinates.                         */
     553                 : /* -------------------------------------------------------------------- */
     554               6 :     GDALDatasetH hSrcDS = GDALGetBandDataset( hSrcBand );
     555               6 :     double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
     556                 : 
     557               6 :     if( hSrcDS )
     558               6 :         GDALGetGeoTransform( hSrcDS, adfGeoTransform );
     559                 : 
     560                 : /* -------------------------------------------------------------------- */
     561                 : /*      The first pass over the raster is only used to build up the     */
     562                 : /*      polygon id map so we will know in advance what polygons are     */
     563                 : /*      what on the second pass.                                        */
     564                 : /* -------------------------------------------------------------------- */
     565                 :     int iY;
     566               6 :     GDALRasterPolygonEnumerator oFirstEnum(nConnectedness);
     567                 : 
     568              83 :     for( iY = 0; eErr == CE_None && iY < nYSize; iY++ )
     569                 :     {
     570                 :         eErr = GDALRasterIO( 
     571                 :             hSrcBand,
     572                 :             GF_Read, 0, iY, nXSize, 1, 
     573              77 :             panThisLineVal, nXSize, 1, GDT_Int32, 0, 0 );
     574                 :         
     575              77 :         if( eErr == CE_None && hMaskBand != NULL )
     576              14 :             eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize, panThisLineVal );
     577                 : 
     578              77 :         if( iY == 0 )
     579                 :             oFirstEnum.ProcessLine( 
     580               6 :                 NULL, panThisLineVal, NULL, panThisLineId, nXSize );
     581                 :         else
     582                 :             oFirstEnum.ProcessLine(
     583                 :                 panLastLineVal, panThisLineVal, 
     584                 :                 panLastLineId,  panThisLineId, 
     585              71 :                 nXSize );
     586                 : 
     587                 :         // swap lines
     588              77 :         GInt32 *panTmp = panLastLineVal;
     589              77 :         panLastLineVal = panThisLineVal;
     590              77 :         panThisLineVal = panTmp;
     591                 : 
     592              77 :         panTmp = panThisLineId;
     593              77 :         panThisLineId = panLastLineId;
     594              77 :         panLastLineId = panTmp;
     595                 : 
     596                 : /* -------------------------------------------------------------------- */
     597                 : /*      Report progress, and support interrupts.                        */
     598                 : /* -------------------------------------------------------------------- */
     599              77 :         if( eErr == CE_None 
     600                 :             && !pfnProgress( 0.10 * ((iY+1) / (double) nYSize), 
     601                 :                              "", pProgressArg ) )
     602                 :         {
     603               0 :             CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
     604               0 :             eErr = CE_Failure;
     605                 :         }
     606                 :     }
     607                 : 
     608                 : /* -------------------------------------------------------------------- */
     609                 : /*      Make a pass through the maps, ensuring every polygon id         */
     610                 : /*      points to the final id it should use, not an intermediate       */
     611                 : /*      value.                                                          */
     612                 : /* -------------------------------------------------------------------- */
     613               6 :     oFirstEnum.CompleteMerges();
     614                 : 
     615                 : /* -------------------------------------------------------------------- */
     616                 : /*      Initialize ids to -1 to serve as a nodata value for the         */
     617                 : /*      previous line, and past the beginning and end of the            */
     618                 : /*      scanlines.                                                      */
     619                 : /* -------------------------------------------------------------------- */
     620                 :     int iX;
     621                 : 
     622               6 :     panThisLineId[0] = -1;
     623               6 :     panThisLineId[nXSize+1] = -1;
     624                 : 
     625              79 :     for( iX = 0; iX < nXSize+2; iX++ )
     626              73 :         panLastLineId[iX] = -1;
     627                 : 
     628                 : /* -------------------------------------------------------------------- */
     629                 : /*      We will use a new enumerator for the second pass primariliy     */
     630                 : /*      so we can preserve the first pass map.                          */
     631                 : /* -------------------------------------------------------------------- */
     632               6 :     GDALRasterPolygonEnumerator oSecondEnum(nConnectedness);
     633                 :     RPolygon **papoPoly = (RPolygon **) 
     634               6 :         CPLCalloc(sizeof(RPolygon*),oFirstEnum.nNextPolygonId);
     635                 : 
     636                 : /* ==================================================================== */
     637                 : /*      Second pass during which we will actually collect polygon       */
     638                 : /*      edges as geometries.                                            */
     639                 : /* ==================================================================== */
     640              89 :     for( iY = 0; eErr == CE_None && iY < nYSize+1; iY++ )
     641                 :     {
     642                 : /* -------------------------------------------------------------------- */
     643                 : /*      Read the image data.                                            */
     644                 : /* -------------------------------------------------------------------- */
     645              83 :         if( iY < nYSize )
     646                 :         {
     647                 :             eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1, 
     648              77 :                                  panThisLineVal, nXSize, 1, GDT_Int32, 0, 0 );
     649                 : 
     650              77 :             if( eErr == CE_None && hMaskBand != NULL )
     651              14 :                 eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize, panThisLineVal );
     652                 :         }
     653                 : 
     654              83 :         if( eErr != CE_None )
     655               0 :             continue;
     656                 : 
     657                 : /* -------------------------------------------------------------------- */
     658                 : /*      Determine what polygon the various pixels belong to (redoing    */
     659                 : /*      the same thing done in the first pass above).                   */
     660                 : /* -------------------------------------------------------------------- */
     661              83 :         if( iY == nYSize )
     662                 :         {
     663              79 :             for( iX = 0; iX < nXSize+2; iX++ )
     664              73 :                 panThisLineId[iX] = -1;
     665                 :         }
     666              77 :         else if( iY == 0 )
     667                 :             oSecondEnum.ProcessLine( 
     668               6 :                 NULL, panThisLineVal, NULL, panThisLineId+1, nXSize );
     669                 :         else
     670                 :             oSecondEnum.ProcessLine(
     671                 :                 panLastLineVal, panThisLineVal, 
     672                 :                 panLastLineId+1,  panThisLineId+1, 
     673              71 :                 nXSize );
     674                 : 
     675                 : /* -------------------------------------------------------------------- */
     676                 : /*      Add polygon edges to our polygon list for the pixel             */
     677                 : /*      boundaries within and above this line.                          */
     678                 : /* -------------------------------------------------------------------- */
     679            1914 :         for( iX = 0; iX < nXSize+1; iX++ )
     680                 :         {
     681                 :             AddEdges( panThisLineId, panLastLineId, 
     682                 :                       oFirstEnum.panPolyIdMap, oFirstEnum.panPolyValue,
     683            1831 :                       papoPoly, iX, iY );
     684                 :         }
     685                 : 
     686                 : /* -------------------------------------------------------------------- */
     687                 : /*      Periodically we scan out polygons and write out those that      */
     688                 : /*      haven't been added to on the last line as we can be sure        */
     689                 : /*      they are complete.                                              */
     690                 : /* -------------------------------------------------------------------- */
     691              83 :         if( iY % 8 == 7 )
     692                 :         {
     693             491 :             for( iX = 0; 
     694                 :                  eErr == CE_None && iX < oSecondEnum.nNextPolygonId; 
     695                 :                  iX++ )
     696                 :             {
     697             481 :                 if( papoPoly[iX] && papoPoly[iX]->nLastLineUpdated < iY-1 )
     698                 :                 {
     699             194 :                     if( hMaskBand == NULL
     700              24 :                         || papoPoly[iX]->nPolyValue != GP_NODATA_MARKER )
     701                 :                     {
     702                 :                         eErr = 
     703                 :                             EmitPolygonToLayer( hOutLayer, iPixValField, 
     704             162 :                                                 papoPoly[iX], adfGeoTransform );
     705                 :                     }
     706             170 :                     delete papoPoly[iX];
     707             170 :                     papoPoly[iX] = NULL;
     708                 :                 }
     709                 :             }
     710                 :         }
     711                 : 
     712                 : /* -------------------------------------------------------------------- */
     713                 : /*      Swap pixel value, and polygon id lines to be ready for the      */
     714                 : /*      next line.                                                      */
     715                 : /* -------------------------------------------------------------------- */
     716              83 :         GInt32 *panTmp = panLastLineVal;
     717              83 :         panLastLineVal = panThisLineVal;
     718              83 :         panThisLineVal = panTmp;
     719                 : 
     720              83 :         panTmp = panThisLineId;
     721              83 :         panThisLineId = panLastLineId;
     722              83 :         panLastLineId = panTmp;
     723                 : 
     724                 : /* -------------------------------------------------------------------- */
     725                 : /*      Report progress, and support interrupts.                        */
     726                 : /* -------------------------------------------------------------------- */
     727              83 :         if( eErr == CE_None 
     728                 :             && !pfnProgress( 0.10 + 0.90 * ((iY+1) / (double) nYSize), 
     729                 :                              "", pProgressArg ) )
     730                 :         {
     731               0 :             CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
     732               0 :             eErr = CE_Failure;
     733                 :         }
     734                 :     }
     735                 : 
     736                 : /* -------------------------------------------------------------------- */
     737                 : /*      Make a cleanup pass for all unflushed polygons.                 */
     738                 : /* -------------------------------------------------------------------- */
     739             236 :     for( iX = 0; eErr == CE_None && iX < oSecondEnum.nNextPolygonId; iX++ )
     740                 :     {
     741             230 :         if( papoPoly[iX] )
     742                 :         {
     743              49 :             if( hMaskBand == NULL
     744              10 :                 || papoPoly[iX]->nPolyValue != GP_NODATA_MARKER )
     745                 :             {
     746                 :                 eErr = 
     747                 :                     EmitPolygonToLayer( hOutLayer, iPixValField, 
     748              39 :                                         papoPoly[iX], adfGeoTransform );
     749                 :             }
     750              39 :             delete papoPoly[iX];
     751              39 :             papoPoly[iX] = NULL;
     752                 :         }
     753                 :     }
     754                 : 
     755                 : /* -------------------------------------------------------------------- */
     756                 : /*      Cleanup                                                         */
     757                 : /* -------------------------------------------------------------------- */
     758               6 :     CPLFree( panThisLineId );
     759               6 :     CPLFree( panLastLineId );
     760               6 :     CPLFree( panThisLineVal );
     761               6 :     CPLFree( panLastLineVal );
     762               6 :     CPLFree( pabyMaskLine );
     763               6 :     CPLFree( papoPoly );
     764                 : 
     765               6 :     return eErr;
     766                 : #endif // OGR_ENABLED
     767                 : }
     768                 : 

Generated by: LCOV version 1.7