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