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