GDAL
polygonize_polygonizer.cpp
1/******************************************************************************
2 * Project: GDAL
3 * Purpose: Implements The Two-Arm Chains EdgeTracing Algorithm
4 * Author: kikitte.lee
5 *
6 ******************************************************************************
7 * Copyright (c) 2023, kikitte.lee <kikitte.lee@gmail.com>
8 *
9 * SPDX-License-Identifier: MIT
10 ****************************************************************************/
11
14#include "polygonize_polygonizer.h"
15
16#include <algorithm>
17
18namespace gdal
19{
20namespace polygonizer
21{
22IndexedArc RPolygon::newArc(bool bFollowRighthand)
23{
24 const std::size_t iArcIndex = oArcs.size();
25 const auto &oNewArc =
26 oArcs.emplace_back(static_cast<unsigned>(iArcIndex), bFollowRighthand);
27 return IndexedArc{oNewArc.poArc.get(), iArcIndex};
28}
29
30void RPolygon::setArcConnection(const IndexedArc &oArc,
31 const IndexedArc &oNextArc)
32{
33 oArcs[oArc.iIndex].nConnection = static_cast<unsigned>(oNextArc.iIndex);
34}
35
36void RPolygon::updateBottomRightPos(IndexType iRow, IndexType iCol)
37{
38 iBottomRightRow = iRow;
39 iBottomRightCol = iCol;
40}
41
45static void ProcessArmConnections(TwoArm *poCurrent, TwoArm *poAbove,
46 TwoArm *poLeft)
47{
48 poCurrent->poPolyInside->updateBottomRightPos(poCurrent->iRow,
49 poCurrent->iCol);
50 poCurrent->bSolidVertical = poCurrent->poPolyInside != poLeft->poPolyInside;
51 poCurrent->bSolidHorizontal =
52 poCurrent->poPolyInside != poAbove->poPolyInside;
53 poCurrent->poPolyAbove = poAbove->poPolyInside;
54 poCurrent->poPolyLeft = poLeft->poPolyInside;
55
56 constexpr int BIT_CUR_HORIZ = 0;
57 constexpr int BIT_CUR_VERT = 1;
58 constexpr int BIT_LEFT = 2;
59 constexpr int BIT_ABOVE = 3;
60
61 const int nArmConnectionType =
62 (static_cast<int>(poAbove->bSolidVertical) << BIT_ABOVE) |
63 (static_cast<int>(poLeft->bSolidHorizontal) << BIT_LEFT) |
64 (static_cast<int>(poCurrent->bSolidVertical) << BIT_CUR_VERT) |
65 (static_cast<int>(poCurrent->bSolidHorizontal) << BIT_CUR_HORIZ);
66
67 constexpr int VIRTUAL = 0;
68 constexpr int SOLID = 1;
69
70 constexpr int ABOVE_VIRTUAL = VIRTUAL << BIT_ABOVE;
71 constexpr int ABOVE_SOLID = SOLID << BIT_ABOVE;
72
73 constexpr int LEFT_VIRTUAL = VIRTUAL << BIT_LEFT;
74 constexpr int LEFT_SOLID = SOLID << BIT_LEFT;
75
76 constexpr int CUR_VERT_VIRTUAL = VIRTUAL << BIT_CUR_VERT;
77 constexpr int CUR_VERT_SOLID = SOLID << BIT_CUR_VERT;
78
79 constexpr int CUR_HORIZ_VIRTUAL = VIRTUAL << BIT_CUR_HORIZ;
80 constexpr int CUR_HORIZ_SOLID = SOLID << BIT_CUR_HORIZ;
81
108 switch (nArmConnectionType)
109 {
110 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
111 CUR_HORIZ_VIRTUAL: // 0
112 // nothing to do
113 break;
114
115 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
116 CUR_HORIZ_SOLID: // 3
117 // add inner arcs
118 poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
119 poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
120 poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
121 poCurrent->oArcVerInner);
122 poCurrent->oArcVerInner.poArc->emplace_back(
123 Point{poCurrent->iRow, poCurrent->iCol});
124
125 // add outer arcs
126 poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
127 poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
128 poAbove->poPolyInside->setArcConnection(poCurrent->oArcVerOuter,
129 poCurrent->oArcHorOuter);
130 poCurrent->oArcHorOuter.poArc->push_back(
131 Point{poCurrent->iRow, poCurrent->iCol});
132
133 break;
134 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
135 CUR_HORIZ_SOLID: // 5
136 // pass arcs
137 poCurrent->oArcHorInner = poLeft->oArcHorInner;
138 poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
139
140 break;
141 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
142 CUR_HORIZ_VIRTUAL: // 6
143 // pass arcs
144 poCurrent->oArcVerInner = poLeft->oArcHorOuter;
145 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
146 poCurrent->oArcVerInner.poArc->push_back(
147 Point{poCurrent->iRow, poCurrent->iCol});
148 poCurrent->oArcVerOuter.poArc->push_back(
149 Point{poCurrent->iRow, poCurrent->iCol});
150
151 break;
152 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
153 CUR_HORIZ_SOLID: // 7
154 // pass arcs
155 poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
156 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
157 poLeft->oArcHorInner.poArc->push_back(
158 Point{poCurrent->iRow, poCurrent->iCol});
159
160 // add inner arcs
161 poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
162 poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
163 poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
164 poCurrent->oArcVerInner);
165 poCurrent->oArcVerInner.poArc->push_back(
166 Point{poCurrent->iRow, poCurrent->iCol});
167
168 break;
169 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
170 CUR_HORIZ_SOLID: // 9
171 // pass arcs
172 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
173 poCurrent->oArcHorInner = poAbove->oArcVerOuter;
174 poCurrent->oArcHorOuter.poArc->push_back(
175 Point{poCurrent->iRow, poCurrent->iCol});
176 poCurrent->oArcHorInner.poArc->push_back(
177 Point{poCurrent->iRow, poCurrent->iCol});
178
179 break;
180 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
181 CUR_HORIZ_VIRTUAL: // 10
182 // pass arcs
183 poCurrent->oArcVerInner = poAbove->oArcVerInner;
184 poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
185
186 break;
187 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
188 CUR_HORIZ_SOLID: // 11
189 // pass arcs
190 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
191 poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
192 poCurrent->oArcHorOuter.poArc->push_back(
193 Point{poCurrent->iRow, poCurrent->iCol});
194 // add inner arcs
195 poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
196 poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
197 poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
198 poCurrent->oArcVerInner);
199 poCurrent->oArcVerInner.poArc->push_back(
200 Point{poCurrent->iRow, poCurrent->iCol});
201
202 break;
203 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
204 CUR_HORIZ_VIRTUAL: // 12
205 // close arcs
206 poLeft->oArcHorOuter.poArc->push_back(
207 Point{poCurrent->iRow, poCurrent->iCol});
208 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
209 poAbove->oArcVerOuter);
210 // close arcs
211 poAbove->oArcVerInner.poArc->push_back(
212 Point{poCurrent->iRow, poCurrent->iCol});
213 poCurrent->poPolyInside->setArcConnection(poAbove->oArcVerInner,
214 poLeft->oArcHorInner);
215
216 break;
217 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
218 CUR_HORIZ_SOLID: // 13
219 // close arcs
220 poLeft->oArcHorOuter.poArc->push_back(
221 Point{poCurrent->iRow, poCurrent->iCol});
222 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
223 poAbove->oArcVerOuter);
224 // pass arcs
225 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
226 poCurrent->oArcHorInner = poLeft->oArcHorInner;
227 poCurrent->oArcHorOuter.poArc->push_back(
228 Point{poCurrent->iRow, poCurrent->iCol});
229
230 break;
231 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID |
232 CUR_HORIZ_VIRTUAL: // 14
233 // close arcs
234 poLeft->oArcHorOuter.poArc->push_back(
235 Point{poCurrent->iRow, poCurrent->iCol});
236 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
237 poAbove->oArcVerOuter);
238 // pass arcs
239 poCurrent->oArcVerInner = poAbove->oArcVerInner;
240 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
241 poCurrent->oArcVerOuter.poArc->push_back(
242 Point{poCurrent->iRow, poCurrent->iCol});
243
244 break;
245 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID | CUR_HORIZ_SOLID: // 15
246 // Tow pixels of the main diagonal belong to the same polygon
247 if (poAbove->poPolyLeft == poCurrent->poPolyInside)
248 {
249 // pass arcs
250 poCurrent->oArcVerInner = poLeft->oArcHorOuter;
251 poCurrent->oArcHorInner = poAbove->oArcVerOuter;
252 poCurrent->oArcVerInner.poArc->push_back(
253 Point{poCurrent->iRow, poCurrent->iCol});
254 poCurrent->oArcHorInner.poArc->push_back(
255 Point{poCurrent->iRow, poCurrent->iCol});
256 }
257 else
258 {
259 // close arcs
260 poLeft->oArcHorOuter.poArc->push_back(
261 Point{poCurrent->iRow, poCurrent->iCol});
262 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
263 poAbove->oArcVerOuter);
264 // add inner arcs
265 poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
266 poCurrent->oArcHorInner =
267 poCurrent->poPolyInside->newArc(false);
268 poCurrent->poPolyInside->setArcConnection(
269 poCurrent->oArcHorInner, poCurrent->oArcVerInner);
270 poCurrent->oArcVerInner.poArc->push_back(
271 Point{poCurrent->iRow, poCurrent->iCol});
272 }
273
274 // Tow pixels of the secondary diagonal belong to the same polygon
275 if (poAbove->poPolyInside == poLeft->poPolyInside)
276 {
277 // close arcs
278 poAbove->poPolyInside->setArcConnection(poAbove->oArcVerInner,
279 poLeft->oArcHorInner);
280 poAbove->oArcVerInner.poArc->push_back(
281 Point{poCurrent->iRow, poCurrent->iCol});
282 // add outer arcs
283 poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
284 poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
285 poCurrent->oArcHorOuter.poArc->push_back(
286 Point{poCurrent->iRow, poCurrent->iCol});
287 poAbove->poPolyInside->setArcConnection(
288 poCurrent->oArcVerOuter, poCurrent->oArcHorOuter);
289 }
290 else
291 {
292 // pass arcs
293 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
294 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
295 poCurrent->oArcHorOuter.poArc->push_back(
296 Point{poCurrent->iRow, poCurrent->iCol});
297 poCurrent->oArcVerOuter.poArc->push_back(
298 Point{poCurrent->iRow, poCurrent->iCol});
299 }
300
301 break;
302
303 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
304 CUR_HORIZ_SOLID: // 1
305 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
306 CUR_HORIZ_VIRTUAL: // 2
307 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
308 CUR_HORIZ_VIRTUAL: // 4
309 default:
310 // Impossible case
311 CPLAssert(false);
312 break;
313 }
314}
315
316template <typename PolyIdType, typename DataType>
317Polygonizer<PolyIdType, DataType>::Polygonizer(
318 PolyIdType nInvalidPolyId, PolygonReceiver<DataType> *poPolygonReceiver)
319 : nInvalidPolyId_(nInvalidPolyId), poPolygonReceiver_(poPolygonReceiver)
320{
321 poTheOuterPolygon_ = createPolygon(THE_OUTER_POLYGON_ID);
322}
323
324template <typename PolyIdType, typename DataType>
325Polygonizer<PolyIdType, DataType>::~Polygonizer()
326{
327 // cppcheck-suppress constVariableReference
328 for (auto &pair : oPolygonMap_)
329 {
330 delete pair.second;
331 }
332}
333
334template <typename PolyIdType, typename DataType>
335RPolygon *Polygonizer<PolyIdType, DataType>::getPolygon(PolyIdType nPolygonId)
336{
337 const auto oIter = oPolygonMap_.find(nPolygonId);
338 if (oIter == oPolygonMap_.end())
339 {
340 return createPolygon(nPolygonId);
341 }
342 else
343 {
344 return oIter->second;
345 }
346}
347
348template <typename PolyIdType, typename DataType>
349RPolygon *
350Polygonizer<PolyIdType, DataType>::createPolygon(PolyIdType nPolygonId)
351{
352 auto polygon = new RPolygon();
353 oPolygonMap_[nPolygonId] = polygon;
354 return polygon;
355}
356
357template <typename PolyIdType, typename DataType>
358void Polygonizer<PolyIdType, DataType>::destroyPolygon(PolyIdType nPolygonId)
359{
360 const auto oIter = oPolygonMap_.find(nPolygonId);
361 CPLAssert(oIter != oPolygonMap_.end());
362 delete oIter->second;
363 oPolygonMap_.erase(oIter);
364}
365
366template <typename PolyIdType, typename DataType>
367bool Polygonizer<PolyIdType, DataType>::processLine(
368 const PolyIdType *panThisLineId, const DataType *panLastLineVal,
369 TwoArm *poThisLineArm, TwoArm *poLastLineArm, const IndexType nCurrentRow,
370 const IndexType nCols)
371{
372 TwoArm *poCurrent, *poAbove, *poLeft;
373
374 try
375 {
376 poCurrent = poThisLineArm + 1;
377 poCurrent->iRow = nCurrentRow;
378 poCurrent->iCol = 0;
379 poCurrent->poPolyInside = getPolygon(panThisLineId[0]);
380 poAbove = poLastLineArm + 1;
381 poLeft = poThisLineArm;
382 poLeft->poPolyInside = poTheOuterPolygon_;
383 ProcessArmConnections(poCurrent, poAbove, poLeft);
384 for (IndexType col = 1; col < nCols; ++col)
385 {
386 IndexType iArmIndex = col + 1;
387 poCurrent = poThisLineArm + iArmIndex;
388 poCurrent->iRow = nCurrentRow;
389 poCurrent->iCol = col;
390 poCurrent->poPolyInside = getPolygon(panThisLineId[col]);
391 poAbove = poLastLineArm + iArmIndex;
392 poLeft = poThisLineArm + iArmIndex - 1;
393 ProcessArmConnections(poCurrent, poAbove, poLeft);
394 }
395 poCurrent = poThisLineArm + nCols + 1;
396 poCurrent->iRow = nCurrentRow;
397 poCurrent->iCol = nCols;
398 poCurrent->poPolyInside = poTheOuterPolygon_;
399 poAbove = poLastLineArm + nCols + 1;
400 poAbove->poPolyInside = poTheOuterPolygon_;
401 poLeft = poThisLineArm + nCols;
402 ProcessArmConnections(poCurrent, poAbove, poLeft);
403
409 std::vector<PolygonMapEntry> oCompletedPolygons;
410 for (auto &entry : oPolygonMap_)
411 {
412 RPolygon *poPolygon = entry.second;
413
414 if (poPolygon->iBottomRightRow + 1 == nCurrentRow)
415 {
416 oCompletedPolygons.push_back(entry);
417 }
418 }
419 // cppcheck-suppress constVariableReference
420 for (auto &entry : oCompletedPolygons)
421 {
422 PolyIdType nPolyId = entry.first;
423 RPolygon *poPolygon = entry.second;
424
425 // emit valid polygon only
426 if (nPolyId != nInvalidPolyId_)
427 {
428 poPolygonReceiver_->receive(
429 poPolygon, panLastLineVal[poPolygon->iBottomRightCol]);
430 }
431
432 destroyPolygon(nPolyId);
433 }
434 return true;
435 }
436 catch (const std::bad_alloc &)
437 {
438 CPLError(CE_Failure, CPLE_OutOfMemory,
439 "Out of memory in Polygonizer::processLine");
440 return false;
441 }
442}
443
444template <typename DataType>
445OGRPolygonWriter<DataType>::OGRPolygonWriter(OGRLayerH hOutLayer,
446 int iPixValField,
447 double *padfGeoTransform)
448 : PolygonReceiver<DataType>(), poOutLayer_(OGRLayer::FromHandle(hOutLayer)),
449 iPixValField_(iPixValField), padfGeoTransform_(padfGeoTransform)
450{
451 poFeature_ = std::make_unique<OGRFeature>(poOutLayer_->GetLayerDefn());
452 poPolygon_ = new OGRPolygon();
453 poFeature_->SetGeometryDirectly(poPolygon_);
454}
455
456template <typename DataType>
457void OGRPolygonWriter<DataType>::receive(RPolygon *poPolygon,
458 DataType nPolygonCellValue)
459{
460 std::vector<bool> oAccessedArc(poPolygon->oArcs.size(), false);
461 double *padfGeoTransform = padfGeoTransform_;
462
463 OGRLinearRing *poFirstRing = poPolygon_->getExteriorRing();
464 if (poFirstRing && poPolygon_->getNumInteriorRings() == 0)
465 {
466 poFirstRing->empty();
467 }
468 else
469 {
470 poFirstRing = nullptr;
471 poPolygon_->empty();
472 }
473
474 auto AddRingToPolygon =
475 [this, &poPolygon, &oAccessedArc,
476 padfGeoTransform](std::size_t iFirstArcIndex, OGRLinearRing *poRing)
477 {
478 std::unique_ptr<OGRLinearRing> poNewRing;
479 if (!poRing)
480 {
481 poNewRing = std::make_unique<OGRLinearRing>();
482 poRing = poNewRing.get();
483 }
484
485 auto AddArcToRing =
486 [&poPolygon, poRing, padfGeoTransform](std::size_t iArcIndex)
487 {
488 const auto &oArc = poPolygon->oArcs[iArcIndex];
489 const bool bArcFollowRighthand = oArc.bFollowRighthand;
490 const int nArcPointCount = static_cast<int>(oArc.poArc->size());
491 int nDstPointIdx = poRing->getNumPoints();
492 poRing->setNumPoints(nDstPointIdx + nArcPointCount,
493 /* bZeroizeNewContent = */ false);
494 if (poRing->getNumPoints() < nDstPointIdx + nArcPointCount)
495 {
496 return false;
497 }
498 for (int i = 0; i < nArcPointCount; ++i)
499 {
500 const Point &oPixel =
501 (*oArc.poArc)[bArcFollowRighthand
502 ? i
503 : (nArcPointCount - i - 1)];
504
505 const double dfX = padfGeoTransform[0] +
506 oPixel[1] * padfGeoTransform[1] +
507 oPixel[0] * padfGeoTransform[2];
508 const double dfY = padfGeoTransform[3] +
509 oPixel[1] * padfGeoTransform[4] +
510 oPixel[0] * padfGeoTransform[5];
511
512 poRing->setPoint(nDstPointIdx, dfX, dfY);
513 ++nDstPointIdx;
514 }
515 return true;
516 };
517
518 if (!AddArcToRing(iFirstArcIndex))
519 {
520 return false;
521 }
522
523 std::size_t iArcIndex = iFirstArcIndex;
524 std::size_t iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
525 oAccessedArc[iArcIndex] = true;
526 while (iNextArcIndex != iFirstArcIndex)
527 {
528 if (!AddArcToRing(iNextArcIndex))
529 {
530 return false;
531 }
532 iArcIndex = iNextArcIndex;
533 iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
534 oAccessedArc[iArcIndex] = true;
535 }
536
537 // close ring manually
538 poRing->closeRings();
539
540 if (poNewRing)
541 poPolygon_->addRingDirectly(poNewRing.release());
542 return true;
543 };
544
545 for (size_t i = 0; i < oAccessedArc.size(); ++i)
546 {
547 if (!oAccessedArc[i])
548 {
549 if (!AddRingToPolygon(i, poFirstRing))
550 {
551 eErr_ = CE_Failure;
552 return;
553 }
554 poFirstRing = nullptr;
555 }
556 }
557
558 // Create the feature object
559 poFeature_->SetFID(OGRNullFID);
560 if (iPixValField_ >= 0)
561 poFeature_->SetField(iPixValField_,
562 static_cast<double>(nPolygonCellValue));
563
564 // Write the to the layer.
565 if (poOutLayer_->CreateFeature(poFeature_.get()) != OGRERR_NONE)
566 eErr_ = CE_Failure;
567
568 // Shouldn't happen for well behaved drivers, but better check...
569 else if (poFeature_->GetGeometryRef() != poPolygon_)
570 {
571 poPolygon_ = new OGRPolygon();
572 poFeature_->SetGeometryDirectly(poPolygon_);
573 }
574}
575
576} // namespace polygonizer
577} // namespace gdal
578
OGRLayer
This class represents a layer of simple features, with access methods.
Definition ogrsf_frmts.h:58
OGRLinearRing
Concrete representation of a closed ring.
Definition ogr_geometry.h:1850
OGRPolygon
Concrete class representing polygons.
Definition ogr_geometry.h:2660
OGRSimpleCurve::empty
virtual void empty() override
Clear geometry information.
Definition ogrlinestring.cpp:116
CPLAssert
#define CPLAssert(expr)
Assert on an expression.
Definition cpl_error.h:209
CPLE_OutOfMemory
#define CPLE_OutOfMemory
Out of memory error.
Definition cpl_error.h:86
OGRLayerH
void * OGRLayerH
Opaque type for a layer (OGRLayer)
Definition ogr_api.h:676
OGRNullFID
#define OGRNullFID
Special value for a unset FID.
Definition ogr_core.h:848
OGRERR_NONE
#define OGRERR_NONE
Success.
Definition ogr_core.h:373
Generated by doxygen 1.9.8