1 /******************************************************************************
2 * Project: GDAL
3 * Purpose: Raster to Polygon Converter
4 * Author: Frank Warmerdam, warmerdam@pobox.com
5 *
6 ******************************************************************************
7 * Copyright (c) 2008, Frank Warmerdam
8 * Copyright (c) 2009-2020, Even Rouault <even dot rouault at spatialys.com>
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 "cpl_port.h"
30 #include "gdal_alg.h"
31
32 #include <stddef.h>
33 #include <stdio.h>
34 #include <cstdlib>
35 #include <string.h>
36
37 #include <algorithm>
38 #include <map>
39 #include <memory>
40 #include <utility>
41 #include <vector>
42
43 #include "gdal_alg_priv.h"
44 #include "gdal.h"
45 #include "ogr_api.h"
46 #include "ogr_core.h"
47 #include "cpl_conv.h"
48 #include "cpl_error.h"
49 #include "cpl_progress.h"
50 #include "cpl_string.h"
51 #include "cpl_vsi.h"
52
53 CPL_CVSID("$Id: polygonize.cpp 0d6994bbbb6aa91af042f65fbb9d24f2155484e7 2021-07-02 23:46:01 +0800 XuXiang $")
54
55 /************************************************************************/
56 /* ==================================================================== */
57 /* RPolygon */
58 /* */
59 /* This is a helper class to hold polygons while they are being */
60 /* formed in memory, and to provide services to coalesce a much */
61 /* of edge sections into complete rings. */
62 /* ==================================================================== */
63 /************************************************************************/
64
65 class RPolygon {
66 public:
67 double dfPolyValue = 0.0;
68 int nLastLineUpdated = -1;
69
70 struct XY
71 {
72 int x;
73 int y;
74
operator <RPolygon::XY75 bool operator< (const XY& other) const
76 {
77 if( x < other.x )
78 return true;
79 if( x > other.x )
80 return false;
81 return y < other.y;
82 }
83
operator ==RPolygon::XY84 bool operator== (const XY& other) const { return x == other.x && y == other.y; }
85 };
86
87 typedef int StringId;
88 typedef std::map< XY, std::pair<StringId, StringId>> MapExtremity;
89
90 std::map< StringId, std::vector<XY> > oMapStrings{};
91 MapExtremity oMapStartStrings{};
92 MapExtremity oMapEndStrings{};
93 StringId iNextStringId = 0;
94
95 static
96 StringId findExtremityNot(const MapExtremity& oMap,
97 const XY& xy,
98 StringId excludedId);
99
100 static void removeExtremity(MapExtremity& oMap,
101 const XY& xy,
102 StringId id);
103
104 static void insertExtremity(MapExtremity& oMap,
105 const XY& xy,
106 StringId id);
107
RPolygon(double dfValue)108 explicit RPolygon( double dfValue ): dfPolyValue(dfValue) {}
109
110 void AddSegment( int x1, int y1, int x2, int y2, int direction );
111 void Dump() const;
112 void Coalesce();
113 void Merge( StringId iBaseString, StringId iSrcString, int iDirection );
114 };
115
116 /************************************************************************/
117 /* Dump() */
118 /************************************************************************/
Dump() const119 void RPolygon::Dump() const
120 {
121 /*ok*/printf( "RPolygon: Value=%g, LastLineUpdated=%d\n",
122 dfPolyValue, nLastLineUpdated );
123
124 for( const auto& oStringIter: oMapStrings )
125 {
126 /*ok*/printf( " String " CPL_FRMT_GIB ":\n", static_cast<GIntBig>(oStringIter.first) );
127 for( const auto& xy: oStringIter.second )
128 {
129 /*ok*/printf( " (%d,%d)\n", xy.x, xy.y );
130 }
131 }
132 }
133
134 /************************************************************************/
135 /* findExtremityNot() */
136 /************************************************************************/
137
findExtremityNot(const MapExtremity & oMap,const XY & xy,StringId excludedId)138 RPolygon::StringId RPolygon::findExtremityNot(const MapExtremity& oMap,
139 const XY& xy,
140 StringId excludedId)
141 {
142 auto oIter = oMap.find(xy);
143 if( oIter == oMap.end() )
144 return -1;
145 if( oIter->second.first != excludedId )
146 return oIter->second.first;
147 if( oIter->second.second != excludedId )
148 return oIter->second.second;
149 return -1;
150 }
151
152 /************************************************************************/
153 /* Coalesce() */
154 /************************************************************************/
155
Coalesce()156 void RPolygon::Coalesce()
157
158 {
159 /* -------------------------------------------------------------------- */
160 /* Iterate over loops starting from the first, trying to merge */
161 /* other segments into them. */
162 /* -------------------------------------------------------------------- */
163 for( auto& oStringIter: oMapStrings )
164 {
165 const auto thisId = oStringIter.first;
166 auto& oString = oStringIter.second;
167
168 /* -------------------------------------------------------------------- */
169 /* Keep trying to merge others strings into our target "base" */
170 /* string while there are matches. */
171 /* -------------------------------------------------------------------- */
172 while(true)
173 {
174 auto nOtherId = findExtremityNot(oMapStartStrings, oString.back(), thisId);
175 if( nOtherId != -1 )
176 {
177 Merge( thisId, nOtherId, 1 );
178 continue;
179 }
180 else
181 {
182 nOtherId = findExtremityNot(oMapEndStrings, oString.back(), thisId);
183 if( nOtherId != -1 )
184 {
185 Merge( thisId, nOtherId, -1 );
186 continue;
187 }
188 }
189 break;
190 }
191
192 // At this point our loop *should* be closed!
193 CPLAssert( oString.front() == oString.back() );
194 }
195 }
196
197 /************************************************************************/
198 /* removeExtremity() */
199 /************************************************************************/
200
removeExtremity(MapExtremity & oMap,const XY & xy,StringId id)201 void RPolygon::removeExtremity(MapExtremity& oMap,
202 const XY& xy,
203 StringId id)
204 {
205 auto oIter = oMap.find(xy);
206 CPLAssert( oIter != oMap.end() );
207 if( oIter->second.first == id )
208 {
209 oIter->second.first = oIter->second.second;
210 oIter->second.second = -1;
211 if( oIter->second.first < 0 )
212 oMap.erase(oIter);
213 }
214 else if( oIter->second.second == id )
215 {
216 oIter->second.second = -1;
217 CPLAssert( oIter->second.first >= 0 );
218 }
219 else
220 {
221 CPLAssert(false);
222 }
223 }
224
225 /************************************************************************/
226 /* insertExtremity() */
227 /************************************************************************/
228
insertExtremity(MapExtremity & oMap,const XY & xy,StringId id)229 void RPolygon::insertExtremity(MapExtremity& oMap,
230 const XY& xy,
231 StringId id)
232 {
233 auto oIter = oMap.find(xy);
234 if( oIter != oMap.end() )
235 {
236 CPLAssert( oIter->second.second == -1 );
237 oIter->second.second = id;
238 }
239 else
240 {
241 oMap[xy] = std::pair<StringId, StringId>(id, -1);
242 }
243 }
244
245 /************************************************************************/
246 /* Merge() */
247 /************************************************************************/
248
Merge(StringId iBaseString,StringId iSrcString,int iDirection)249 void RPolygon::Merge( StringId iBaseString, StringId iSrcString, int iDirection )
250
251 {
252 auto &anBase = oMapStrings.find(iBaseString)->second;
253 auto anStringIter = oMapStrings.find(iSrcString);
254 auto& anString = anStringIter->second;
255 int iStart = 1;
256 int iEnd = -1;
257
258 if( iDirection == 1 )
259 {
260 iEnd = static_cast<int>(anString.size());
261 }
262 else
263 {
264 iStart = static_cast<int>(anString.size()) - 2;
265 }
266
267 removeExtremity(oMapEndStrings, anBase.back(), iBaseString);
268
269 anBase.reserve(anBase.size() + anString.size() - 1);
270 for( int i = iStart; i != iEnd; i += iDirection )
271 {
272 anBase.push_back( anString[i] );
273 }
274
275 removeExtremity(oMapStartStrings, anString.front(), iSrcString);
276 removeExtremity(oMapEndStrings, anString.back(), iSrcString);
277 oMapStrings.erase(anStringIter);
278 insertExtremity(oMapEndStrings, anBase.back(), iBaseString);
279 }
280
281 /************************************************************************/
282 /* AddSegment() */
283 /************************************************************************/
284
AddSegment(int x1,int y1,int x2,int y2,int direction)285 void RPolygon::AddSegment( int x1, int y1, int x2, int y2, int direction )
286
287 {
288 nLastLineUpdated = std::max(y1, y2);
289
290 /* -------------------------------------------------------------------- */
291 /* Is there an existing string ending with this? */
292 /* -------------------------------------------------------------------- */
293
294 XY xy1 = { x1, y1 };
295 XY xy2 = { x2, y2 };
296
297 StringId iExistingString = findExtremityNot(oMapEndStrings, xy1, -1);
298 if( iExistingString >= 0 )
299 {
300 StringId iExistingString2 = findExtremityNot(oMapEndStrings, xy1, iExistingString);
301 if ( iExistingString2 >= 0 )
302 {
303 // If there are two strings that ending with this segment
304 // choose the string which is in the same pixel with this segment
305 auto& anString2 = oMapStrings[iExistingString2];
306 const size_t nSSize2 = anString2.size();
307
308 if ( x1 == x2 )
309 {
310 // vertical segment input, choose the horizontal string
311 if ( anString2[nSSize2 - 2].y == anString2[nSSize2 - 1].y )
312 {
313 iExistingString = iExistingString2;
314 }
315 }
316 else
317 {
318 // horizontal segment input, choose the vertical string
319 if ( anString2[nSSize2 - 2].x == anString2[nSSize2 - 1].x )
320 {
321 iExistingString = iExistingString2;
322 }
323 }
324 }
325 auto& anString = oMapStrings[iExistingString];
326 const size_t nSSize = anString.size();
327
328 // We are going to add a segment, but should we just extend
329 // an existing segment already going in the right direction?
330
331 const int nLastLen =
332 std::max(std::abs(anString[nSSize - 2].x - anString[nSSize - 1].x),
333 std::abs(anString[nSSize - 2].y - anString[nSSize - 1].y));
334
335 removeExtremity(oMapEndStrings, anString.back(), iExistingString);
336
337 if( (anString[nSSize - 2].x - anString[nSSize - 1].x
338 == (anString[nSSize - 1].x - xy2.x) * nLastLen)
339 && (anString[nSSize - 2].y - anString[nSSize - 1].y
340 == (anString[nSSize - 1].y - xy2.y) * nLastLen) )
341 {
342 anString[nSSize - 1] = xy2;
343 }
344 else
345 {
346 anString.push_back( xy2 );
347 }
348 insertExtremity(oMapEndStrings, anString.back(), iExistingString);
349 }
350 else
351 {
352 oMapStrings[iNextStringId] = std::vector<XY>{ xy1, xy2 };
353 insertExtremity(oMapStartStrings, xy1, iNextStringId);
354 insertExtremity(oMapEndStrings, xy2, iNextStringId);
355
356 iExistingString = iNextStringId;
357 iNextStringId ++;
358 }
359
360 // merge rings if possible
361 if ( direction == 1 ) {
362 StringId iExistingString2 = findExtremityNot(oMapEndStrings, xy2, iExistingString);
363 if ( iExistingString2 >= 0 )
364 {
365 this->Merge(iExistingString, iExistingString2, -1);
366 }
367 }
368 }
369
370 /************************************************************************/
371 /* ==================================================================== */
372 /* End of RPolygon */
373 /* ==================================================================== */
374 /************************************************************************/
375
376 /************************************************************************/
377 /* AddEdges() */
378 /* */
379 /* Examine one pixel and compare to its neighbour above */
380 /* (previous) and right. If they are different polygon ids */
381 /* then add the pixel edge to this polygon and the one on the */
382 /* other side of the edge. */
383 /************************************************************************/
384
385 template<class DataType>
AddEdges(GInt32 * panThisLineId,GInt32 * panLastLineId,GInt32 * panPolyIdMap,DataType * panPolyValue,RPolygon ** papoPoly,int iX,int iY)386 static void AddEdges( GInt32 *panThisLineId, GInt32 *panLastLineId,
387 GInt32 *panPolyIdMap, DataType *panPolyValue,
388 RPolygon **papoPoly, int iX, int iY )
389
390 {
391 // TODO(schwehr): Simplify these three vars.
392 int nThisId = panThisLineId[iX];
393 if( nThisId != -1 )
394 nThisId = panPolyIdMap[nThisId];
395 int nRightId = panThisLineId[iX+1];
396 if( nRightId != -1 )
397 nRightId = panPolyIdMap[nRightId];
398 int nPreviousId = panLastLineId[iX];
399 if( nPreviousId != -1 )
400 nPreviousId = panPolyIdMap[nPreviousId];
401
402 const int iXReal = iX - 1;
403
404 if( nThisId != nPreviousId )
405 {
406 if( nThisId != -1 )
407 {
408 if( papoPoly[nThisId] == nullptr )
409 papoPoly[nThisId] = new RPolygon( panPolyValue[nThisId] );
410
411 papoPoly[nThisId]->AddSegment( iXReal, iY, iXReal+1, iY, 1 );
412 }
413 if( nPreviousId != -1 )
414 {
415 if( papoPoly[nPreviousId] == nullptr )
416 papoPoly[nPreviousId] = new RPolygon(panPolyValue[nPreviousId]);
417
418 papoPoly[nPreviousId]->AddSegment( iXReal, iY, iXReal+1, iY, 0 );
419 }
420 }
421
422 if( nThisId != nRightId )
423 {
424 if( nThisId != -1 )
425 {
426 if( papoPoly[nThisId] == nullptr )
427 papoPoly[nThisId] = new RPolygon(panPolyValue[nThisId]);
428
429 papoPoly[nThisId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1, 1 );
430 }
431
432 if( nRightId != -1 )
433 {
434 if( papoPoly[nRightId] == nullptr )
435 papoPoly[nRightId] = new RPolygon(panPolyValue[nRightId]);
436
437 papoPoly[nRightId]->AddSegment( iXReal+1, iY, iXReal+1, iY+1, 0 );
438 }
439 }
440 }
441
442 /************************************************************************/
443 /* EmitPolygonToLayer() */
444 /************************************************************************/
445
446 static CPLErr
EmitPolygonToLayer(OGRLayerH hOutLayer,int iPixValField,RPolygon * poRPoly,double * padfGeoTransform)447 EmitPolygonToLayer( OGRLayerH hOutLayer, int iPixValField,
448 RPolygon *poRPoly, double *padfGeoTransform )
449
450 {
451 /* -------------------------------------------------------------------- */
452 /* Turn bits of lines into coherent rings. */
453 /* -------------------------------------------------------------------- */
454 poRPoly->Coalesce();
455
456 /* -------------------------------------------------------------------- */
457 /* Create the polygon geometry. */
458 /* -------------------------------------------------------------------- */
459 OGRGeometryH hPolygon = OGR_G_CreateGeometry( wkbPolygon );
460
461 for( const auto& oIter: poRPoly->oMapStrings )
462 {
463 const auto &anString = oIter.second;
464 OGRGeometryH hRing = OGR_G_CreateGeometry( wkbLinearRing );
465
466 // We go last to first to ensure the linestring is allocated to
467 // the proper size on the first try.
468 for( int iVert = static_cast<int>(anString.size()) - 1;
469 iVert >= 0;
470 iVert-- )
471 {
472 const int nPixelX = anString[iVert].x;
473 const int nPixelY = anString[iVert].y;
474
475 const double dfX =
476 padfGeoTransform[0]
477 + nPixelX * padfGeoTransform[1]
478 + nPixelY * padfGeoTransform[2];
479 const double dfY =
480 padfGeoTransform[3]
481 + nPixelX * padfGeoTransform[4]
482 + nPixelY * padfGeoTransform[5];
483
484 OGR_G_SetPoint_2D( hRing, iVert, dfX, dfY );
485 }
486
487 OGR_G_AddGeometryDirectly( hPolygon, hRing );
488 }
489
490 /* -------------------------------------------------------------------- */
491 /* Create the feature object. */
492 /* -------------------------------------------------------------------- */
493 OGRFeatureH hFeat = OGR_F_Create( OGR_L_GetLayerDefn( hOutLayer ) );
494
495 OGR_F_SetGeometryDirectly( hFeat, hPolygon );
496
497 if( iPixValField >= 0 )
498 OGR_F_SetFieldDouble( hFeat, iPixValField, poRPoly->dfPolyValue );
499
500 /* -------------------------------------------------------------------- */
501 /* Write the to the layer. */
502 /* -------------------------------------------------------------------- */
503 CPLErr eErr = CE_None;
504
505 if( OGR_L_CreateFeature( hOutLayer, hFeat ) != OGRERR_NONE )
506 eErr = CE_Failure;
507
508 OGR_F_Destroy( hFeat );
509
510 return eErr;
511 }
512
513 /************************************************************************/
514 /* GPMaskImageData() */
515 /* */
516 /* Mask out image pixels to a special nodata value if the mask */
517 /* band is zero. */
518 /************************************************************************/
519
520 template<class DataType>
521 static CPLErr
GPMaskImageData(GDALRasterBandH hMaskBand,GByte * pabyMaskLine,int iY,int nXSize,DataType * panImageLine)522 GPMaskImageData( GDALRasterBandH hMaskBand, GByte* pabyMaskLine,
523 int iY, int nXSize,
524 DataType *panImageLine )
525
526 {
527 const CPLErr eErr =
528 GDALRasterIO( hMaskBand, GF_Read, 0, iY, nXSize, 1,
529 pabyMaskLine, nXSize, 1, GDT_Byte, 0, 0 );
530 if( eErr != CE_None )
531 return eErr;
532
533 for( int i = 0; i < nXSize; i++ )
534 {
535 if( pabyMaskLine[i] == 0 )
536 panImageLine[i] = GP_NODATA_MARKER;
537 }
538
539 return CE_None;
540 }
541
542 /************************************************************************/
543 /* GDALPolygonizeT() */
544 /************************************************************************/
545
546 template<class DataType, class EqualityTest>
547 static CPLErr
GDALPolygonizeT(GDALRasterBandH hSrcBand,GDALRasterBandH hMaskBand,OGRLayerH hOutLayer,int iPixValField,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg,GDALDataType eDT)548 GDALPolygonizeT( GDALRasterBandH hSrcBand,
549 GDALRasterBandH hMaskBand,
550 OGRLayerH hOutLayer, int iPixValField,
551 char **papszOptions,
552 GDALProgressFunc pfnProgress,
553 void * pProgressArg,
554 GDALDataType eDT)
555
556 {
557 VALIDATE_POINTER1( hSrcBand, "GDALPolygonize", CE_Failure );
558 VALIDATE_POINTER1( hOutLayer, "GDALPolygonize", CE_Failure );
559
560 if( pfnProgress == nullptr )
561 pfnProgress = GDALDummyProgress;
562
563 const int nConnectedness =
564 CSLFetchNameValue( papszOptions, "8CONNECTED" ) ? 8 : 4;
565
566 /* -------------------------------------------------------------------- */
567 /* Confirm our output layer will support feature creation. */
568 /* -------------------------------------------------------------------- */
569 if( !OGR_L_TestCapability( hOutLayer, OLCSequentialWrite ) )
570 {
571 CPLError(CE_Failure, CPLE_AppDefined,
572 "Output feature layer does not appear to support creation "
573 "of features in GDALPolygonize().");
574 return CE_Failure;
575 }
576
577 /* -------------------------------------------------------------------- */
578 /* Allocate working buffers. */
579 /* -------------------------------------------------------------------- */
580 const int nXSize = GDALGetRasterBandXSize( hSrcBand );
581 const int nYSize = GDALGetRasterBandYSize( hSrcBand );
582
583 DataType *panLastLineVal = static_cast<DataType *>(
584 VSI_MALLOC2_VERBOSE(sizeof(DataType), nXSize + 2));
585 DataType *panThisLineVal = static_cast<DataType *>(
586 VSI_MALLOC2_VERBOSE(sizeof(DataType), nXSize + 2));
587 GInt32 *panLastLineId = static_cast<GInt32 *>(
588 VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize + 2));
589 GInt32 *panThisLineId = static_cast<GInt32 *>(
590 VSI_MALLOC2_VERBOSE(sizeof(GInt32), nXSize + 2));
591
592 GByte *pabyMaskLine =
593 hMaskBand != nullptr
594 ? static_cast<GByte *>(VSI_MALLOC_VERBOSE(nXSize))
595 : nullptr;
596
597 if( panLastLineVal == nullptr || panThisLineVal == nullptr ||
598 panLastLineId == nullptr || panThisLineId == nullptr ||
599 (hMaskBand != nullptr && pabyMaskLine == nullptr) )
600 {
601 CPLFree( panThisLineId );
602 CPLFree( panLastLineId );
603 CPLFree( panThisLineVal );
604 CPLFree( panLastLineVal );
605 CPLFree( pabyMaskLine );
606 return CE_Failure;
607 }
608
609 /* -------------------------------------------------------------------- */
610 /* Get the geotransform, if there is one, so we can convert the */
611 /* vectors into georeferenced coordinates. */
612 /* -------------------------------------------------------------------- */
613 double adfGeoTransform[6] = { 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 };
614 bool bGotGeoTransform = false;
615 const char* pszDatasetForGeoRef = CSLFetchNameValue(papszOptions,
616 "DATASET_FOR_GEOREF");
617 if( pszDatasetForGeoRef )
618 {
619 GDALDatasetH hSrcDS = GDALOpen(pszDatasetForGeoRef, GA_ReadOnly);
620 if( hSrcDS )
621 {
622 bGotGeoTransform = GDALGetGeoTransform( hSrcDS, adfGeoTransform ) == CE_None;
623 GDALClose(hSrcDS);
624 }
625 }
626 else
627 {
628 GDALDatasetH hSrcDS = GDALGetBandDataset( hSrcBand );
629 if( hSrcDS )
630 bGotGeoTransform = GDALGetGeoTransform( hSrcDS, adfGeoTransform ) == CE_None;
631 }
632 if( !bGotGeoTransform )
633 {
634 adfGeoTransform[0] = 0;
635 adfGeoTransform[1] = 1;
636 adfGeoTransform[2] = 0;
637 adfGeoTransform[3] = 0;
638 adfGeoTransform[4] = 0;
639 adfGeoTransform[5] = 1;
640 }
641
642 /* -------------------------------------------------------------------- */
643 /* The first pass over the raster is only used to build up the */
644 /* polygon id map so we will know in advance what polygons are */
645 /* what on the second pass. */
646 /* -------------------------------------------------------------------- */
647 GDALRasterPolygonEnumeratorT<DataType,
648 EqualityTest> oFirstEnum(nConnectedness);
649
650 CPLErr eErr = CE_None;
651
652 for( int iY = 0; eErr == CE_None && iY < nYSize; iY++ )
653 {
654 eErr = GDALRasterIO(
655 hSrcBand,
656 GF_Read, 0, iY, nXSize, 1,
657 panThisLineVal, nXSize, 1, eDT, 0, 0 );
658
659 if( eErr == CE_None && hMaskBand != nullptr )
660 eErr = GPMaskImageData(hMaskBand, pabyMaskLine, iY, nXSize,
661 panThisLineVal);
662
663 if( iY == 0 )
664 oFirstEnum.ProcessLine(
665 nullptr, panThisLineVal, nullptr, panThisLineId, nXSize );
666 else
667 oFirstEnum.ProcessLine(
668 panLastLineVal, panThisLineVal,
669 panLastLineId, panThisLineId,
670 nXSize );
671
672 // Swap lines.
673 std::swap(panLastLineVal, panThisLineVal);
674 std::swap(panLastLineId, panThisLineId);
675
676 /* -------------------------------------------------------------------- */
677 /* Report progress, and support interrupts. */
678 /* -------------------------------------------------------------------- */
679 if( eErr == CE_None
680 && !pfnProgress( 0.10 * ((iY+1) / static_cast<double>(nYSize)),
681 "", pProgressArg ) )
682 {
683 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
684 eErr = CE_Failure;
685 }
686 }
687
688 /* -------------------------------------------------------------------- */
689 /* Make a pass through the maps, ensuring every polygon id */
690 /* points to the final id it should use, not an intermediate */
691 /* value. */
692 /* -------------------------------------------------------------------- */
693 oFirstEnum.CompleteMerges();
694
695 /* -------------------------------------------------------------------- */
696 /* Initialize ids to -1 to serve as a nodata value for the */
697 /* previous line, and past the beginning and end of the */
698 /* scanlines. */
699 /* -------------------------------------------------------------------- */
700 panThisLineId[0] = -1;
701 panThisLineId[nXSize+1] = -1;
702
703 for( int iX = 0; iX < nXSize+2; iX++ )
704 panLastLineId[iX] = -1;
705
706 /* -------------------------------------------------------------------- */
707 /* We will use a new enumerator for the second pass primarily */
708 /* so we can preserve the first pass map. */
709 /* -------------------------------------------------------------------- */
710 GDALRasterPolygonEnumeratorT<DataType,
711 EqualityTest> oSecondEnum(nConnectedness);
712 RPolygon **papoPoly = static_cast<RPolygon **>(
713 CPLCalloc(sizeof(RPolygon*), oFirstEnum.nNextPolygonId));
714
715 /* ==================================================================== */
716 /* Second pass during which we will actually collect polygon */
717 /* edges as geometries. */
718 /* ==================================================================== */
719 for( int iY = 0; eErr == CE_None && iY < nYSize+1; iY++ )
720 {
721 /* -------------------------------------------------------------------- */
722 /* Read the image data. */
723 /* -------------------------------------------------------------------- */
724 if( iY < nYSize )
725 {
726 eErr = GDALRasterIO( hSrcBand, GF_Read, 0, iY, nXSize, 1,
727 panThisLineVal, nXSize, 1, eDT, 0, 0 );
728
729 if( eErr == CE_None && hMaskBand != nullptr )
730 eErr = GPMaskImageData( hMaskBand, pabyMaskLine, iY, nXSize,
731 panThisLineVal );
732 }
733
734 if( eErr != CE_None )
735 continue;
736
737 /* -------------------------------------------------------------------- */
738 /* Determine what polygon the various pixels belong to (redoing */
739 /* the same thing done in the first pass above). */
740 /* -------------------------------------------------------------------- */
741 if( iY == nYSize )
742 {
743 for( int iX = 0; iX < nXSize+2; iX++ )
744 panThisLineId[iX] = -1;
745 }
746 else if( iY == 0 )
747 {
748 oSecondEnum.ProcessLine(
749 nullptr, panThisLineVal, nullptr, panThisLineId+1, nXSize );
750 }
751 else
752 {
753 oSecondEnum.ProcessLine(
754 panLastLineVal, panThisLineVal,
755 panLastLineId+1, panThisLineId+1,
756 nXSize );
757 }
758
759 /* -------------------------------------------------------------------- */
760 /* Add polygon edges to our polygon list for the pixel */
761 /* boundaries within and above this line. */
762 /* -------------------------------------------------------------------- */
763 for( int iX = 0; iX < nXSize+1; iX++ )
764 {
765 AddEdges( panThisLineId, panLastLineId,
766 oFirstEnum.panPolyIdMap, oFirstEnum.panPolyValue,
767 papoPoly, iX, iY );
768 }
769
770 /* -------------------------------------------------------------------- */
771 /* Periodically we scan out polygons and write out those that */
772 /* haven't been added to on the last line as we can be sure */
773 /* they are complete. */
774 /* -------------------------------------------------------------------- */
775 if( iY % 8 == 7 )
776 {
777 for( int iX = 0;
778 eErr == CE_None && iX < oSecondEnum.nNextPolygonId;
779 iX++ )
780 {
781 if( papoPoly[iX] && papoPoly[iX]->nLastLineUpdated < iY-1 )
782 {
783 eErr =
784 EmitPolygonToLayer( hOutLayer, iPixValField,
785 papoPoly[iX], adfGeoTransform );
786
787 delete papoPoly[iX];
788 papoPoly[iX] = nullptr;
789 }
790 }
791 }
792
793 /* -------------------------------------------------------------------- */
794 /* Swap pixel value, and polygon id lines to be ready for the */
795 /* next line. */
796 /* -------------------------------------------------------------------- */
797 std::swap(panLastLineVal, panThisLineVal);
798 std::swap(panLastLineId, panThisLineId);
799
800 /* -------------------------------------------------------------------- */
801 /* Report progress, and support interrupts. */
802 /* -------------------------------------------------------------------- */
803 if( eErr == CE_None
804 && !pfnProgress( 0.10 + 0.90 * ((iY + 1) /
805 static_cast<double>(nYSize)),
806 "", pProgressArg ) )
807 {
808 CPLError( CE_Failure, CPLE_UserInterrupt, "User terminated" );
809 eErr = CE_Failure;
810 }
811 }
812
813 /* -------------------------------------------------------------------- */
814 /* Make a cleanup pass for all unflushed polygons. */
815 /* -------------------------------------------------------------------- */
816 for( int iX = 0; eErr == CE_None && iX < oSecondEnum.nNextPolygonId; iX++ )
817 {
818 if( papoPoly[iX] )
819 {
820 eErr = EmitPolygonToLayer( hOutLayer, iPixValField,
821 papoPoly[iX], adfGeoTransform );
822
823 delete papoPoly[iX];
824 papoPoly[iX] = nullptr;
825 }
826 }
827
828 /* -------------------------------------------------------------------- */
829 /* Cleanup */
830 /* -------------------------------------------------------------------- */
831 CPLFree( panThisLineId );
832 CPLFree( panLastLineId );
833 CPLFree( panThisLineVal );
834 CPLFree( panLastLineVal );
835 CPLFree( pabyMaskLine );
836 CPLFree( papoPoly );
837
838 return eErr;
839 }
840
841 /******************************************************************************/
842 /* GDALFloatEquals() */
843 /* Code from: */
844 /* http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm */
845 /******************************************************************************/
GDALFloatEquals(float A,float B)846 GBool GDALFloatEquals( float A, float B )
847 {
848 // This function will allow maxUlps-1 floats between A and B.
849 const int maxUlps = MAX_ULPS;
850
851 // Make sure maxUlps is non-negative and small enough that the default NAN
852 // won't compare as equal to anything.
853 #if MAX_ULPS <= 0 || MAX_ULPS >= 4 * 1024 * 1024
854 #error "Invalid MAX_ULPS"
855 #endif
856
857 // This assignation could violate strict aliasing. It causes a warning with
858 // gcc -O2. Use of memcpy preferred. Credits for Even Rouault. Further info
859 // at http://trac.osgeo.org/gdal/ticket/4005#comment:6
860 int aInt = 0;
861 memcpy(&aInt, &A, 4);
862
863 // Make aInt lexicographically ordered as a twos-complement int.
864 if( aInt < 0 )
865 aInt = INT_MIN - aInt;
866
867 // Make bInt lexicographically ordered as a twos-complement int.
868 int bInt = 0;
869 memcpy(&bInt, &B, 4);
870
871 if( bInt < 0 )
872 bInt = INT_MIN - bInt;
873 #ifdef COMPAT_WITH_ICC_CONVERSION_CHECK
874 const int intDiff =
875 abs(static_cast<int>(static_cast<GUIntBig>(
876 static_cast<GIntBig>(aInt) - static_cast<GIntBig>(bInt))
877 & 0xFFFFFFFFU));
878 #else
879 // To make -ftrapv happy we compute the diff on larger type and
880 // cast down later.
881 const int intDiff = abs(static_cast<int>(
882 static_cast<GIntBig>(aInt) - static_cast<GIntBig>(bInt)));
883 #endif
884 if( intDiff <= maxUlps )
885 return true;
886 return false;
887 }
888
889 /************************************************************************/
890 /* GDALPolygonize() */
891 /************************************************************************/
892
893 /**
894 * Create polygon coverage from raster data.
895 *
896 * This function creates vector polygons for all connected regions of pixels in
897 * the raster sharing a common pixel value. Optionally each polygon may be
898 * labeled with the pixel value in an attribute. Optionally a mask band
899 * can be provided to determine which pixels are eligible for processing.
900 *
901 * Note that currently the source pixel band values are read into a
902 * signed 32bit integer buffer (Int32), so floating point or complex
903 * bands will be implicitly truncated before processing. If you want to use a
904 * version using 32bit float buffers, see GDALFPolygonize().
905 *
906 * Polygon features will be created on the output layer, with polygon
907 * geometries representing the polygons. The polygon geometries will be
908 * in the georeferenced coordinate system of the image (based on the
909 * geotransform of the source dataset). It is acceptable for the output
910 * layer to already have features. Note that GDALPolygonize() does not
911 * set the coordinate system on the output layer. Application code should
912 * do this when the layer is created, presumably matching the raster
913 * coordinate system.
914 *
915 * The algorithm used attempts to minimize memory use so that very large
916 * rasters can be processed. However, if the raster has many polygons
917 * or very large/complex polygons, the memory use for holding polygon
918 * enumerations and active polygon geometries may grow to be quite large.
919 *
920 * The algorithm will generally produce very dense polygon geometries, with
921 * edges that follow exactly on pixel boundaries for all non-interior pixels.
922 * For non-thematic raster data (such as satellite images) the result will
923 * essentially be one small polygon per pixel, and memory and output layer
924 * sizes will be substantial. The algorithm is primarily intended for
925 * relatively simple thematic imagery, masks, and classification results.
926 *
927 * @param hSrcBand the source raster band to be processed.
928 * @param hMaskBand an optional mask band. All pixels in the mask band with a
929 * value other than zero will be considered suitable for collection as
930 * polygons.
931 * @param hOutLayer the vector feature layer to which the polygons should
932 * be written.
933 * @param iPixValField the attribute field index indicating the feature
934 * attribute into which the pixel value of the polygon should be written. Or
935 * -1 to indicate that the pixel value must not be written.
936 * @param papszOptions a name/value list of additional options
937 * <ul>
938 * <li>8CONNECTED=8: May be set to "8" to use 8 connectedness.
939 * Otherwise 4 connectedness will be applied to the algorithm</li>
940 * </ul>
941 * @param pfnProgress callback for reporting algorithm progress matching the
942 * GDALProgressFunc() semantics. May be NULL.
943 * @param pProgressArg callback argument passed to pfnProgress.
944 *
945 * @return CE_None on success or CE_Failure on a failure.
946 */
947
948 CPLErr CPL_STDCALL
GDALPolygonize(GDALRasterBandH hSrcBand,GDALRasterBandH hMaskBand,OGRLayerH hOutLayer,int iPixValField,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)949 GDALPolygonize( GDALRasterBandH hSrcBand,
950 GDALRasterBandH hMaskBand,
951 OGRLayerH hOutLayer, int iPixValField,
952 char **papszOptions,
953 GDALProgressFunc pfnProgress,
954 void * pProgressArg )
955
956 {
957 return GDALPolygonizeT<GInt32, IntEqualityTest>(hSrcBand,
958 hMaskBand,
959 hOutLayer,
960 iPixValField,
961 papszOptions,
962 pfnProgress,
963 pProgressArg,
964 GDT_Int32);
965 }
966
967 /************************************************************************/
968 /* GDALFPolygonize() */
969 /************************************************************************/
970
971 /**
972 * Create polygon coverage from raster data.
973 *
974 * This function creates vector polygons for all connected regions of pixels in
975 * the raster sharing a common pixel value. Optionally each polygon may be
976 * labeled with the pixel value in an attribute. Optionally a mask band
977 * can be provided to determine which pixels are eligible for processing.
978 *
979 * The source pixel band values are read into a 32bit float buffer. If you want
980 * to use a (probably faster) version using signed 32bit integer buffer, see
981 * GDALPolygonize().
982 *
983 * Polygon features will be created on the output layer, with polygon
984 * geometries representing the polygons. The polygon geometries will be
985 * in the georeferenced coordinate system of the image (based on the
986 * geotransform of the source dataset). It is acceptable for the output
987 * layer to already have features. Note that GDALFPolygonize() does not
988 * set the coordinate system on the output layer. Application code should
989 * do this when the layer is created, presumably matching the raster
990 * coordinate system.
991 *
992 * The algorithm used attempts to minimize memory use so that very large
993 * rasters can be processed. However, if the raster has many polygons
994 * or very large/complex polygons, the memory use for holding polygon
995 * enumerations and active polygon geometries may grow to be quite large.
996 *
997 * The algorithm will generally produce very dense polygon geometries, with
998 * edges that follow exactly on pixel boundaries for all non-interior pixels.
999 * For non-thematic raster data (such as satellite images) the result will
1000 * essentially be one small polygon per pixel, and memory and output layer
1001 * sizes will be substantial. The algorithm is primarily intended for
1002 * relatively simple thematic imagery, masks, and classification results.
1003 *
1004 * @param hSrcBand the source raster band to be processed.
1005 * @param hMaskBand an optional mask band. All pixels in the mask band with a
1006 * value other than zero will be considered suitable for collection as
1007 * polygons.
1008 * @param hOutLayer the vector feature layer to which the polygons should
1009 * be written.
1010 * @param iPixValField the attribute field index indicating the feature
1011 * attribute into which the pixel value of the polygon should be written. Or
1012 * -1 to indicate that the pixel value must not be written.
1013 * @param papszOptions a name/value list of additional options
1014 * <ul>
1015 * <li>8CONNECTED=8: May be set to "8" to use 8 connectedness.
1016 * Otherwise 4 connectedness will be applied to the algorithm</li>
1017 * </ul>
1018 * @param pfnProgress callback for reporting algorithm progress matching the
1019 * GDALProgressFunc() semantics. May be NULL.
1020 * @param pProgressArg callback argument passed to pfnProgress.
1021 *
1022 * @return CE_None on success or CE_Failure on a failure.
1023 *
1024 * @since GDAL 1.9.0
1025 */
1026
1027 CPLErr CPL_STDCALL
GDALFPolygonize(GDALRasterBandH hSrcBand,GDALRasterBandH hMaskBand,OGRLayerH hOutLayer,int iPixValField,char ** papszOptions,GDALProgressFunc pfnProgress,void * pProgressArg)1028 GDALFPolygonize( GDALRasterBandH hSrcBand,
1029 GDALRasterBandH hMaskBand,
1030 OGRLayerH hOutLayer, int iPixValField,
1031 char **papszOptions,
1032 GDALProgressFunc pfnProgress,
1033 void * pProgressArg )
1034
1035 {
1036 return GDALPolygonizeT<float, FloatEqualityTest>(hSrcBand,
1037 hMaskBand,
1038 hOutLayer,
1039 iPixValField,
1040 papszOptions,
1041 pfnProgress,
1042 pProgressArg,
1043 GDT_Float32);
1044 }
1045