14#include "polygonize_polygonizer.h"
22IndexedArc RPolygon::newArc(
bool bFollowRighthand)
24 const std::size_t iArcIndex = oArcs.size();
26 oArcs.emplace_back(
static_cast<unsigned>(iArcIndex), bFollowRighthand);
27 return IndexedArc{oNewArc.poArc.get(), iArcIndex};
30void RPolygon::setArcConnection(
const IndexedArc &oArc,
31 const IndexedArc &oNextArc)
33 oArcs[oArc.iIndex].nConnection =
static_cast<unsigned>(oNextArc.iIndex);
36void RPolygon::updateBottomRightPos(IndexType iRow, IndexType iCol)
38 iBottomRightRow = iRow;
39 iBottomRightCol = iCol;
45static void ProcessArmConnections(TwoArm *poCurrent, TwoArm *poAbove,
48 poCurrent->poPolyInside->updateBottomRightPos(poCurrent->iRow,
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;
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;
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);
67 constexpr int VIRTUAL = 0;
68 constexpr int SOLID = 1;
70 constexpr int ABOVE_VIRTUAL = VIRTUAL << BIT_ABOVE;
71 constexpr int ABOVE_SOLID = SOLID << BIT_ABOVE;
73 constexpr int LEFT_VIRTUAL = VIRTUAL << BIT_LEFT;
74 constexpr int LEFT_SOLID = SOLID << BIT_LEFT;
76 constexpr int CUR_VERT_VIRTUAL = VIRTUAL << BIT_CUR_VERT;
77 constexpr int CUR_VERT_SOLID = SOLID << BIT_CUR_VERT;
79 constexpr int CUR_HORIZ_VIRTUAL = VIRTUAL << BIT_CUR_HORIZ;
80 constexpr int CUR_HORIZ_SOLID = SOLID << BIT_CUR_HORIZ;
108 switch (nArmConnectionType)
110 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
115 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
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});
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});
134 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
137 poCurrent->oArcHorInner = poLeft->oArcHorInner;
138 poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
141 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
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});
152 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
155 poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
156 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
157 poLeft->oArcHorInner.poArc->push_back(
158 Point{poCurrent->iRow, poCurrent->iCol});
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});
169 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
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});
180 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
183 poCurrent->oArcVerInner = poAbove->oArcVerInner;
184 poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
187 case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
190 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
191 poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
192 poCurrent->oArcHorOuter.poArc->push_back(
193 Point{poCurrent->iRow, poCurrent->iCol});
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});
203 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
206 poLeft->oArcHorOuter.poArc->push_back(
207 Point{poCurrent->iRow, poCurrent->iCol});
208 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
209 poAbove->oArcVerOuter);
211 poAbove->oArcVerInner.poArc->push_back(
212 Point{poCurrent->iRow, poCurrent->iCol});
213 poCurrent->poPolyInside->setArcConnection(poAbove->oArcVerInner,
214 poLeft->oArcHorInner);
217 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
220 poLeft->oArcHorOuter.poArc->push_back(
221 Point{poCurrent->iRow, poCurrent->iCol});
222 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
223 poAbove->oArcVerOuter);
225 poCurrent->oArcHorOuter = poAbove->oArcVerInner;
226 poCurrent->oArcHorInner = poLeft->oArcHorInner;
227 poCurrent->oArcHorOuter.poArc->push_back(
228 Point{poCurrent->iRow, poCurrent->iCol});
231 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID |
234 poLeft->oArcHorOuter.poArc->push_back(
235 Point{poCurrent->iRow, poCurrent->iCol});
236 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
237 poAbove->oArcVerOuter);
239 poCurrent->oArcVerInner = poAbove->oArcVerInner;
240 poCurrent->oArcVerOuter = poLeft->oArcHorInner;
241 poCurrent->oArcVerOuter.poArc->push_back(
242 Point{poCurrent->iRow, poCurrent->iCol});
245 case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID | CUR_HORIZ_SOLID:
247 if (poAbove->poPolyLeft == poCurrent->poPolyInside)
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});
260 poLeft->oArcHorOuter.poArc->push_back(
261 Point{poCurrent->iRow, poCurrent->iCol});
262 poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
263 poAbove->oArcVerOuter);
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});
275 if (poAbove->poPolyInside == poLeft->poPolyInside)
278 poAbove->poPolyInside->setArcConnection(poAbove->oArcVerInner,
279 poLeft->oArcHorInner);
280 poAbove->oArcVerInner.poArc->push_back(
281 Point{poCurrent->iRow, poCurrent->iCol});
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);
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});
303 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
305 case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
307 case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
316template <
typename PolyIdType,
typename DataType>
317Polygonizer<PolyIdType, DataType>::Polygonizer(
318 PolyIdType nInvalidPolyId, PolygonReceiver<DataType> *poPolygonReceiver)
319 : nInvalidPolyId_(nInvalidPolyId), poPolygonReceiver_(poPolygonReceiver)
321 poTheOuterPolygon_ = createPolygon(THE_OUTER_POLYGON_ID);
324template <
typename PolyIdType,
typename DataType>
325Polygonizer<PolyIdType, DataType>::~Polygonizer()
328 for (
auto &pair : oPolygonMap_)
334template <
typename PolyIdType,
typename DataType>
335RPolygon *Polygonizer<PolyIdType, DataType>::getPolygon(PolyIdType nPolygonId)
337 const auto oIter = oPolygonMap_.find(nPolygonId);
338 if (oIter == oPolygonMap_.end())
340 return createPolygon(nPolygonId);
344 return oIter->second;
348template <
typename PolyIdType,
typename DataType>
350Polygonizer<PolyIdType, DataType>::createPolygon(PolyIdType nPolygonId)
352 auto polygon =
new RPolygon();
353 oPolygonMap_[nPolygonId] = polygon;
357template <
typename PolyIdType,
typename DataType>
358void Polygonizer<PolyIdType, DataType>::destroyPolygon(PolyIdType nPolygonId)
360 const auto oIter = oPolygonMap_.find(nPolygonId);
362 delete oIter->second;
363 oPolygonMap_.erase(oIter);
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)
372 TwoArm *poCurrent, *poAbove, *poLeft;
376 poCurrent = poThisLineArm + 1;
377 poCurrent->iRow = nCurrentRow;
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)
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);
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);
409 std::vector<PolygonMapEntry> oCompletedPolygons;
410 for (
auto &entry : oPolygonMap_)
412 RPolygon *poPolygon = entry.second;
414 if (poPolygon->iBottomRightRow + 1 == nCurrentRow)
416 oCompletedPolygons.push_back(entry);
420 for (
auto &entry : oCompletedPolygons)
422 PolyIdType nPolyId = entry.first;
423 RPolygon *poPolygon = entry.second;
426 if (nPolyId != nInvalidPolyId_)
428 poPolygonReceiver_->receive(
429 poPolygon, panLastLineVal[poPolygon->iBottomRightCol]);
432 destroyPolygon(nPolyId);
436 catch (
const std::bad_alloc &)
439 "Out of memory in Polygonizer::processLine");
444template <
typename DataType>
445OGRPolygonWriter<DataType>::OGRPolygonWriter(
OGRLayerH hOutLayer,
449 : PolygonReceiver<DataType>(), poOutLayer_(
OGRLayer::FromHandle(hOutLayer)),
450 poOutDS_(poOutLayer_->GetDataset()), iPixValField_(iPixValField), gt_(gt),
451 nCommitInterval_(nCommitInterval)
453 if (nCommitInterval_ != 0)
459 nCommitInterval_ = 0;
462 poFeature_ = std::make_unique<OGRFeature>(poOutLayer_->GetLayerDefn());
464 poFeature_->SetGeometryDirectly(poPolygon_);
467template <
typename DataType> OGRPolygonWriter<DataType>::~OGRPolygonWriter()
469 OGRPolygonWriter::Finalize();
472template <
typename DataType>
bool OGRPolygonWriter<DataType>::Finalize()
474 if (nFeaturesWrittenInTransaction_)
476 nFeaturesWrittenInTransaction_ = 0;
485template <
typename DataType>
486void OGRPolygonWriter<DataType>::receive(RPolygon *poPolygon,
487 DataType nPolygonCellValue)
489 std::vector<bool> oAccessedArc(poPolygon->oArcs.size(),
false);
492 if (poFirstRing && poPolygon_->getNumInteriorRings() == 0)
494 poFirstRing->
empty();
498 poFirstRing =
nullptr;
502 auto AddRingToPolygon =
503 [
this, &poPolygon, &oAccessedArc](std::size_t iFirstArcIndex,
506 std::unique_ptr<OGRLinearRing> poNewRing;
509 poNewRing = std::make_unique<OGRLinearRing>();
510 poRing = poNewRing.get();
513 auto AddArcToRing = [
this, &poPolygon, poRing](std::size_t iArcIndex)
515 const auto &oArc = poPolygon->oArcs[iArcIndex];
516 const bool bArcFollowRighthand = oArc.bFollowRighthand;
517 const int nArcPointCount =
static_cast<int>(oArc.poArc->size());
518 int nDstPointIdx = poRing->getNumPoints();
519 poRing->setNumPoints(nDstPointIdx + nArcPointCount,
521 if (poRing->getNumPoints() < nDstPointIdx + nArcPointCount)
525 for (
int i = 0; i < nArcPointCount; ++i)
527 const Point &oPixel =
528 (*oArc.poArc)[bArcFollowRighthand
530 : (nArcPointCount - i - 1)];
532 const auto oGeoreferenced = gt_.Apply(oPixel[1], oPixel[0]);
533 poRing->setPoint(nDstPointIdx, oGeoreferenced.first,
534 oGeoreferenced.second);
540 if (!AddArcToRing(iFirstArcIndex))
545 std::size_t iArcIndex = iFirstArcIndex;
546 std::size_t iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
547 oAccessedArc[iArcIndex] =
true;
548 while (iNextArcIndex != iFirstArcIndex)
550 if (!AddArcToRing(iNextArcIndex))
554 iArcIndex = iNextArcIndex;
555 iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
556 oAccessedArc[iArcIndex] =
true;
560 poRing->closeRings();
563 poPolygon_->addRingDirectly(poNewRing.release());
567 for (
size_t i = 0; i < oAccessedArc.size(); ++i)
569 if (!oAccessedArc[i])
571 if (!AddRingToPolygon(i, poFirstRing))
576 poFirstRing =
nullptr;
582 if (iPixValField_ >= 0)
583 poFeature_->SetField(iPixValField_,
584 static_cast<double>(nPolygonCellValue));
587 if (poOutLayer_->CreateFeature(poFeature_.get()) !=
OGRERR_NONE)
592 if (nCommitInterval_ != 0)
594 ++nFeaturesWrittenInTransaction_;
595 if (nFeaturesWrittenInTransaction_ == nCommitInterval_)
597 if (poOutDS_->CommitTransaction() !=
OGRERR_NONE ||
602 nFeaturesWrittenInTransaction_ = 0;
607 if (poFeature_->GetGeometryRef() != poPolygon_)
610 poFeature_->SetGeometryDirectly(poPolygon_);
Class that saves the error state on construction, and restores it on destruction.
Definition cpl_error.h:410
This class represents a layer of simple features, with access methods.
Definition ogrsf_frmts.h:61
Concrete representation of a closed ring.
Definition ogr_geometry.h:1820
Concrete class representing polygons.
Definition ogr_geometry.h:2635
void empty() override
Clear geometry information.
Definition ogrlinestring.cpp:168
#define CPLAssert(expr)
Assert on an expression.
Definition cpl_error.h:330
#define CPLE_OutOfMemory
Out of memory error.
Definition cpl_error.h:87
void * OGRLayerH
Opaque type for a layer (OGRLayer).
Definition gdal_fwd.h:157
#define OGRNullFID
Special value for a unset FID.
Definition ogr_core.h:830
#define OGRERR_NONE
Success.
Definition ogr_core.h:372
#define ODsCTransactions
Dataset capability for dataset transactions.
Definition ogr_core.h:1032