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