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 : * Permission is hereby granted, free of charge, to any person obtaining a
10 : * copy of this software and associated documentation files (the "Software"),
11 : * to deal in the Software without restriction, including without limitation
12 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13 : * and/or sell copies of the Software, and to permit persons to whom the
14 : * Software is furnished to do so, subject to the following conditions:
15 : *
16 : * The above copyright notice and this permission notice shall be included
17 : * in all copies or substantial portions of the Software.
18 : *
19 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 : * DEALINGS IN THE SOFTWARE.
26 : ****************************************************************************/
27 :
28 : /*! @cond Doxygen_Suppress */
29 :
30 : #include "polygonize_polygonizer.h"
31 :
32 : #include <algorithm>
33 :
34 : namespace gdal
35 : {
36 : namespace polygonizer
37 : {
38 :
39 3639 : RPolygon::~RPolygon()
40 : {
41 39629 : for (auto &arc : oArcs)
42 : {
43 35990 : delete arc;
44 : }
45 3639 : }
46 :
47 35990 : IndexedArc RPolygon::newArc(bool bFollowRighthand)
48 : {
49 35990 : Arc *poArc = new Arc();
50 35990 : std::size_t iArcIndex = oArcs.size();
51 35990 : oArcs.push_back(poArc);
52 35990 : oArcRighthandFollow.push_back(bFollowRighthand);
53 35990 : oArcConnections.push_back(iArcIndex);
54 35990 : return IndexedArc{poArc, iArcIndex};
55 : }
56 :
57 35990 : void RPolygon::setArcConnection(const IndexedArc &oArc,
58 : const IndexedArc &oNextArc)
59 : {
60 35990 : oArcConnections[oArc.iIndex] = oNextArc.iIndex;
61 35990 : }
62 :
63 423251 : void RPolygon::updateBottomRightPos(IndexType iRow, IndexType iCol)
64 : {
65 423251 : iBottomRightRow = iRow;
66 423251 : iBottomRightCol = iCol;
67 423251 : }
68 :
69 : /**
70 : * Process different kinds of Arm connections.
71 : */
72 423251 : static void ProcessArmConnections(TwoArm *poCurrent, TwoArm *poAbove,
73 : TwoArm *poLeft)
74 : {
75 423251 : poCurrent->poPolyInside->updateBottomRightPos(poCurrent->iRow,
76 : poCurrent->iCol);
77 423251 : poCurrent->bSolidVertical = poCurrent->poPolyInside != poLeft->poPolyInside;
78 423251 : poCurrent->bSolidHorizontal =
79 423251 : poCurrent->poPolyInside != poAbove->poPolyInside;
80 423251 : poCurrent->poPolyAbove = poAbove->poPolyInside;
81 423251 : poCurrent->poPolyLeft = poLeft->poPolyInside;
82 :
83 423251 : constexpr int BIT_CUR_HORIZ = 0;
84 423251 : constexpr int BIT_CUR_VERT = 1;
85 423251 : constexpr int BIT_LEFT = 2;
86 423251 : constexpr int BIT_ABOVE = 3;
87 :
88 423251 : const int nArmConnectionType =
89 423251 : (static_cast<int>(poAbove->bSolidVertical) << BIT_ABOVE) |
90 423251 : (static_cast<int>(poLeft->bSolidHorizontal) << BIT_LEFT) |
91 423251 : (static_cast<int>(poCurrent->bSolidVertical) << BIT_CUR_VERT) |
92 423251 : (static_cast<int>(poCurrent->bSolidHorizontal) << BIT_CUR_HORIZ);
93 :
94 423251 : constexpr int VIRTUAL = 0;
95 423251 : constexpr int SOLID = 1;
96 :
97 423251 : constexpr int ABOVE_VIRTUAL = VIRTUAL << BIT_ABOVE;
98 423251 : constexpr int ABOVE_SOLID = SOLID << BIT_ABOVE;
99 :
100 423251 : constexpr int LEFT_VIRTUAL = VIRTUAL << BIT_LEFT;
101 423251 : constexpr int LEFT_SOLID = SOLID << BIT_LEFT;
102 :
103 423251 : constexpr int CUR_VERT_VIRTUAL = VIRTUAL << BIT_CUR_VERT;
104 423251 : constexpr int CUR_VERT_SOLID = SOLID << BIT_CUR_VERT;
105 :
106 423251 : constexpr int CUR_HORIZ_VIRTUAL = VIRTUAL << BIT_CUR_HORIZ;
107 423251 : constexpr int CUR_HORIZ_SOLID = SOLID << BIT_CUR_HORIZ;
108 :
109 : /**
110 : * There are 12 valid connection types depending on the arm types(virtual or solid)
111 : * The following diagram illustrates these kinds of connection types, ⇢⇣ means virtual arm, →↓ means solid arm.
112 : * ⇣ ⇣ ⇣ ⇣ ↓
113 : * ⇢ → → → → ⇢ → → ⇢ →
114 : * ↓ ⇣ ↓ ↓ ⇣
115 : * type=3 type=5 type=6 type=7 type=9
116 : *
117 : * ↓ ↓ ↓ ↓ ↓
118 : * ⇢ ⇢ ⇢ → → ⇢ → → → ⇢
119 : * ↓ ↓ ⇣ ⇣ ↓
120 : * type=10 type=11 type=12 type=13 type=14
121 : *
122 : * ↓ ⇣
123 : * → → ⇢ ⇢
124 : * ↓ ⇣
125 : * type=15 type=0
126 : *
127 : * For each connection type, we may create new arc, ,
128 : * Depending on the connection type, we may do the following things:
129 : * 1. Create new arc. If the arc is closed to the inner polygon, it is called "Inner Arc", otherwise "Outer Arc"
130 : * 2. Pass an arc to the next arm.
131 : * 3. "Close" two arcs. If two arcs meet at the bottom right corner of a cell, close them by recording the arc connection.
132 : * 4. Add grid position(row, col) to an arc.
133 : */
134 :
135 423251 : switch (nArmConnectionType)
136 : {
137 358585 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
138 : CUR_HORIZ_VIRTUAL: // 0
139 : // nothing to do
140 358585 : break;
141 :
142 7175 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
143 : CUR_HORIZ_SOLID: // 3
144 : // add inner arcs
145 7175 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
146 7175 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
147 7175 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
148 7175 : poCurrent->oArcVerInner);
149 7175 : poCurrent->oArcVerInner.poArc->emplace_back(
150 7175 : Point{poCurrent->iRow, poCurrent->iCol});
151 :
152 : // add outer arcs
153 7175 : poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
154 7175 : poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
155 7175 : poAbove->poPolyInside->setArcConnection(poCurrent->oArcVerOuter,
156 7175 : poCurrent->oArcHorOuter);
157 7175 : poCurrent->oArcHorOuter.poArc->push_back(
158 7175 : Point{poCurrent->iRow, poCurrent->iCol});
159 :
160 7175 : break;
161 15831 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
162 : CUR_HORIZ_SOLID: // 5
163 : // pass arcs
164 15831 : poCurrent->oArcHorInner = poLeft->oArcHorInner;
165 15831 : poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
166 :
167 15831 : break;
168 8688 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
169 : CUR_HORIZ_VIRTUAL: // 6
170 : // pass arcs
171 8688 : poCurrent->oArcVerInner = poLeft->oArcHorOuter;
172 8688 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
173 8688 : poCurrent->oArcVerInner.poArc->push_back(
174 8688 : Point{poCurrent->iRow, poCurrent->iCol});
175 8688 : poCurrent->oArcVerOuter.poArc->push_back(
176 8688 : Point{poCurrent->iRow, poCurrent->iCol});
177 :
178 8688 : break;
179 942 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_SOLID |
180 : CUR_HORIZ_SOLID: // 7
181 : // pass arcs
182 942 : poCurrent->oArcHorOuter = poLeft->oArcHorOuter;
183 942 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
184 942 : poLeft->oArcHorInner.poArc->push_back(
185 942 : Point{poCurrent->iRow, poCurrent->iCol});
186 :
187 : // add inner arcs
188 942 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
189 942 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
190 942 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
191 942 : poCurrent->oArcVerInner);
192 942 : poCurrent->oArcVerInner.poArc->push_back(
193 942 : Point{poCurrent->iRow, poCurrent->iCol});
194 :
195 942 : break;
196 8671 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
197 : CUR_HORIZ_SOLID: // 9
198 : // pass arcs
199 8671 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
200 8671 : poCurrent->oArcHorInner = poAbove->oArcVerOuter;
201 8671 : poCurrent->oArcHorOuter.poArc->push_back(
202 8671 : Point{poCurrent->iRow, poCurrent->iCol});
203 8671 : poCurrent->oArcHorInner.poArc->push_back(
204 8671 : Point{poCurrent->iRow, poCurrent->iCol});
205 :
206 8671 : break;
207 11609 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
208 : CUR_HORIZ_VIRTUAL: // 10
209 : // pass arcs
210 11609 : poCurrent->oArcVerInner = poAbove->oArcVerInner;
211 11609 : poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
212 :
213 11609 : break;
214 976 : case ABOVE_SOLID | LEFT_VIRTUAL | CUR_VERT_SOLID |
215 : CUR_HORIZ_SOLID: // 11
216 : // pass arcs
217 976 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
218 976 : poCurrent->oArcVerOuter = poAbove->oArcVerOuter;
219 976 : poCurrent->oArcHorOuter.poArc->push_back(
220 976 : Point{poCurrent->iRow, poCurrent->iCol});
221 : // add inner arcs
222 976 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
223 976 : poCurrent->oArcHorInner = poCurrent->poPolyInside->newArc(false);
224 976 : poCurrent->poPolyInside->setArcConnection(poCurrent->oArcHorInner,
225 976 : poCurrent->oArcVerInner);
226 976 : poCurrent->oArcVerInner.poArc->push_back(
227 976 : Point{poCurrent->iRow, poCurrent->iCol});
228 :
229 976 : break;
230 7195 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
231 : CUR_HORIZ_VIRTUAL: // 12
232 : // close arcs
233 7195 : poLeft->oArcHorOuter.poArc->push_back(
234 7195 : Point{poCurrent->iRow, poCurrent->iCol});
235 7195 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
236 7195 : poAbove->oArcVerOuter);
237 : // close arcs
238 7195 : poAbove->oArcVerInner.poArc->push_back(
239 7195 : Point{poCurrent->iRow, poCurrent->iCol});
240 7195 : poCurrent->poPolyInside->setArcConnection(poAbove->oArcVerInner,
241 7195 : poLeft->oArcHorInner);
242 :
243 7195 : break;
244 939 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_VIRTUAL |
245 : CUR_HORIZ_SOLID: // 13
246 : // close arcs
247 939 : poLeft->oArcHorOuter.poArc->push_back(
248 939 : Point{poCurrent->iRow, poCurrent->iCol});
249 939 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
250 939 : poAbove->oArcVerOuter);
251 : // pass arcs
252 939 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
253 939 : poCurrent->oArcHorInner = poLeft->oArcHorInner;
254 939 : poCurrent->oArcHorOuter.poArc->push_back(
255 939 : Point{poCurrent->iRow, poCurrent->iCol});
256 :
257 939 : break;
258 939 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID |
259 : CUR_HORIZ_VIRTUAL: // 14
260 : // close arcs
261 939 : poLeft->oArcHorOuter.poArc->push_back(
262 939 : Point{poCurrent->iRow, poCurrent->iCol});
263 939 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
264 939 : poAbove->oArcVerOuter);
265 : // pass arcs
266 939 : poCurrent->oArcVerInner = poAbove->oArcVerInner;
267 939 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
268 939 : poCurrent->oArcVerOuter.poArc->push_back(
269 939 : Point{poCurrent->iRow, poCurrent->iCol});
270 :
271 939 : break;
272 1701 : case ABOVE_SOLID | LEFT_SOLID | CUR_VERT_SOLID | CUR_HORIZ_SOLID: // 15
273 : // Tow pixels of the main diagonal belong to the same polygon
274 1701 : if (poAbove->poPolyLeft == poCurrent->poPolyInside)
275 : {
276 : // pass arcs
277 42 : poCurrent->oArcVerInner = poLeft->oArcHorOuter;
278 42 : poCurrent->oArcHorInner = poAbove->oArcVerOuter;
279 42 : poCurrent->oArcVerInner.poArc->push_back(
280 42 : Point{poCurrent->iRow, poCurrent->iCol});
281 42 : poCurrent->oArcHorInner.poArc->push_back(
282 42 : Point{poCurrent->iRow, poCurrent->iCol});
283 : }
284 : else
285 : {
286 : // close arcs
287 1659 : poLeft->oArcHorOuter.poArc->push_back(
288 1659 : Point{poCurrent->iRow, poCurrent->iCol});
289 1659 : poLeft->poPolyAbove->setArcConnection(poLeft->oArcHorOuter,
290 1659 : poAbove->oArcVerOuter);
291 : // add inner arcs
292 1659 : poCurrent->oArcVerInner = poCurrent->poPolyInside->newArc(true);
293 : poCurrent->oArcHorInner =
294 1659 : poCurrent->poPolyInside->newArc(false);
295 1659 : poCurrent->poPolyInside->setArcConnection(
296 1659 : poCurrent->oArcHorInner, poCurrent->oArcVerInner);
297 1659 : poCurrent->oArcVerInner.poArc->push_back(
298 1659 : Point{poCurrent->iRow, poCurrent->iCol});
299 : }
300 :
301 : // Tow pixels of the secondary diagonal belong to the same polygon
302 1701 : if (poAbove->poPolyInside == poLeft->poPolyInside)
303 : {
304 : // close arcs
305 68 : poAbove->poPolyInside->setArcConnection(poAbove->oArcVerInner,
306 68 : poLeft->oArcHorInner);
307 68 : poAbove->oArcVerInner.poArc->push_back(
308 68 : Point{poCurrent->iRow, poCurrent->iCol});
309 : // add outer arcs
310 68 : poCurrent->oArcHorOuter = poAbove->poPolyInside->newArc(true);
311 68 : poCurrent->oArcVerOuter = poAbove->poPolyInside->newArc(false);
312 68 : poCurrent->oArcHorOuter.poArc->push_back(
313 68 : Point{poCurrent->iRow, poCurrent->iCol});
314 68 : poAbove->poPolyInside->setArcConnection(
315 68 : poCurrent->oArcVerOuter, poCurrent->oArcHorOuter);
316 : }
317 : else
318 : {
319 : // pass arcs
320 1633 : poCurrent->oArcHorOuter = poAbove->oArcVerInner;
321 1633 : poCurrent->oArcVerOuter = poLeft->oArcHorInner;
322 1633 : poCurrent->oArcHorOuter.poArc->push_back(
323 1633 : Point{poCurrent->iRow, poCurrent->iCol});
324 1633 : poCurrent->oArcVerOuter.poArc->push_back(
325 1633 : Point{poCurrent->iRow, poCurrent->iCol});
326 : }
327 :
328 1701 : break;
329 :
330 0 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_VIRTUAL |
331 : CUR_HORIZ_SOLID: // 1
332 : case ABOVE_VIRTUAL | LEFT_VIRTUAL | CUR_VERT_SOLID |
333 : CUR_HORIZ_VIRTUAL: // 2
334 : case ABOVE_VIRTUAL | LEFT_SOLID | CUR_VERT_VIRTUAL |
335 : CUR_HORIZ_VIRTUAL: // 4
336 : default:
337 : // Impossible case
338 0 : CPLAssert(false);
339 : break;
340 : }
341 423251 : }
342 :
343 : template <typename PolyIdType, typename DataType>
344 52 : Polygonizer<PolyIdType, DataType>::Polygonizer(
345 : PolyIdType nInvalidPolyId, PolygonReceiver<DataType> *poPolygonReceiver)
346 52 : : nInvalidPolyId_(nInvalidPolyId), poPolygonReceiver_(poPolygonReceiver)
347 : {
348 52 : poTheOuterPolygon_ = createPolygon(THE_OUTER_POLYGON_ID);
349 52 : }
350 :
351 : template <typename PolyIdType, typename DataType>
352 52 : Polygonizer<PolyIdType, DataType>::~Polygonizer()
353 : {
354 : // cppcheck-suppress constVariableReference
355 104 : for (auto &pair : oPolygonMap_)
356 : {
357 52 : delete pair.second;
358 : }
359 52 : }
360 :
361 : template <typename PolyIdType, typename DataType>
362 421761 : RPolygon *Polygonizer<PolyIdType, DataType>::getPolygon(PolyIdType nPolygonId)
363 : {
364 421761 : if (oPolygonMap_.count(nPolygonId) == 0)
365 : {
366 3587 : return createPolygon(nPolygonId);
367 : }
368 : else
369 : {
370 418174 : return oPolygonMap_[nPolygonId];
371 : }
372 : }
373 :
374 : template <typename PolyIdType, typename DataType>
375 : RPolygon *
376 3639 : Polygonizer<PolyIdType, DataType>::createPolygon(PolyIdType nPolygonId)
377 : {
378 3639 : auto polygon = new RPolygon();
379 3639 : oPolygonMap_[nPolygonId] = polygon;
380 3639 : return polygon;
381 : }
382 :
383 : template <typename PolyIdType, typename DataType>
384 3587 : void Polygonizer<PolyIdType, DataType>::destroyPolygon(PolyIdType nPolygonId)
385 : {
386 3587 : delete oPolygonMap_[nPolygonId];
387 3587 : oPolygonMap_.erase(nPolygonId);
388 3587 : }
389 :
390 : template <typename PolyIdType, typename DataType>
391 1490 : void Polygonizer<PolyIdType, DataType>::processLine(
392 : const PolyIdType *panThisLineId, const DataType *panLastLineVal,
393 : TwoArm *poThisLineArm, TwoArm *poLastLineArm, const IndexType nCurrentRow,
394 : const IndexType nCols)
395 : {
396 : TwoArm *poCurrent, *poAbove, *poLeft;
397 :
398 1490 : poCurrent = poThisLineArm + 1;
399 1490 : poCurrent->iRow = nCurrentRow;
400 1490 : poCurrent->iCol = 0;
401 1490 : poCurrent->poPolyInside = getPolygon(panThisLineId[0]);
402 1490 : poAbove = poLastLineArm + 1;
403 1490 : poLeft = poThisLineArm;
404 1490 : poLeft->poPolyInside = poTheOuterPolygon_;
405 1490 : ProcessArmConnections(poCurrent, poAbove, poLeft);
406 421761 : for (IndexType col = 1; col < nCols; ++col)
407 : {
408 420271 : IndexType iArmIndex = col + 1;
409 420271 : poCurrent = poThisLineArm + iArmIndex;
410 420271 : poCurrent->iRow = nCurrentRow;
411 420271 : poCurrent->iCol = col;
412 420271 : poCurrent->poPolyInside = getPolygon(panThisLineId[col]);
413 420271 : poAbove = poLastLineArm + iArmIndex;
414 420271 : poLeft = poThisLineArm + iArmIndex - 1;
415 420271 : ProcessArmConnections(poCurrent, poAbove, poLeft);
416 : }
417 1490 : poCurrent = poThisLineArm + nCols + 1;
418 1490 : poCurrent->iRow = nCurrentRow;
419 1490 : poCurrent->iCol = nCols;
420 1490 : poCurrent->poPolyInside = poTheOuterPolygon_;
421 1490 : poAbove = poLastLineArm + nCols + 1;
422 1490 : poAbove->poPolyInside = poTheOuterPolygon_;
423 1490 : poLeft = poThisLineArm + nCols;
424 1490 : ProcessArmConnections(poCurrent, poAbove, poLeft);
425 :
426 : /**
427 : *
428 : * Find those polygons haven't been processed on this line as we can be sure they are completed
429 : *
430 : */
431 2980 : std::vector<PolygonMapEntry> oCompletedPolygons;
432 34097 : for (auto &entry : oPolygonMap_)
433 : {
434 32607 : RPolygon *poPolygon = entry.second;
435 :
436 32607 : if (poPolygon->iBottomRightRow + 1 == nCurrentRow)
437 : {
438 3587 : oCompletedPolygons.push_back(entry);
439 : }
440 : }
441 : // cppcheck-suppress constVariableReference
442 5077 : for (auto &entry : oCompletedPolygons)
443 : {
444 3587 : PolyIdType nPolyId = entry.first;
445 3587 : RPolygon *poPolygon = entry.second;
446 :
447 : // emit valid polygon only
448 3587 : if (nPolyId != nInvalidPolyId_)
449 : {
450 3568 : poPolygonReceiver_->receive(
451 3568 : poPolygon, panLastLineVal[poPolygon->iBottomRightCol]);
452 : }
453 :
454 3587 : destroyPolygon(nPolyId);
455 : }
456 1490 : }
457 :
458 : template <typename DataType>
459 52 : OGRPolygonWriter<DataType>::OGRPolygonWriter(OGRLayerH hOutLayer,
460 : int iPixValField,
461 : double *padfGeoTransform)
462 : : PolygonReceiver<DataType>(), hOutLayer_(hOutLayer),
463 52 : iPixValField_(iPixValField), padfGeoTransform_(padfGeoTransform)
464 : {
465 52 : }
466 :
467 : template <typename DataType>
468 3568 : void OGRPolygonWriter<DataType>::receive(RPolygon *poPolygon,
469 : DataType nPolygonCellValue)
470 : {
471 7136 : std::vector<bool> oAccessedArc(poPolygon->oArcConnections.size(), false);
472 3568 : double *padfGeoTransform = padfGeoTransform_;
473 :
474 3568 : OGRGeometryH hPolygon = OGR_G_CreateGeometry(wkbPolygon);
475 :
476 7159 : auto AddRingToPolygon = [&poPolygon, &oAccessedArc, &hPolygon,
477 : padfGeoTransform](std::size_t iFirstArcIndex)
478 : {
479 3591 : OGRGeometryH hRing = OGR_G_CreateGeometry(wkbLinearRing);
480 :
481 3591 : auto AddArcToRing =
482 465442 : [&poPolygon, &hRing, padfGeoTransform](std::size_t iArcIndex)
483 : {
484 33504 : auto oArc = poPolygon->oArcs[iArcIndex];
485 33504 : bool bArcFollowRighthand =
486 33504 : poPolygon->oArcRighthandFollow[iArcIndex];
487 106490 : for (std::size_t i = 0; i < oArc->size(); ++i)
488 : {
489 92619 : const Point &oPixel =
490 19633 : (*oArc)[bArcFollowRighthand ? i : (oArc->size() - i - 1)];
491 :
492 218958 : const double dfX = padfGeoTransform[0] +
493 72986 : oPixel[1] * padfGeoTransform[1] +
494 72986 : oPixel[0] * padfGeoTransform[2];
495 218958 : const double dfY = padfGeoTransform[3] +
496 72986 : oPixel[1] * padfGeoTransform[4] +
497 72986 : oPixel[0] * padfGeoTransform[5];
498 :
499 72986 : OGR_G_AddPoint_2D(hRing, dfX, dfY);
500 : }
501 : };
502 :
503 3591 : AddArcToRing(iFirstArcIndex);
504 :
505 3591 : std::size_t iArcIndex = iFirstArcIndex;
506 3591 : std::size_t iNextArcIndex = poPolygon->oArcConnections[iArcIndex];
507 3591 : oAccessedArc[iArcIndex] = true;
508 33504 : while (iNextArcIndex != iFirstArcIndex)
509 : {
510 29913 : AddArcToRing(iNextArcIndex);
511 29913 : iArcIndex = iNextArcIndex;
512 29913 : iNextArcIndex = poPolygon->oArcConnections[iArcIndex];
513 29913 : oAccessedArc[iArcIndex] = true;
514 : }
515 :
516 : // close ring manually
517 3591 : OGR_G_AddPoint_2D(hRing, OGR_G_GetX(hRing, 0), OGR_G_GetY(hRing, 0));
518 :
519 3591 : OGR_G_AddGeometryDirectly(hPolygon, hRing);
520 : };
521 :
522 3568 : std::vector<bool>::iterator ite;
523 7159 : while ((ite = std::find_if_not(oAccessedArc.begin(), oAccessedArc.end(),
524 44394 : [](bool accessed) { return accessed; })) !=
525 14318 : oAccessedArc.end())
526 : {
527 3591 : AddRingToPolygon(ite - oAccessedArc.begin());
528 : }
529 :
530 : // Create the feature object
531 3568 : OGRFeatureH hFeat = OGR_F_Create(OGR_L_GetLayerDefn(hOutLayer_));
532 :
533 3568 : OGR_F_SetGeometryDirectly(hFeat, hPolygon);
534 :
535 3568 : if (iPixValField_ >= 0)
536 2975 : OGR_F_SetFieldDouble(hFeat, iPixValField_,
537 : static_cast<double>(nPolygonCellValue));
538 :
539 : // Write the to the layer.
540 3568 : if (OGR_L_CreateFeature(hOutLayer_, hFeat) != OGRERR_NONE)
541 0 : eErr_ = CE_Failure;
542 :
543 3568 : OGR_F_Destroy(hFeat);
544 3568 : }
545 :
546 : } // namespace polygonizer
547 : } // namespace gdal
548 :
549 : /*! @endcond */
|
