← Back to branch summary

~ubuntu-branches/debian/sid/gdal/sid

« back to all changes in this revision

Viewing changes to alg/fpolygonize.cpp

  • Committer: Package Import Robot
  • Author(s): Francesco Paolo Lovergine
  • Date: 2012-05-07 15:04:42 UTC
  • mfrom: (5.5.16 experimental)
  • Revision ID: package-import@ubuntu.com-20120507150442-2eks97loeh6rq005
Tags: 1.9.0-1
* Ready for sid, starting transition.
* All symfiles updated to latest builds.
* Added dh_numpy call in debian/rules to depend on numpy ABI.
* Policy bumped to 3.9.3, no changes required.

Show diffs side-by-side

added added

removed removed

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