Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: SDTS Translator
4 : * Purpose: Implementation of SDTSPolygonReader and SDTSRawPolygon classes.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1999, Frank Warmerdam
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "sdts_al.h"
14 :
15 : #include <cmath>
16 :
17 : /************************************************************************/
18 : /* ==================================================================== */
19 : /* SDTSRawPolygon */
20 : /* */
21 : /* This is a simple class for holding the data related with a */
22 : /* polygon feature. */
23 : /* ==================================================================== */
24 : /************************************************************************/
25 :
26 : /************************************************************************/
27 : /* SDTSRawPolygon() */
28 : /************************************************************************/
29 :
30 35 : SDTSRawPolygon::SDTSRawPolygon()
31 : : nEdges(0), papoEdges(nullptr), nRings(0), nVertices(0),
32 35 : panRingStart(nullptr), padfX(nullptr), padfY(nullptr), padfZ(nullptr)
33 : {
34 35 : nAttributes = 0;
35 35 : }
36 :
37 : /************************************************************************/
38 : /* ~SDTSRawPolygon() */
39 : /************************************************************************/
40 :
41 70 : SDTSRawPolygon::~SDTSRawPolygon()
42 :
43 : {
44 35 : CPLFree(papoEdges);
45 35 : CPLFree(panRingStart);
46 35 : CPLFree(padfX);
47 35 : CPLFree(padfY);
48 35 : CPLFree(padfZ);
49 70 : }
50 :
51 : /************************************************************************/
52 : /* Read() */
53 : /* */
54 : /* Read a record from the passed SDTSPolygonReader, and assign the */
55 : /* values from that record to this object. This is the bulk of */
56 : /* the work in this whole file. */
57 : /************************************************************************/
58 :
59 35 : int SDTSRawPolygon::Read(DDFRecord *poRecord)
60 :
61 : {
62 : /* ==================================================================== */
63 : /* Loop over fields in this record, looking for those we */
64 : /* recognise, and need. */
65 : /* ==================================================================== */
66 105 : for (int iField = 0; iField < poRecord->GetFieldCount(); iField++)
67 : {
68 70 : DDFField *poField = poRecord->GetField(iField);
69 70 : if (poField == nullptr)
70 0 : return FALSE;
71 70 : DDFFieldDefn *poFieldDefn = poField->GetFieldDefn();
72 70 : if (poFieldDefn == nullptr)
73 0 : return FALSE;
74 :
75 70 : const char *pszFieldName = poFieldDefn->GetName();
76 :
77 70 : if (EQUAL(pszFieldName, "POLY"))
78 : {
79 35 : oModId.Set(poField);
80 : }
81 :
82 35 : else if (EQUAL(pszFieldName, "ATID"))
83 : {
84 0 : ApplyATID(poField);
85 : }
86 : }
87 :
88 35 : return TRUE;
89 : }
90 :
91 : /************************************************************************/
92 : /* AddEdge() */
93 : /************************************************************************/
94 :
95 50 : void SDTSRawPolygon::AddEdge(SDTSRawLine *poNewLine)
96 :
97 : {
98 50 : nEdges++;
99 :
100 50 : papoEdges = reinterpret_cast<SDTSRawLine **>(
101 50 : CPLRealloc(papoEdges, sizeof(void *) * nEdges));
102 50 : papoEdges[nEdges - 1] = poNewLine;
103 50 : }
104 :
105 : /************************************************************************/
106 : /* AddEdgeToRing() */
107 : /************************************************************************/
108 :
109 52 : void SDTSRawPolygon::AddEdgeToRing(int nVertToAdd, double *padfXToAdd,
110 : double *padfYToAdd, double *padfZToAdd,
111 : int bReverse, int bDropVertex)
112 :
113 : {
114 52 : int iStart = 0;
115 52 : int iEnd = nVertToAdd - 1;
116 52 : int iStep = 1;
117 :
118 52 : if (bDropVertex && bReverse)
119 : {
120 4 : iStart = nVertToAdd - 2;
121 4 : iEnd = 0;
122 4 : iStep = -1;
123 : }
124 48 : else if (bDropVertex && !bReverse)
125 : {
126 29 : iStart = 1;
127 29 : iEnd = nVertToAdd - 1;
128 29 : iStep = 1;
129 : }
130 19 : else if (!bDropVertex && !bReverse)
131 : {
132 19 : iStart = 0;
133 19 : iEnd = nVertToAdd - 1;
134 19 : iStep = 1;
135 : }
136 0 : else if (!bDropVertex && bReverse)
137 : {
138 0 : iStart = nVertToAdd - 1;
139 0 : iEnd = 0;
140 0 : iStep = -1;
141 : }
142 :
143 1211 : for (int i = iStart; i != (iEnd + iStep); i += iStep)
144 : {
145 1159 : padfX[nVertices] = padfXToAdd[i];
146 1159 : padfY[nVertices] = padfYToAdd[i];
147 1159 : padfZ[nVertices] = padfZToAdd[i];
148 :
149 1159 : nVertices++;
150 : }
151 52 : }
152 :
153 : /************************************************************************/
154 : /* AssembleRings() */
155 : /************************************************************************/
156 :
157 : /**
158 : * Form border lines (arcs) into outer and inner rings.
159 : *
160 : * See SDTSPolygonReader::AssemblePolygons() for a simple one step process
161 : * to assembling geometry for all polygons in a transfer.
162 : *
163 : * This method will assemble the lines attached to a polygon into
164 : * an outer ring, and zero or more inner rings. Before calling it is
165 : * necessary that all the lines associated with this polygon have already
166 : * been attached. Normally this is accomplished by calling
167 : * SDTSLineReader::AttachToPolygons() on all line layers that might
168 : * contain edges related to this layer.
169 : *
170 : * This method then forms the lines into rings. Rings are formed by:
171 : * <ol>
172 : * <li> Take a previously unconsumed line, and start a ring with it. Mark
173 : * it as consumed, and keep track of its start and end node ids as
174 : * being the start and end node ids of the ring.
175 : * <li> If the rings start id is the same as the end node id then this ring
176 : * is completely formed, return to step 1.
177 : * <li> Search all unconsumed lines for a line with the same start or end
178 : * node id as the rings current node id. If none are found then the
179 : * assembly has failed. Return to step 1 but report failure on
180 : * completion.
181 : * <li> Once found, add the line to the current ring, dropping the duplicated
182 : * vertex and reverse order if necessary. Mark the line as consumed,
183 : * and update the rings end node id accordingly.
184 : * <li> go to step 2.
185 : * </ol>
186 : *
187 : * Once ring assembly from lines is complete, another pass is made to
188 : * order the rings such that the exterior ring is first, the first ring
189 : * has counter-clockwise vertex ordering and the inner rings have clockwise
190 : * vertex ordering. This is accomplished based on the assumption that the
191 : * outer ring has the largest area, and using the +/- sign of area to establish
192 : * direction of rings.
193 : *
194 : * @return TRUE if all rings assembled without problems or FALSE if a problem
195 : * occurred. If a problem occurs rings are still formed from all lines, but
196 : * some of the rings will not be closed, and rings will have no particular
197 : * order or direction.
198 : */
199 :
200 35 : int SDTSRawPolygon::AssembleRings()
201 :
202 : {
203 35 : if (nRings > 0)
204 0 : return TRUE;
205 :
206 35 : if (nEdges == 0)
207 22 : return FALSE;
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Setup array of line markers indicating if they have been */
211 : /* added to a ring yet. */
212 : /* -------------------------------------------------------------------- */
213 13 : int nRemainingEdges = nEdges;
214 :
215 : int *panEdgeConsumed =
216 13 : reinterpret_cast<int *>(CPLCalloc(sizeof(int), nEdges));
217 :
218 : /* -------------------------------------------------------------------- */
219 : /* Allocate ring arrays. */
220 : /* -------------------------------------------------------------------- */
221 13 : panRingStart = reinterpret_cast<int *>(CPLMalloc(sizeof(int) * nEdges));
222 :
223 13 : nVertices = 0;
224 63 : for (int iEdge = 0; iEdge < nEdges; iEdge++)
225 : {
226 50 : if (papoEdges[iEdge]->nVertices < 2)
227 : {
228 0 : panEdgeConsumed[iEdge] = TRUE;
229 0 : nRemainingEdges--;
230 : }
231 : else
232 : {
233 50 : nVertices += papoEdges[iEdge]->nVertices;
234 : }
235 : }
236 :
237 13 : padfX = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
238 13 : padfY = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
239 13 : padfZ = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
240 :
241 13 : nVertices = 0;
242 :
243 : /* ==================================================================== */
244 : /* Loop generating rings. */
245 : /* ==================================================================== */
246 13 : bool bSuccess = true;
247 :
248 30 : while (nRemainingEdges > 0)
249 : {
250 : /* --------------------------------------------------------------------
251 : */
252 : /* Find the first unconsumed edge. */
253 : /* --------------------------------------------------------------------
254 : */
255 17 : int iEdge = 0;
256 25 : for (; panEdgeConsumed[iEdge]; iEdge++)
257 : {
258 : }
259 :
260 17 : SDTSRawLine *poEdge = papoEdges[iEdge];
261 :
262 : /* --------------------------------------------------------------------
263 : */
264 : /* Start a new ring, copying in the current line directly */
265 : /* --------------------------------------------------------------------
266 : */
267 17 : panRingStart[nRings++] = nVertices;
268 :
269 17 : AddEdgeToRing(poEdge->nVertices, poEdge->padfX, poEdge->padfY,
270 : poEdge->padfZ, FALSE, FALSE);
271 :
272 17 : panEdgeConsumed[iEdge] = TRUE;
273 17 : nRemainingEdges--;
274 :
275 17 : const int nStartNode = poEdge->oStartNode.nRecord;
276 17 : int nLinkNode = poEdge->oEndNode.nRecord;
277 :
278 : /* ====================================================================
279 : */
280 : /* Loop adding edges to this ring until we make a whole pass */
281 : /* within finding anything to add. */
282 : /* ====================================================================
283 : */
284 17 : bool bWorkDone = true;
285 :
286 38 : while (nLinkNode != nStartNode && nRemainingEdges > 0 && bWorkDone)
287 : {
288 21 : bWorkDone = false;
289 :
290 319 : for (iEdge = 0; iEdge < nEdges; iEdge++)
291 : {
292 298 : if (panEdgeConsumed[iEdge])
293 182 : continue;
294 :
295 116 : poEdge = papoEdges[iEdge];
296 116 : if (poEdge->oStartNode.nRecord == nLinkNode)
297 : {
298 29 : AddEdgeToRing(poEdge->nVertices, poEdge->padfX,
299 : poEdge->padfY, poEdge->padfZ, FALSE, TRUE);
300 29 : nLinkNode = poEdge->oEndNode.nRecord;
301 : }
302 87 : else if (poEdge->oEndNode.nRecord == nLinkNode)
303 : {
304 4 : AddEdgeToRing(poEdge->nVertices, poEdge->padfX,
305 : poEdge->padfY, poEdge->padfZ, TRUE, TRUE);
306 4 : nLinkNode = poEdge->oStartNode.nRecord;
307 : }
308 : else
309 : {
310 83 : continue;
311 : }
312 :
313 33 : panEdgeConsumed[iEdge] = TRUE;
314 33 : nRemainingEdges--;
315 33 : bWorkDone = true;
316 : }
317 : }
318 :
319 : /* --------------------------------------------------------------------
320 : */
321 : /* Did we fail to complete the ring? */
322 : /* --------------------------------------------------------------------
323 : */
324 17 : if (nLinkNode != nStartNode)
325 15 : bSuccess = false;
326 : } /* next ring */
327 :
328 13 : CPLFree(panEdgeConsumed);
329 :
330 13 : if (!bSuccess)
331 11 : return bSuccess;
332 :
333 : /* ==================================================================== */
334 : /* Compute the area of each ring. The sign will be positive */
335 : /* for counter clockwise rings, otherwise negative. */
336 : /* */
337 : /* The algorithm used in this function was taken from _Graphics */
338 : /* Gems II_, James Arvo, 1991, Academic Press, Inc., section 1.1, */
339 : /* "The Area of a Simple Polygon", Jon Rokne, pp. 5-6. */
340 : /* ==================================================================== */
341 2 : double dfMaxArea = 0.0;
342 2 : int iBiggestRing = -1;
343 :
344 : double *padfRingArea =
345 2 : reinterpret_cast<double *>(CPLCalloc(sizeof(double), nRings));
346 :
347 4 : for (int iRing = 0; iRing < nRings; iRing++)
348 : {
349 : int nRingVertices;
350 2 : if (iRing == nRings - 1)
351 2 : nRingVertices = nVertices - panRingStart[iRing];
352 : else
353 0 : nRingVertices = panRingStart[iRing + 1] - panRingStart[iRing];
354 :
355 2 : double dfSum1 = 0.0;
356 2 : double dfSum2 = 0.0;
357 2 : for (int i = panRingStart[iRing];
358 382 : i < panRingStart[iRing] + nRingVertices - 1; i++)
359 : {
360 380 : dfSum1 += padfX[i] * padfY[i + 1];
361 380 : dfSum2 += padfY[i] * padfX[i + 1];
362 : }
363 :
364 2 : padfRingArea[iRing] = (dfSum1 - dfSum2) / 2;
365 :
366 2 : if (std::abs(padfRingArea[iRing]) > dfMaxArea)
367 : {
368 2 : dfMaxArea = std::abs(padfRingArea[iRing]);
369 2 : iBiggestRing = iRing;
370 : }
371 : }
372 :
373 2 : if (iBiggestRing < 0)
374 : {
375 0 : CPLFree(padfRingArea);
376 0 : return FALSE;
377 : }
378 :
379 : /* ==================================================================== */
380 : /* Make a new set of vertices, and copy the largest ring into */
381 : /* it, adjusting the direction if necessary to ensure that this */
382 : /* outer ring is counter clockwise. */
383 : /* ==================================================================== */
384 2 : double *padfXRaw = padfX;
385 2 : double *padfYRaw = padfY;
386 2 : double *padfZRaw = padfZ;
387 2 : int *panRawRingStart = panRingStart;
388 2 : int nRawVertices = nVertices;
389 2 : int nRawRings = nRings;
390 :
391 2 : padfX = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
392 2 : padfY = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
393 2 : padfZ = reinterpret_cast<double *>(CPLMalloc(sizeof(double) * nVertices));
394 2 : panRingStart = reinterpret_cast<int *>(CPLMalloc(sizeof(int) * nRawRings));
395 2 : nVertices = 0;
396 2 : nRings = 0;
397 :
398 : int nRingVertices;
399 2 : if (iBiggestRing == nRawRings - 1)
400 2 : nRingVertices = nRawVertices - panRawRingStart[iBiggestRing];
401 : else
402 0 : nRingVertices =
403 0 : panRawRingStart[iBiggestRing + 1] - panRawRingStart[iBiggestRing];
404 :
405 2 : panRingStart[nRings++] = 0;
406 2 : AddEdgeToRing(nRingVertices, padfXRaw + panRawRingStart[iBiggestRing],
407 2 : padfYRaw + panRawRingStart[iBiggestRing],
408 2 : padfZRaw + panRawRingStart[iBiggestRing],
409 2 : padfRingArea[iBiggestRing] < 0.0, FALSE);
410 :
411 : /* ==================================================================== */
412 : /* Add the rest of the rings, which must be holes, in clockwise */
413 : /* order. */
414 : /* ==================================================================== */
415 4 : for (int iRing = 0; iRing < nRawRings; iRing++)
416 : {
417 2 : if (iRing == iBiggestRing)
418 2 : continue;
419 :
420 0 : if (iRing == nRawRings - 1)
421 0 : nRingVertices = nRawVertices - panRawRingStart[iRing];
422 : else
423 0 : nRingVertices = panRawRingStart[iRing + 1] - panRawRingStart[iRing];
424 :
425 0 : panRingStart[nRings++] = nVertices;
426 0 : AddEdgeToRing(nRingVertices, padfXRaw + panRawRingStart[iRing],
427 0 : padfYRaw + panRawRingStart[iRing],
428 0 : padfZRaw + panRawRingStart[iRing],
429 0 : padfRingArea[iRing] > 0.0, FALSE);
430 : }
431 :
432 : /* -------------------------------------------------------------------- */
433 : /* Cleanup */
434 : /* -------------------------------------------------------------------- */
435 2 : CPLFree(padfXRaw);
436 2 : CPLFree(padfYRaw);
437 2 : CPLFree(padfZRaw);
438 2 : CPLFree(padfRingArea);
439 2 : CPLFree(panRawRingStart);
440 :
441 2 : CPLFree(papoEdges);
442 2 : papoEdges = nullptr;
443 2 : nEdges = 0;
444 :
445 2 : return TRUE;
446 : }
447 :
448 : /************************************************************************/
449 : /* Dump() */
450 : /************************************************************************/
451 :
452 0 : void SDTSRawPolygon::Dump(FILE *fp)
453 :
454 : {
455 0 : fprintf(fp, "SDTSRawPolygon %s: ", oModId.GetName());
456 :
457 0 : for (int i = 0; i < nAttributes; i++)
458 0 : fprintf(fp, " ATID[%d]=%s", i, paoATID[i].GetName());
459 :
460 0 : fprintf(fp, "\n");
461 0 : }
462 :
463 : /************************************************************************/
464 : /* ==================================================================== */
465 : /* SDTSPolygonReader */
466 : /* */
467 : /* This is the class used to read a Polygon module. */
468 : /* ==================================================================== */
469 : /************************************************************************/
470 :
471 : /************************************************************************/
472 : /* SDTSPolygonReader() */
473 : /************************************************************************/
474 :
475 1 : SDTSPolygonReader::SDTSPolygonReader() : bRingsAssembled(FALSE)
476 : {
477 1 : }
478 :
479 : /************************************************************************/
480 : /* ~SDTSPolygonReader() */
481 : /************************************************************************/
482 :
483 2 : SDTSPolygonReader::~SDTSPolygonReader()
484 : {
485 2 : }
486 :
487 : /************************************************************************/
488 : /* Close() */
489 : /************************************************************************/
490 :
491 0 : void SDTSPolygonReader::Close()
492 :
493 : {
494 0 : oDDFModule.Close();
495 0 : }
496 :
497 : /************************************************************************/
498 : /* Open() */
499 : /* */
500 : /* Open the requested line file, and prepare to start reading */
501 : /* data records. */
502 : /************************************************************************/
503 :
504 1 : int SDTSPolygonReader::Open(const char *pszFilename)
505 :
506 : {
507 1 : return oDDFModule.Open(pszFilename);
508 : }
509 :
510 : /************************************************************************/
511 : /* GetNextPolygon() */
512 : /* */
513 : /* Fetch the next feature as an STDSRawPolygon. */
514 : /************************************************************************/
515 :
516 36 : SDTSRawPolygon *SDTSPolygonReader::GetNextPolygon()
517 :
518 : {
519 : /* -------------------------------------------------------------------- */
520 : /* Read a record. */
521 : /* -------------------------------------------------------------------- */
522 36 : if (oDDFModule.GetFP() == nullptr)
523 0 : return nullptr;
524 :
525 36 : DDFRecord *poRecord = oDDFModule.ReadRecord();
526 :
527 36 : if (poRecord == nullptr)
528 1 : return nullptr;
529 :
530 : /* -------------------------------------------------------------------- */
531 : /* Transform into a Polygon feature. */
532 : /* -------------------------------------------------------------------- */
533 35 : SDTSRawPolygon *poRawPolygon = new SDTSRawPolygon();
534 :
535 35 : if (poRawPolygon->Read(poRecord))
536 : {
537 35 : return poRawPolygon;
538 : }
539 :
540 0 : delete poRawPolygon;
541 0 : return nullptr;
542 : }
543 :
544 : /************************************************************************/
545 : /* AssembleRings() */
546 : /************************************************************************/
547 :
548 : /**
549 : * Assemble geometry for a polygon transfer.
550 : *
551 : * This method takes care of attaching lines from all the line layers in
552 : * this transfer to this polygon layer, assembling the lines into rings on
553 : * the polygons, and then cleaning up unnecessary intermediate results.
554 : *
555 : * Currently this method will leave the line layers rewound to the beginning
556 : * but indexed, and the polygon layer rewound but indexed. In the future
557 : * it may restore reading positions, and possibly flush line indexes if they
558 : * were not previously indexed.
559 : *
560 : * This method does nothing if the rings have already been assembled on
561 : * this layer using this method.
562 : *
563 : * See SDTSRawPolygon::AssembleRings() for more information on how the lines
564 : * are assembled into rings.
565 : *
566 : * @param poTransfer the SDTSTransfer that this reader is a part of. Used
567 : * to get a list of line layers that might be needed.
568 : * @param iPolyLayer the polygon reader instance number, used to avoid
569 : * processing lines for other layers.
570 : */
571 :
572 38 : void SDTSPolygonReader::AssembleRings(SDTSTransfer *poTransfer, int iPolyLayer)
573 :
574 : {
575 38 : if (bRingsAssembled)
576 37 : return;
577 :
578 1 : bRingsAssembled = TRUE;
579 :
580 : /* -------------------------------------------------------------------- */
581 : /* To write polygons we need to build them from their related */
582 : /* arcs. We don't know off hand which arc (line) layers */
583 : /* contribute so we process all line layers, attaching them to */
584 : /* polygons as appropriate. */
585 : /* -------------------------------------------------------------------- */
586 9 : for (int iLineLayer = 0; iLineLayer < poTransfer->GetLayerCount();
587 : iLineLayer++)
588 : {
589 8 : if (poTransfer->GetLayerType(iLineLayer) != SLTLine)
590 7 : continue;
591 :
592 : SDTSLineReader *poLineReader = reinterpret_cast<SDTSLineReader *>(
593 1 : poTransfer->GetLayerIndexedReader(iLineLayer));
594 1 : if (poLineReader == nullptr)
595 0 : continue;
596 :
597 1 : poLineReader->AttachToPolygons(poTransfer, iPolyLayer);
598 1 : poLineReader->Rewind();
599 : }
600 :
601 1 : if (!IsIndexed())
602 0 : return;
603 :
604 : /* -------------------------------------------------------------------- */
605 : /* Scan all polygons indexed on this reader, and assemble their */
606 : /* rings. */
607 : /* -------------------------------------------------------------------- */
608 1 : Rewind();
609 :
610 1 : SDTSFeature *poFeature = nullptr;
611 36 : while ((poFeature = GetNextFeature()) != nullptr)
612 : {
613 35 : SDTSRawPolygon *poPoly = reinterpret_cast<SDTSRawPolygon *>(poFeature);
614 :
615 35 : poPoly->AssembleRings();
616 : }
617 :
618 1 : Rewind();
619 : }
|
