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