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