1 : /******************************************************************************
2 : * $Id: ograssemblepolygon.cpp 19790 2010-06-01 15:46:20Z warmerdam $
3 : *
4 : * Project: S-57 Reader
5 : * Purpose: Implements polygon assembly from a bunch of arcs.
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 "ogr_geometry.h"
31 : #include "ogr_api.h"
32 : #include "cpl_conv.h"
33 :
34 : CPL_CVSID("$Id: ograssemblepolygon.cpp 19790 2010-06-01 15:46:20Z warmerdam $");
35 :
36 : /************************************************************************/
37 : /* CheckPoints() */
38 : /* */
39 : /* Check if two points are closer than the current best */
40 : /* distance. Update the current best distance if they are. */
41 : /************************************************************************/
42 :
43 : static int CheckPoints( OGRLineString *poLine1, int iPoint1,
44 : OGRLineString *poLine2, int iPoint2,
45 250 : double *pdfDistance )
46 :
47 : {
48 : double dfDeltaX, dfDeltaY, dfDistance;
49 :
50 250 : if( pdfDistance == NULL || *pdfDistance == 0 )
51 : return poLine1->getX(iPoint1) == poLine2->getX(iPoint2)
52 250 : && poLine1->getY(iPoint1) == poLine2->getY(iPoint2);
53 :
54 0 : dfDeltaX = poLine1->getX(iPoint1) - poLine2->getX(iPoint2);
55 0 : dfDeltaY = poLine1->getY(iPoint1) - poLine2->getY(iPoint2);
56 :
57 0 : dfDeltaX = ABS(dfDeltaX);
58 0 : dfDeltaY = ABS(dfDeltaY);
59 :
60 0 : if( dfDeltaX > *pdfDistance || dfDeltaY > *pdfDistance )
61 0 : return FALSE;
62 :
63 0 : dfDistance = sqrt(dfDeltaX*dfDeltaX + dfDeltaY*dfDeltaY);
64 :
65 0 : if( dfDistance < *pdfDistance )
66 : {
67 0 : *pdfDistance = dfDistance;
68 0 : return TRUE;
69 : }
70 : else
71 0 : return FALSE;
72 : }
73 :
74 : /************************************************************************/
75 : /* AddEdgeToRing() */
76 : /************************************************************************/
77 :
78 : static void AddEdgeToRing( OGRLinearRing * poRing, OGRLineString * poLine,
79 25 : int bReverse )
80 :
81 : {
82 : /* -------------------------------------------------------------------- */
83 : /* Establish order and range of traverse. */
84 : /* -------------------------------------------------------------------- */
85 25 : int iStart=0, iEnd=0, iStep=0;
86 25 : int nVertToAdd = poLine->getNumPoints();
87 :
88 25 : if( !bReverse )
89 : {
90 19 : iStart = 0;
91 19 : iEnd = nVertToAdd - 1;
92 19 : iStep = 1;
93 : }
94 : else
95 : {
96 6 : iStart = nVertToAdd - 1;
97 6 : iEnd = 0;
98 6 : iStep = -1;
99 : }
100 :
101 : /* -------------------------------------------------------------------- */
102 : /* Should we skip a repeating vertex? */
103 : /* -------------------------------------------------------------------- */
104 25 : if( poRing->getNumPoints() > 0
105 : && CheckPoints( poRing, poRing->getNumPoints()-1,
106 : poLine, iStart, NULL ) )
107 : {
108 22 : iStart += iStep;
109 : }
110 :
111 25 : poRing->addSubLineString( poLine, iStart, iEnd );
112 25 : }
113 :
114 : /************************************************************************/
115 : /* OGRBuildPolygonFromEdges() */
116 : /************************************************************************/
117 :
118 : /**
119 : * Build a ring from a bunch of arcs.
120 : *
121 : * @param hLines handle to an OGRGeometryCollection (or OGRMultiLineString) containing the line string geometries to be built into rings.
122 : * @param bBestEffort not yet implemented???.
123 : * @param bAutoClose indicates if the ring should be close when first and
124 : * last points of the ring are the same.
125 : * @param dfTolerance tolerance into which two arcs are considered
126 : * close enough to be joined.
127 : * @param peErr OGRERR_NONE on success, or OGRERR_FAILURE on failure.
128 : * @return an handle to the new geometry, a polygon.
129 : *
130 : */
131 :
132 : OGRGeometryH OGRBuildPolygonFromEdges( OGRGeometryH hLines,
133 : int bBestEffort,
134 : int bAutoClose,
135 : double dfTolerance,
136 5 : OGRErr * peErr )
137 :
138 : {
139 5 : if( hLines == NULL )
140 : {
141 0 : if (peErr != NULL)
142 0 : *peErr = OGRERR_NONE;
143 0 : return NULL;
144 : }
145 :
146 : /* -------------------------------------------------------------------- */
147 : /* Check for the case of a geometrycollection that can be */
148 : /* promoted to MultiLineString. */
149 : /* -------------------------------------------------------------------- */
150 5 : OGRGeometry* poGeom = (OGRGeometry*) hLines;
151 5 : if( wkbFlatten(poGeom->getGeometryType()) == wkbGeometryCollection )
152 : {
153 : int iGeom;
154 3 : OGRGeometryCollection *poGC = (OGRGeometryCollection *) poGeom;
155 :
156 23 : for( iGeom = 0; iGeom < poGC->getNumGeometries(); iGeom++ )
157 : {
158 21 : if( wkbFlatten(poGC->getGeometryRef(iGeom)->getGeometryType())
159 : != wkbLineString )
160 : {
161 1 : if (peErr != NULL)
162 1 : *peErr = OGRERR_FAILURE;
163 : CPLError(CE_Failure, CPLE_NotSupported,
164 1 : "The geometry collection contains non line string geometries");
165 1 : return NULL;
166 : }
167 : }
168 : }
169 2 : else if (wkbFlatten(poGeom->getGeometryType()) != wkbMultiLineString )
170 : {
171 1 : if (peErr != NULL)
172 1 : *peErr = OGRERR_FAILURE;
173 : CPLError(CE_Failure, CPLE_NotSupported,
174 : "The passed geometry is not an OGRGeometryCollection (or OGRMultiLineString) "
175 1 : "containing line string geometries");
176 1 : return NULL;
177 : }
178 :
179 3 : int bSuccess = TRUE;
180 3 : OGRGeometryCollection *poLines = (OGRGeometryCollection *) hLines;
181 3 : OGRPolygon *poPolygon = new OGRPolygon();
182 :
183 : (void) bBestEffort;
184 :
185 : /* -------------------------------------------------------------------- */
186 : /* Setup array of line markers indicating if they have been */
187 : /* added to a ring yet. */
188 : /* -------------------------------------------------------------------- */
189 3 : int nEdges = poLines->getNumGeometries();
190 3 : int *panEdgeConsumed, nRemainingEdges = nEdges;
191 :
192 3 : panEdgeConsumed = (int *) CPLCalloc(sizeof(int),nEdges);
193 :
194 : /* ==================================================================== */
195 : /* Loop generating rings. */
196 : /* ==================================================================== */
197 9 : while( nRemainingEdges > 0 )
198 : {
199 : int iEdge;
200 : OGRLineString *poLine;
201 :
202 : /* -------------------------------------------------------------------- */
203 : /* Find the first unconsumed edge. */
204 : /* -------------------------------------------------------------------- */
205 3 : for( iEdge = 0; panEdgeConsumed[iEdge]; iEdge++ ) {}
206 :
207 3 : poLine = (OGRLineString *) poLines->getGeometryRef(iEdge);
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Start a new ring, copying in the current line directly */
211 : /* -------------------------------------------------------------------- */
212 3 : OGRLinearRing *poRing = new OGRLinearRing();
213 :
214 3 : AddEdgeToRing( poRing, poLine, FALSE );
215 :
216 3 : panEdgeConsumed[iEdge] = TRUE;
217 3 : nRemainingEdges--;
218 :
219 : /* ==================================================================== */
220 : /* Loop adding edges to this ring until we make a whole pass */
221 : /* within finding anything to add. */
222 : /* ==================================================================== */
223 3 : int bWorkDone = TRUE;
224 3 : double dfBestDist = dfTolerance;
225 :
226 28 : while( !CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,NULL)
227 : && nRemainingEdges > 0
228 : && bWorkDone )
229 : {
230 22 : int iBestEdge = -1, bReverse = FALSE;
231 :
232 22 : bWorkDone = FALSE;
233 22 : dfBestDist = dfTolerance;
234 :
235 : // We consider linking the end to the beginning. If this is
236 : // closer than any other option we will just close the loop.
237 :
238 : //CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,&dfBestDist);
239 :
240 : // Find unused edge with end point closest to our loose end.
241 130 : for( iEdge = 0; iEdge < nEdges; iEdge++ )
242 : {
243 130 : if( panEdgeConsumed[iEdge] )
244 30 : continue;
245 :
246 100 : poLine = (OGRLineString *) poLines->getGeometryRef(iEdge);
247 :
248 100 : if( CheckPoints(poLine,0,poRing,poRing->getNumPoints()-1,
249 : &dfBestDist) )
250 : {
251 16 : iBestEdge = iEdge;
252 16 : bReverse = FALSE;
253 : }
254 100 : if( CheckPoints(poLine,poLine->getNumPoints()-1,
255 : poRing,poRing->getNumPoints()-1,
256 : &dfBestDist) )
257 : {
258 6 : iBestEdge = iEdge;
259 6 : bReverse = TRUE;
260 : }
261 :
262 : // if we use exact comparison, jump now
263 100 : if (dfTolerance == 0.0 && iBestEdge != -1) break;
264 : }
265 :
266 : // We found one within tolerance - add it.
267 22 : if( iBestEdge != -1 )
268 : {
269 : poLine = (OGRLineString *)
270 22 : poLines->getGeometryRef(iBestEdge);
271 :
272 22 : AddEdgeToRing( poRing, poLine, bReverse );
273 :
274 22 : panEdgeConsumed[iBestEdge] = TRUE;
275 22 : nRemainingEdges--;
276 22 : bWorkDone = TRUE;
277 : }
278 : }
279 :
280 : /* -------------------------------------------------------------------- */
281 : /* Did we fail to complete the ring? */
282 : /* -------------------------------------------------------------------- */
283 3 : dfBestDist = dfTolerance;
284 :
285 3 : if( !CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,
286 : &dfBestDist) )
287 : {
288 : CPLDebug( "OGR",
289 : "Failed to close ring %d.\n"
290 : "End Points are: (%.8f,%.7f) and (%.7f,%.7f)\n",
291 : poPolygon->getNumInteriorRings()+1,
292 : poRing->getX(0), poRing->getY(0),
293 : poRing->getX(poRing->getNumPoints()-1),
294 0 : poRing->getY(poRing->getNumPoints()-1) );
295 :
296 0 : bSuccess = FALSE;
297 : }
298 :
299 : /* -------------------------------------------------------------------- */
300 : /* Do we need to auto-close this ring? */
301 : /* -------------------------------------------------------------------- */
302 3 : if( bAutoClose
303 : && !CheckPoints(poRing,0,poRing,poRing->getNumPoints()-1,NULL) )
304 : {
305 : poRing->addPoint( poRing->getX(0),
306 : poRing->getY(0),
307 0 : poRing->getZ(0));
308 : }
309 :
310 3 : poPolygon->addRingDirectly( poRing );
311 : } /* next ring */
312 :
313 : /* -------------------------------------------------------------------- */
314 : /* Cleanup. */
315 : /* -------------------------------------------------------------------- */
316 3 : CPLFree( panEdgeConsumed );
317 :
318 : /* -------------------------------------------------------------------- */
319 : /* Identify exterior ring - it will be the largest. #3610 */
320 : /* -------------------------------------------------------------------- */
321 : double maxarea, tarea;
322 3 : int maxring = -1, rn, rcount;
323 3 : OGREnvelope tenv;
324 : OGRLinearRing *tring;
325 :
326 3 : tring = poPolygon->getExteriorRing();
327 3 : if (tring) tring->getEnvelope(&tenv);
328 3 : maxarea = (tenv.MaxX - tenv.MinX) * (tenv.MaxY - tenv.MinY);
329 :
330 3 : rcount = poPolygon->getNumInteriorRings();
331 3 : for (rn = 0; rn < rcount; ++rn) {
332 0 : tring = poPolygon->getInteriorRing(rn);
333 0 : tring->getEnvelope(&tenv);
334 0 : tarea = (tenv.MaxX - tenv.MinX) * (tenv.MaxY - tenv.MinY);
335 0 : if (tarea > maxarea) {
336 0 : maxarea = tarea;
337 0 : maxring = rn;
338 : }
339 : }
340 :
341 3 : if (maxring != -1) {
342 0 : OGRPolygon *poNewPoly = new OGRPolygon();
343 :
344 0 : poNewPoly->addRing(poPolygon->getInteriorRing(maxring));
345 0 : poNewPoly->addRing(poPolygon->getExteriorRing());
346 0 : for (rn = 0; rn < rcount; ++rn) {
347 0 : if (rn == maxring) continue;
348 0 : poNewPoly->addRing(poPolygon->getInteriorRing(rn));
349 : }
350 :
351 0 : delete poPolygon;
352 0 : poPolygon = poNewPoly;
353 : }
354 :
355 3 : if( peErr != NULL )
356 : {
357 3 : if( bSuccess )
358 3 : *peErr = OGRERR_NONE;
359 : else
360 0 : *peErr = OGRERR_FAILURE;
361 : }
362 :
363 3 : return (OGRGeometryH) poPolygon;
364 : }
365 :
366 :
367 :
|
