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 :
|
