Line data Source code
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 :
12 : /*! @cond Doxygen_Suppress */
13 :
14 : #include "polygonize_polygonizer.h"
15 :
16 : #include <algorithm>
17 :
18 : namespace gdal
19 : {
20 : namespace polygonizer
21 : {
22 35990 : IndexedArc RPolygon::newArc(bool bFollowRighthand)
23 : {
24 35990 : const std::size_t iArcIndex = oArcs.size();
25 : const auto &oNewArc =
26 35990 : oArcs.emplace_back(static_cast<unsigned>(iArcIndex), bFollowRighthand);
27 35990 : return IndexedArc{oNewArc.poArc.get(), iArcIndex};
28 : }
29 :
30 35990 : void RPolygon::setArcConnection(const IndexedArc &oArc,
31 : const IndexedArc &oNextArc)
32 : {
33 35990 : oArcs[oArc.iIndex].nConnection = static_cast<unsigned>(oNextArc.iIndex);
34 35990 : }
35 :
36 423251 : void RPolygon::updateBottomRightPos(IndexType iRow, IndexType iCol)
37 : {
38 423251 : iBottomRightRow = iRow;
39 423251 : iBottomRightCol = iCol;
40 423251 : }
41 :
42 : /**
43 : * Process different kinds of Arm connections.
44 : */
45 423251 : static void ProcessArmConnections(TwoArm *poCurrent, TwoArm *poAbove,
46 : TwoArm *poLeft)
47 : {
48 423251 : poCurrent->poPolyInside->updateBottomRightPos(poCurrent->iRow,
49 : poCurrent->iCol);
50 423251 : poCurrent->bSolidVertical = poCurrent->poPolyInside != poLeft->poPolyInside;
51 423251 : poCurrent->bSolidHorizontal =
52 423251 : poCurrent->poPolyInside != poAbove->poPolyInside;
53 423251 : poCurrent->poPolyAbove = poAbove->poPolyInside;
54 423251 : poCurrent->poPolyLeft = poLeft->poPolyInside;
55 :
56 423251 : constexpr int BIT_CUR_HORIZ = 0;
57 423251 : constexpr int BIT_CUR_VERT = 1;
58 423251 : constexpr int BIT_LEFT = 2;
59 423251 : constexpr int BIT_ABOVE = 3;
60 :
61 423251 : const int nArmConnectionType =
62 423251 : (static_cast<int>(poAbove->bSolidVertical) << BIT_ABOVE) |
63 423251 : (static_cast<int>(poLeft->bSolidHorizontal) << BIT_LEFT) |
64 423251 : (static_cast<int>(poCurrent->bSolidVertical) << BIT_CUR_VERT) |
65 423251 : (static_cast<int>(poCurrent->bSolidHorizontal) << BIT_CUR_HORIZ);
66 :
67 423251 : constexpr int VIRTUAL = 0;
68 423251 : constexpr int SOLID = 1;
69 :
70 423251 : constexpr int ABOVE_VIRTUAL = VIRTUAL << BIT_ABOVE;
71 423251 : constexpr int ABOVE_SOLID = SOLID << BIT_ABOVE;
72 :
73 423251 : constexpr int LEFT_VIRTUAL = VIRTUAL << BIT_LEFT;
74 423251 : constexpr int LEFT_SOLID = SOLID << BIT_LEFT;
75 :
76 423251 : constexpr int CUR_VERT_VIRTUAL = VIRTUAL << BIT_CUR_VERT;
77 423251 : constexpr int CUR_VERT_SOLID = SOLID << BIT_CUR_VERT;
78 :
79 423251 : constexpr int CUR_HORIZ_VIRTUAL = VIRTUAL << BIT_CUR_HORIZ;
80 423251 : constexpr int CUR_HORIZ_SOLID = SOLID << BIT_CUR_HORIZ;
81 :
82 : /**
83 : * There are 12 valid connection types depending on the arm types(virtual or solid)
84 : * The following diagram illustrates these kinds of connection types, ⇢⇣ means virtual arm, →↓ means solid arm.
85 : * ⇣ ⇣ ⇣ ⇣ ↓
86 : * ⇢ → → → → ⇢ → → ⇢ →
87 : * ↓ ⇣ ↓ ↓ ⇣
88 : * type=3 type=5 type=6 type=7 type=9
89 : *
90 : * ↓ ↓ ↓ ↓ ↓
91 : * ⇢ ⇢ ⇢ → → ⇢ → → → ⇢
92 : * ↓ ↓ ⇣ ⇣ ↓
93 : * type=10 type=11 type=12 type=13 type=14
94 : *
95 : * ↓ ⇣
96 : * → → ⇢ ⇢
97 : * ↓ ⇣
98 : * type=15 type=0
99 : *
100 : * For each connection type, we may create new arc, ,
101 : * Depending on the connection type, we may do the following things:
102 : * 1. Create new arc. If the arc is closed to the inner polygon, it is called "Inner Arc", otherwise "Outer Arc"
103 : * 2. Pass an arc to the next arm.
104 : * 3. "Close" two arcs. If two arcs meet at the bottom right corner of a cell, close them by recording the arc connection.
105 : * 4. Add grid position(row, col) to an arc.
106 : */
107 :
108 423251 : switch (nArmConnectionType)
109 : {
110 358585 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
111 : CUR_HORIZ_VIRTUAL: // 0
112 : // nothing to do
113 358585 : break;
114 :
115 7175 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
116 : CUR_HORIZ_SOLID: // 3
117 : // add inner arcs
118 7175 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
119 7175 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
120 7175 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
121 7175 : poCurrent->oArcVerInner);
122 7175 : poCurrent->oArcVerInner.poArc->emplace_back(
123 7175 : Point{poCurrent->iRow, poCurrent->iCol});
124 :
125 : // add outer arcs
126 7175 : poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
127 7175 : poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
128 7175 : poAbove->poPolyInside->setArcConnection(poCurrent->oArcVerOuter,
129 7175 : poCurrent->oArcHorOuter);
130 7175 : poCurrent->oArcHorOuter.poArc->push_back(
131 7175 : Point{poCurrent->iRow, poCurrent->iCol});
132 :
133 7175 : break;
134 15831 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
135 : CUR_HORIZ_SOLID: // 5
136 : // pass arcs
137 15831 : poCurrent->oArcHorInner = poLeft->oArcHorInner;
138 15831 : poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
139 :
140 15831 : break;
141 8688 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
142 : CUR_HORIZ_VIRTUAL: // 6
143 : // pass arcs
144 8688 : poCurrent->oArcVerInner = poLeft->oArcHorOuter;
145 8688 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
146 8688 : poCurrent->oArcVerInner.poArc->push_back(
147 8688 : Point{poCurrent->iRow, poCurrent->iCol});
148 8688 : poCurrent->oArcVerOuter.poArc->push_back(
149 8688 : Point{poCurrent->iRow, poCurrent->iCol});
150 :
151 8688 : break;
152 942 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
153 : CUR_HORIZ_SOLID: // 7
154 : // pass arcs
155 942 : poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
156 942 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
157 942 : poLeft->oArcHorInner.poArc->push_back(
158 942 : Point{poCurrent->iRow, poCurrent->iCol});
159 :
160 : // add inner arcs
161 942 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
162 942 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
163 942 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
164 942 : poCurrent->oArcVerInner);
165 942 : poCurrent->oArcVerInner.poArc->push_back(
166 942 : Point{poCurrent->iRow, poCurrent->iCol});
167 :
168 942 : break;
169 8671 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
170 : CUR_HORIZ_SOLID: // 9
171 : // pass arcs
172 8671 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
173 8671 : poCurrent->oArcHorInner = poAbove->oArcVerOuter;
174 8671 : poCurrent->oArcHorOuter.poArc->push_back(
175 8671 : Point{poCurrent->iRow, poCurrent->iCol});
176 8671 : poCurrent->oArcHorInner.poArc->push_back(
177 8671 : Point{poCurrent->iRow, poCurrent->iCol});
178 :
179 8671 : break;
180 11609 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
181 : CUR_HORIZ_VIRTUAL: // 10
182 : // pass arcs
183 11609 : poCurrent->oArcVerInner = poAbove->oArcVerInner;
184 11609 : poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
185 :
186 11609 : break;
187 976 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
188 : CUR_HORIZ_SOLID: // 11
189 : // pass arcs
190 976 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
191 976 : poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
192 976 : poCurrent->oArcHorOuter.poArc->push_back(
193 976 : Point{poCurrent->iRow, poCurrent->iCol});
194 : // add inner arcs
195 976 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
196 976 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
197 976 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
198 976 : poCurrent->oArcVerInner);
199 976 : poCurrent->oArcVerInner.poArc->push_back(
200 976 : Point{poCurrent->iRow, poCurrent->iCol});
201 :
202 976 : break;
203 7195 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
204 : CUR_HORIZ_VIRTUAL: // 12
205 : // close arcs
206 7195 : poLeft->oArcHorOuter.poArc->push_back(
207 7195 : Point{poCurrent->iRow, poCurrent->iCol});
208 7195 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
209 7195 : poAbove->oArcVerOuter);
210 : // close arcs
211 7195 : poAbove->oArcVerInner.poArc->push_back(
212 7195 : Point{poCurrent->iRow, poCurrent->iCol});
213 7195 : poCurrent->poPolyInside->setArcConnection(poAbove->oArcVerInner,
214 7195 : poLeft->oArcHorInner);
215 :
216 7195 : break;
217 939 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
218 : CUR_HORIZ_SOLID: // 13
219 : // close arcs
220 939 : poLeft->oArcHorOuter.poArc->push_back(
221 939 : Point{poCurrent->iRow, poCurrent->iCol});
222 939 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
223 939 : poAbove->oArcVerOuter);
224 : // pass arcs
225 939 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
226 939 : poCurrent->oArcHorInner = poLeft->oArcHorInner;
227 939 : poCurrent->oArcHorOuter.poArc->push_back(
228 939 : Point{poCurrent->iRow, poCurrent->iCol});
229 :
230 939 : break;
231 939 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID |
232 : CUR_HORIZ_VIRTUAL: // 14
233 : // close arcs
234 939 : poLeft->oArcHorOuter.poArc->push_back(
235 939 : Point{poCurrent->iRow, poCurrent->iCol});
236 939 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
237 939 : poAbove->oArcVerOuter);
238 : // pass arcs
239 939 : poCurrent->oArcVerInner = poAbove->oArcVerInner;
240 939 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
241 939 : poCurrent->oArcVerOuter.poArc->push_back(
242 939 : Point{poCurrent->iRow, poCurrent->iCol});
243 :
244 939 : break;
245 1701 : 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 1701 : if (poAbove->poPolyLeft == poCurrent->poPolyInside)
248 : {
249 : // pass arcs
250 42 : poCurrent->oArcVerInner = poLeft->oArcHorOuter;
251 42 : poCurrent->oArcHorInner = poAbove->oArcVerOuter;
252 42 : poCurrent->oArcVerInner.poArc->push_back(
253 42 : Point{poCurrent->iRow, poCurrent->iCol});
254 42 : poCurrent->oArcHorInner.poArc->push_back(
255 42 : Point{poCurrent->iRow, poCurrent->iCol});
256 : }
257 : else
258 : {
259 : // close arcs
260 1659 : poLeft->oArcHorOuter.poArc->push_back(
261 1659 : Point{poCurrent->iRow, poCurrent->iCol});
262 1659 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
263 1659 : poAbove->oArcVerOuter);
264 : // add inner arcs
265 1659 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
266 : poCurrent->oArcHorInner =
267 1659 : poCurrent->poPolyInside->newArc(false);
268 1659 : poCurrent->poPolyInside->setArcConnection(
269 1659 : poCurrent->oArcHorInner, poCurrent->oArcVerInner);
270 1659 : poCurrent->oArcVerInner.poArc->push_back(
271 1659 : Point{poCurrent->iRow, poCurrent->iCol});
272 : }
273 :
274 : // Tow pixels of the secondary diagonal belong to the same polygon
275 1701 : if (poAbove->poPolyInside == poLeft->poPolyInside)
276 : {
277 : // close arcs
278 68 : poAbove->poPolyInside->setArcConnection(poAbove->oArcVerInner,
279 68 : poLeft->oArcHorInner);
280 68 : poAbove->oArcVerInner.poArc->push_back(
281 68 : Point{poCurrent->iRow, poCurrent->iCol});
282 : // add outer arcs
283 68 : poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
284 68 : poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
285 68 : poCurrent->oArcHorOuter.poArc->push_back(
286 68 : Point{poCurrent->iRow, poCurrent->iCol});
287 68 : poAbove->poPolyInside->setArcConnection(
288 68 : poCurrent->oArcVerOuter, poCurrent->oArcHorOuter);
289 : }
290 : else
291 : {
292 : // pass arcs
293 1633 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
294 1633 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
295 1633 : poCurrent->oArcHorOuter.poArc->push_back(
296 1633 : Point{poCurrent->iRow, poCurrent->iCol});
297 1633 : poCurrent->oArcVerOuter.poArc->push_back(
298 1633 : Point{poCurrent->iRow, poCurrent->iCol});
299 : }
300 :
301 1701 : break;
302 :
303 0 : 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 0 : CPLAssert(false);
312 : break;
313 : }
314 423251 : }
315 :
316 : template <typename PolyIdType, typename DataType>
317 52 : Polygonizer<PolyIdType, DataType>::Polygonizer(
318 : PolyIdType nInvalidPolyId, PolygonReceiver<DataType> *poPolygonReceiver)
319 52 : : nInvalidPolyId_(nInvalidPolyId), poPolygonReceiver_(poPolygonReceiver)
320 : {
321 52 : poTheOuterPolygon_ = createPolygon(THE_OUTER_POLYGON_ID);
322 52 : }
323 :
324 : template <typename PolyIdType, typename DataType>
325 52 : Polygonizer<PolyIdType, DataType>::~Polygonizer()
326 : {
327 : // cppcheck-suppress constVariableReference
328 104 : for (auto &pair : oPolygonMap_)
329 : {
330 52 : delete pair.second;
331 : }
332 52 : }
333 :
334 : template <typename PolyIdType, typename DataType>
335 421761 : RPolygon *Polygonizer<PolyIdType, DataType>::getPolygon(PolyIdType nPolygonId)
336 : {
337 421761 : const auto oIter = oPolygonMap_.find(nPolygonId);
338 421761 : if (oIter == oPolygonMap_.end())
339 : {
340 3587 : return createPolygon(nPolygonId);
341 : }
342 : else
343 : {
344 418174 : return oIter->second;
345 : }
346 : }
347 :
348 : template <typename PolyIdType, typename DataType>
349 : RPolygon *
350 3639 : Polygonizer<PolyIdType, DataType>::createPolygon(PolyIdType nPolygonId)
351 : {
352 3639 : auto polygon = new RPolygon();
353 3639 : oPolygonMap_[nPolygonId] = polygon;
354 3639 : return polygon;
355 : }
356 :
357 : template <typename PolyIdType, typename DataType>
358 3587 : void Polygonizer<PolyIdType, DataType>::destroyPolygon(PolyIdType nPolygonId)
359 : {
360 3587 : const auto oIter = oPolygonMap_.find(nPolygonId);
361 3587 : CPLAssert(oIter != oPolygonMap_.end());
362 3587 : delete oIter->second;
363 3587 : oPolygonMap_.erase(oIter);
364 3587 : }
365 :
366 : template <typename PolyIdType, typename DataType>
367 1490 : bool 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 1490 : poCurrent = poThisLineArm + 1;
377 1490 : poCurrent->iRow = nCurrentRow;
378 1490 : poCurrent->iCol = 0;
379 1490 : poCurrent->poPolyInside = getPolygon(panThisLineId[0]);
380 1490 : poAbove = poLastLineArm + 1;
381 1490 : poLeft = poThisLineArm;
382 1490 : poLeft->poPolyInside = poTheOuterPolygon_;
383 1490 : ProcessArmConnections(poCurrent, poAbove, poLeft);
384 421761 : for (IndexType col = 1; col < nCols; ++col)
385 : {
386 420271 : IndexType iArmIndex = col + 1;
387 420271 : poCurrent = poThisLineArm + iArmIndex;
388 420271 : poCurrent->iRow = nCurrentRow;
389 420271 : poCurrent->iCol = col;
390 420271 : poCurrent->poPolyInside = getPolygon(panThisLineId[col]);
391 420271 : poAbove = poLastLineArm + iArmIndex;
392 420271 : poLeft = poThisLineArm + iArmIndex - 1;
393 420271 : ProcessArmConnections(poCurrent, poAbove, poLeft);
394 : }
395 1490 : poCurrent = poThisLineArm + nCols + 1;
396 1490 : poCurrent->iRow = nCurrentRow;
397 1490 : poCurrent->iCol = nCols;
398 1490 : poCurrent->poPolyInside = poTheOuterPolygon_;
399 1490 : poAbove = poLastLineArm + nCols + 1;
400 1490 : poAbove->poPolyInside = poTheOuterPolygon_;
401 1490 : poLeft = poThisLineArm + nCols;
402 1490 : ProcessArmConnections(poCurrent, poAbove, poLeft);
403 :
404 : /**
405 : *
406 : * Find those polygons haven't been processed on this line as we can be sure they are completed
407 : *
408 : */
409 1490 : std::vector<PolygonMapEntry> oCompletedPolygons;
410 34097 : for (auto &entry : oPolygonMap_)
411 : {
412 32607 : RPolygon *poPolygon = entry.second;
413 :
414 32607 : if (poPolygon->iBottomRightRow + 1 == nCurrentRow)
415 : {
416 3587 : oCompletedPolygons.push_back(entry);
417 : }
418 : }
419 : // cppcheck-suppress constVariableReference
420 5077 : for (auto &entry : oCompletedPolygons)
421 : {
422 3587 : PolyIdType nPolyId = entry.first;
423 3587 : RPolygon *poPolygon = entry.second;
424 :
425 : // emit valid polygon only
426 3587 : if (nPolyId != nInvalidPolyId_)
427 : {
428 3568 : poPolygonReceiver_->receive(
429 3568 : poPolygon, panLastLineVal[poPolygon->iBottomRightCol]);
430 : }
431 :
432 3587 : destroyPolygon(nPolyId);
433 : }
434 1490 : return true;
435 : }
436 0 : catch (const std::bad_alloc &)
437 : {
438 0 : CPLError(CE_Failure, CPLE_OutOfMemory,
439 : "Out of memory in Polygonizer::processLine");
440 0 : return false;
441 : }
442 : }
443 :
444 : template <typename DataType>
445 52 : OGRPolygonWriter<DataType>::OGRPolygonWriter(OGRLayerH hOutLayer,
446 : int iPixValField,
447 : double *padfGeoTransform)
448 52 : : PolygonReceiver<DataType>(), poOutLayer_(OGRLayer::FromHandle(hOutLayer)),
449 52 : iPixValField_(iPixValField), padfGeoTransform_(padfGeoTransform)
450 : {
451 52 : poFeature_ = std::make_unique<OGRFeature>(poOutLayer_->GetLayerDefn());
452 52 : poPolygon_ = new OGRPolygon();
453 52 : poFeature_->SetGeometryDirectly(poPolygon_);
454 52 : }
455 :
456 : template <typename DataType>
457 3568 : void OGRPolygonWriter<DataType>::receive(RPolygon *poPolygon,
458 : DataType nPolygonCellValue)
459 : {
460 3568 : std::vector<bool> oAccessedArc(poPolygon->oArcs.size(), false);
461 3568 : double *padfGeoTransform = padfGeoTransform_;
462 :
463 3568 : OGRLinearRing *poFirstRing = poPolygon_->getExteriorRing();
464 3568 : if (poFirstRing && poPolygon_->getNumInteriorRings() == 0)
465 : {
466 3504 : poFirstRing->empty();
467 : }
468 : else
469 : {
470 64 : poFirstRing = nullptr;
471 64 : poPolygon_->empty();
472 : }
473 :
474 3568 : auto AddRingToPolygon =
475 74277 : [this, &poPolygon, &oAccessedArc,
476 : padfGeoTransform](std::size_t iFirstArcIndex, OGRLinearRing *poRing)
477 : {
478 3591 : std::unique_ptr<OGRLinearRing> poNewRing;
479 3591 : if (!poRing)
480 : {
481 87 : poNewRing = std::make_unique<OGRLinearRing>();
482 87 : poRing = poNewRing.get();
483 : }
484 :
485 3591 : auto AddArcToRing =
486 498946 : [&poPolygon, poRing, padfGeoTransform](std::size_t iArcIndex)
487 : {
488 33504 : const auto &oArc = poPolygon->oArcs[iArcIndex];
489 33504 : const bool bArcFollowRighthand = oArc.bFollowRighthand;
490 33504 : const int nArcPointCount = static_cast<int>(oArc.poArc->size());
491 33504 : int nDstPointIdx = poRing->getNumPoints();
492 33504 : poRing->setNumPoints(nDstPointIdx + nArcPointCount,
493 : /* bZeroizeNewContent = */ false);
494 33504 : if (poRing->getNumPoints() < nDstPointIdx + nArcPointCount)
495 : {
496 0 : return false;
497 : }
498 106490 : for (int i = 0; i < nArcPointCount; ++i)
499 : {
500 92619 : const Point &oPixel =
501 72986 : (*oArc.poArc)[bArcFollowRighthand
502 : ? i
503 19633 : : (nArcPointCount - i - 1)];
504 :
505 218958 : const double dfX = padfGeoTransform[0] +
506 72986 : oPixel[1] * padfGeoTransform[1] +
507 72986 : oPixel[0] * padfGeoTransform[2];
508 218958 : const double dfY = padfGeoTransform[3] +
509 72986 : oPixel[1] * padfGeoTransform[4] +
510 72986 : oPixel[0] * padfGeoTransform[5];
511 :
512 72986 : poRing->setPoint(nDstPointIdx, dfX, dfY);
513 72986 : ++nDstPointIdx;
514 : }
515 33504 : return true;
516 : };
517 :
518 3591 : if (!AddArcToRing(iFirstArcIndex))
519 : {
520 0 : return false;
521 : }
522 :
523 3591 : std::size_t iArcIndex = iFirstArcIndex;
524 3591 : std::size_t iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
525 3591 : oAccessedArc[iArcIndex] = true;
526 33504 : while (iNextArcIndex != iFirstArcIndex)
527 : {
528 29913 : if (!AddArcToRing(iNextArcIndex))
529 : {
530 0 : return false;
531 : }
532 29913 : iArcIndex = iNextArcIndex;
533 29913 : iNextArcIndex = poPolygon->oArcs[iArcIndex].nConnection;
534 29913 : oAccessedArc[iArcIndex] = true;
535 : }
536 :
537 : // close ring manually
538 3591 : poRing->closeRings();
539 :
540 3591 : if (poNewRing)
541 87 : poPolygon_->addRingDirectly(poNewRing.release());
542 3591 : return true;
543 : };
544 :
545 37072 : for (size_t i = 0; i < oAccessedArc.size(); ++i)
546 : {
547 33504 : if (!oAccessedArc[i])
548 : {
549 3591 : if (!AddRingToPolygon(i, poFirstRing))
550 : {
551 0 : eErr_ = CE_Failure;
552 0 : return;
553 : }
554 3591 : poFirstRing = nullptr;
555 : }
556 : }
557 :
558 : // Create the feature object
559 3568 : poFeature_->SetFID(OGRNullFID);
560 3568 : if (iPixValField_ >= 0)
561 2975 : poFeature_->SetField(iPixValField_,
562 : static_cast<double>(nPolygonCellValue));
563 :
564 : // Write the to the layer.
565 3568 : if (poOutLayer_->CreateFeature(poFeature_.get()) != OGRERR_NONE)
566 0 : eErr_ = CE_Failure;
567 :
568 : // Shouldn't happen for well behaved drivers, but better check...
569 3568 : else if (poFeature_->GetGeometryRef() != poPolygon_)
570 : {
571 0 : poPolygon_ = new OGRPolygon();
572 0 : poFeature_->SetGeometryDirectly(poPolygon_);
573 : }
574 : }
575 :
576 : } // namespace polygonizer
577 : } // namespace gdal
578 :
579 : /*! @endcond */
|
