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