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