1 : /**********************************************************************
2 : * $Id: mitab_geometry.cpp,v 1.5 2004/06/30 20:29:04 dmorissette Exp $
3 : *
4 : * Name: mitab_geometry.cpp
5 : * Project: MapInfo TAB Read/Write library
6 : * Language: C++
7 : * Purpose: Geometry manipulation functions.
8 : * Author: Daniel Morissette, dmorissette@dmsolutions.ca
9 : * Based on functions from mapprimitive.c/mapsearch.c in the source
10 : * of UMN MapServer by Stephen Lime (http://mapserver.gis.umn.edu/)
11 : **********************************************************************
12 : * Copyright (c) 1999-2001, Daniel Morissette
13 : *
14 : * Permission is hereby granted, free of charge, to any person obtaining a
15 : * copy of this software and associated documentation files (the "Software"),
16 : * to deal in the Software without restriction, including without limitation
17 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 : * and/or sell copies of the Software, and to permit persons to whom the
19 : * Software is furnished to do so, subject to the following conditions:
20 : *
21 : * The above copyright notice and this permission notice shall be included
22 : * in all copies or substantial portions of the Software.
23 : *
24 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
25 : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 : * DEALINGS IN THE SOFTWARE.
31 : **********************************************************************
32 : *
33 : * $Log: mitab_geometry.cpp,v $
34 : * Revision 1.5 2004/06/30 20:29:04 dmorissette
35 : * Fixed refs to old address danmo@videotron.ca
36 : *
37 : * Revision 1.4 2001/12/18 23:42:28 daniel
38 : * Added a test in OGRPolygonLabelPoint() to prevent returning a point
39 : * outside of the polygon MBR (bug 673).
40 : *
41 : * Revision 1.3 2001/01/22 16:03:58 warmerda
42 : * expanded tabs
43 : *
44 : * Revision 1.2 2000/09/28 16:39:44 warmerda
45 : * avoid warnings for unused, and unitialized variables
46 : *
47 : * Revision 1.1 2000/09/19 17:19:40 daniel
48 : * Initial Revision
49 : *
50 : **********************************************************************/
51 :
52 : #include "ogr_geometry.h"
53 :
54 : #define OGR_NUM_RINGS(poly) (poly->getNumInteriorRings()+1)
55 : #define OGR_GET_RING(poly, i) (i==0?poly->getExteriorRing():poly->getInteriorRing(i-1))
56 :
57 :
58 : /**********************************************************************
59 : * OGRPointInRing()
60 : *
61 : * Returns TRUE is point is inside ring, FALSE otherwise
62 : *
63 : * Adapted version of msPointInPolygon() from MapServer's mapsearch.c
64 : **********************************************************************/
65 11 : GBool OGRPointInRing(OGRPoint *poPoint, OGRLineString *poRing)
66 : {
67 : int i, j, numpoints;
68 11 : GBool status = FALSE;
69 : double x, y;
70 :
71 11 : numpoints = poRing->getNumPoints();
72 11 : x = poPoint->getX();
73 11 : y = poPoint->getY();
74 :
75 261 : for (i = 0, j = numpoints-1; i < numpoints; j = i++)
76 : {
77 250 : if ((((poRing->getY(i)<=y) && (y<poRing->getY(j))) ||
78 : ((poRing->getY(j)<=y) && (y<poRing->getY(i)))) &&
79 : (x < (poRing->getX(j) - poRing->getX(i)) * (y - poRing->getY(i)) /
80 : (poRing->getY(j) - poRing->getY(i)) + poRing->getX(i)))
81 13 : status = !status;
82 : }
83 :
84 11 : return status;
85 : }
86 :
87 : /**********************************************************************
88 : * OGRIntersectPointPolygon()
89 : *
90 : * Instead of using ring orientation we count the number of parts the
91 : * point falls in. If odd the point is in the polygon, if 0 or even
92 : * then the point is in a hole or completely outside.
93 : *
94 : * Returns TRUE is point is inside polygon, FALSE otherwise
95 : *
96 : * Adapted version of msIntersectPointPolygon() from MapServer's mapsearch.c
97 : **********************************************************************/
98 11 : GBool OGRIntersectPointPolygon(OGRPoint *poPoint, OGRPolygon *poPoly)
99 : {
100 : int i;
101 11 : GBool status = FALSE;
102 :
103 22 : for(i=0; i<OGR_NUM_RINGS(poPoly); i++)
104 : {
105 11 : if (OGRPointInRing(poPoint, OGR_GET_RING(poPoly, i)))
106 : {
107 : /* ok, the point is in a line */
108 9 : status = !status;
109 : }
110 : }
111 :
112 11 : return status;
113 : }
114 :
115 :
116 : /**********************************************************************
117 : * OGRPolygonLabelPoint()
118 : *
119 : * Generate a label point on the surface of a polygon.
120 : *
121 : * The function is based on a scanline conversion routine used for polygon
122 : * fills. Instead of processing each line the as with drawing, the
123 : * polygon is sampled. The center of the longest sample is chosen for the
124 : * label point. The label point is guaranteed to be in the polygon even if
125 : * it has holes assuming the polygon is properly formed.
126 : *
127 : * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
128 : *
129 : * Adapted version of msPolygonLabelPoint() from MapServer's mapprimitive.c
130 : **********************************************************************/
131 :
132 : typedef enum {CLIP_LEFT, CLIP_MIDDLE, CLIP_RIGHT} CLIP_STATE;
133 : #define CLIP_CHECK(min, a, max) ((a) < (min) ? CLIP_LEFT : ((a) > (max) ? CLIP_RIGHT : CLIP_MIDDLE));
134 : #define ROUND(a) ( (a) + 0.5 )
135 : #define SWAP( a, b, t) ( (t) = (a), (a) = (b), (b) = (t) )
136 : #define EDGE_CHECK( x0, x, x1) ((x) < MIN( (x0), (x1)) ? CLIP_LEFT : ((x) > MAX( (x0), (x1)) ? CLIP_RIGHT : CLIP_MIDDLE ))
137 :
138 : #define NUM_SCANLINES 5
139 :
140 11 : int OGRPolygonLabelPoint(OGRPolygon *poPoly, OGRPoint *poLabelPoint)
141 : {
142 : double slope;
143 11 : OGRRawPoint point1, point2;
144 : int i, j, k, nfound;
145 : double x, y, *xintersect, temp;
146 : double hi_y, lo_y;
147 : int wrong_order, n;
148 11 : double len, max_len=0;
149 : double skip;
150 11 : OGREnvelope oEnv;
151 :
152 11 : if (poPoly == NULL)
153 0 : return OGRERR_FAILURE;
154 :
155 11 : poPoly->getEnvelope(&oEnv);
156 :
157 11 : poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
158 11 : poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);
159 :
160 : //if(get_centroid(p, lp, &miny, &maxy) == -1) return(-1);
161 :
162 11 : if(OGRIntersectPointPolygon(poLabelPoint, poPoly) == TRUE) /* cool, done */
163 9 : return OGRERR_NONE;
164 :
165 : /* do it the hard way - scanline */
166 :
167 2 : skip = (oEnv.MaxY - oEnv.MinY)/NUM_SCANLINES;
168 :
169 2 : n=0;
170 4 : for(j=0; j<OGR_NUM_RINGS(poPoly); j++)
171 : {
172 : /* count total number of points */
173 2 : n += OGR_GET_RING(poPoly, j)->getNumPoints();
174 : }
175 :
176 2 : xintersect = (double *)calloc(n, sizeof(double));
177 2 : if (xintersect == NULL)
178 0 : return OGRERR_FAILURE;
179 :
180 12 : for(k=1; k<=NUM_SCANLINES; k++)
181 : {
182 : /* sample the shape in the y direction */
183 :
184 10 : y = oEnv.MaxY - k*skip;
185 :
186 : /* need to find a y that won't intersect any vertices exactly */
187 10 : hi_y = y - 1; /* first initializing lo_y, hi_y to be any 2 pnts on either side of lp->y */
188 10 : lo_y = y + 1;
189 20 : for(j=0; j<OGR_NUM_RINGS(poPoly); j++)
190 : {
191 10 : OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);
192 :
193 10 : if((lo_y < y) && (hi_y >= y))
194 0 : break; /* already initialized */
195 128 : for(i=0; i < poRing->getNumPoints(); i++)
196 : {
197 126 : if((lo_y < y) && (hi_y >= y))
198 8 : break; /* already initialized */
199 118 : if(poRing->getY(i) < y)
200 12 : lo_y = poRing->getY(i);
201 118 : if(poRing->getY(i) >= y)
202 106 : hi_y = poRing->getY(i);
203 : }
204 : }
205 :
206 10 : n=0;
207 20 : for(j=0; j<OGR_NUM_RINGS(poPoly); j++)
208 : {
209 10 : OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);
210 :
211 255 : for(i=0; i < poRing->getNumPoints(); i++)
212 : {
213 245 : if((poRing->getY(i) < y) &&
214 : ((y - poRing->getY(i)) < (y - lo_y)))
215 5 : lo_y = poRing->getY(i);
216 245 : if((poRing->getY(i) >= y) &&
217 : ((poRing->getY(i) - y) < (hi_y - y)))
218 22 : hi_y = poRing->getY(i);
219 : }
220 : }
221 :
222 10 : if(lo_y == hi_y)
223 0 : return OGRERR_FAILURE;
224 : else
225 10 : y = (hi_y + lo_y)/2.0;
226 :
227 10 : nfound = 0;
228 20 : for(j=0; j<OGR_NUM_RINGS(poPoly); j++) /* for each line */
229 : {
230 10 : OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);
231 10 : point1.x = poRing->getX(poRing->getNumPoints()-1);
232 10 : point1.y = poRing->getY(poRing->getNumPoints()-1);
233 :
234 255 : for(i=0; i < poRing->getNumPoints(); i++)
235 : {
236 245 : point2.x = poRing->getX(i);
237 245 : point2.y = poRing->getY(i);
238 :
239 245 : if(EDGE_CHECK(point1.y, y, point2.y) == CLIP_MIDDLE)
240 : {
241 22 : if(point1.y == point2.y)
242 0 : continue; /* ignore horizontal edges */
243 : else
244 22 : slope = (point2.x - point1.x) / (point2.y - point1.y);
245 :
246 22 : x = point1.x + (y - point1.y)*slope;
247 22 : xintersect[nfound++] = x;
248 : } /* End of checking this edge */
249 :
250 245 : point1 = point2; /* Go on to next edge */
251 : }
252 : } /* Finished the scanline */
253 :
254 : /* First, sort the intersections */
255 19 : do
256 : {
257 19 : wrong_order = 0;
258 44 : for(i=0; i < nfound-1; i++)
259 : {
260 25 : if(xintersect[i] > xintersect[i+1])
261 : {
262 10 : wrong_order = 1;
263 10 : SWAP(xintersect[i], xintersect[i+1], temp);
264 : }
265 : }
266 : } while(wrong_order);
267 :
268 : /* Great, now find longest span */
269 10 : point1.y = point2.y = y;
270 21 : for(i=0; i < nfound; i += 2)
271 : {
272 11 : point1.x = xintersect[i];
273 11 : point2.x = xintersect[i+1];
274 : /* len = length(point1, point2); */
275 11 : len = ABS((point2.x - point1.x));
276 11 : if(len > max_len)
277 : {
278 4 : max_len = len;
279 4 : poLabelPoint->setX( (point1.x + point2.x)/2 );
280 4 : poLabelPoint->setY( y );
281 : }
282 : }
283 : }
284 :
285 2 : free(xintersect);
286 :
287 : /* __TODO__ Bug 673
288 : * There seem to be some polygons for which the label is returned
289 : * completely outside of the polygon's MBR and this messes the
290 : * file bounds, etc.
291 : * Until we find the source of the problem, we'll at least validate
292 : * the label point to make sure that it overlaps the polygon MBR.
293 : */
294 2 : if( poLabelPoint->getX() < oEnv.MinX
295 : || poLabelPoint->getY() < oEnv.MinY
296 : || poLabelPoint->getX() > oEnv.MaxX
297 : || poLabelPoint->getY() > oEnv.MaxY )
298 : {
299 : // Reset label coordinates to center of MBR, just in case
300 0 : poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
301 0 : poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);
302 :
303 : // And return an error
304 0 : return OGRERR_FAILURE;
305 : }
306 :
307 2 : if(max_len > 0)
308 2 : return OGRERR_NONE;
309 : else
310 0 : return OGRERR_FAILURE;
311 : }
312 :
313 :
314 : /**********************************************************************
315 : * OGRGetCentroid()
316 : *
317 : * Calculate polygon gravity center.
318 : *
319 : * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
320 : *
321 : * Adapted version of get_centroid() from MapServer's mapprimitive.c
322 : **********************************************************************/
323 :
324 0 : int OGRGetCentroid(OGRPolygon *poPoly, OGRPoint *poCentroid)
325 : {
326 : int i,j;
327 0 : double cent_weight_x=0.0, cent_weight_y=0.0;
328 0 : double len, total_len=0;
329 :
330 0 : for(i=0; i<OGR_NUM_RINGS(poPoly); i++)
331 : {
332 : double x1, y1, x2, y2;
333 0 : OGRLinearRing *poRing = OGR_GET_RING(poPoly, i);
334 :
335 0 : x2 = poRing->getX(0);
336 0 : y2 = poRing->getY(0);
337 :
338 0 : for(j=1; j<poRing->getNumPoints(); j++)
339 : {
340 0 : x1 = x2;
341 0 : y1 = y2;
342 0 : x2 = poRing->getX(j);
343 0 : y2 = poRing->getY(j);
344 :
345 0 : len = sqrt( pow((x2-x1),2) + pow((y2-y1),2) );
346 0 : cent_weight_x += len * ((x1 + x2)/2.0);
347 0 : cent_weight_y += len * ((y1 + y2)/2.0);
348 0 : total_len += len;
349 : }
350 : }
351 :
352 0 : if(total_len == 0)
353 0 : return(OGRERR_FAILURE);
354 :
355 0 : poCentroid->setX( cent_weight_x / total_len );
356 0 : poCentroid->setY( cent_weight_y / total_len );
357 :
358 0 : return OGRERR_NONE;
359 : }
360 :
361 :
362 :
363 : /**********************************************************************
364 : * OGRPolylineCenterPoint()
365 : *
366 : * Return the center point of a polyline.
367 : *
368 : * In MapInfo, for a simple or multiple polyline (pline), the center point
369 : * in the object definition is supposed to be either the center point of
370 : * the pline or the first section of a multiple pline (if an odd number of
371 : * points in the pline or first section), or the midway point between the
372 : * two central points (if an even number of points involved).
373 : *
374 : * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
375 : **********************************************************************/
376 :
377 0 : int OGRPolylineCenterPoint(OGRLineString *poLine, OGRPoint *poLabelPoint)
378 : {
379 0 : if (poLine == NULL || poLine->getNumPoints() < 2)
380 0 : return OGRERR_FAILURE;
381 :
382 0 : if (poLine->getNumPoints() % 2 == 0)
383 : {
384 : // Return the midway between the 2 center points
385 0 : int i = poLine->getNumPoints()/2;
386 0 : poLabelPoint->setX( (poLine->getX(i-1) + poLine->getX(i))/2.0 );
387 0 : poLabelPoint->setY( (poLine->getY(i-1) + poLine->getY(i))/2.0 );
388 : }
389 : else
390 : {
391 : // Return the center point
392 0 : poLine->getPoint(poLine->getNumPoints()/2, poLabelPoint);
393 : }
394 :
395 0 : return OGRERR_NONE;
396 : }
397 :
398 :
399 : /**********************************************************************
400 : * OGRPolylineLabelPoint()
401 : *
402 : * Generate a label point on a polyline: The center of the longest segment.
403 : *
404 : * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
405 : **********************************************************************/
406 :
407 0 : int OGRPolylineLabelPoint(OGRLineString *poLine, OGRPoint *poLabelPoint)
408 : {
409 0 : double segment_length, max_segment_length = 0.0;
410 : double x1, y1, x2, y2;
411 :
412 0 : if (poLine == NULL || poLine->getNumPoints() < 2)
413 0 : return OGRERR_FAILURE;
414 :
415 0 : max_segment_length = -1.0;
416 :
417 0 : x2 = poLine->getX(0);
418 0 : y2 = poLine->getY(0);
419 :
420 0 : for(int i=1; i<poLine->getNumPoints(); i++)
421 : {
422 0 : x1 = x2;
423 0 : y1 = y2;
424 0 : x2 = poLine->getX(i);
425 0 : y2 = poLine->getY(i);
426 :
427 0 : segment_length = pow((x2-x1),2) + pow((y2-y1),2);
428 0 : if (segment_length > max_segment_length)
429 : {
430 0 : max_segment_length = segment_length;
431 0 : poLabelPoint->setX( (x1 + x2)/2.0 );
432 0 : poLabelPoint->setY( (y1 + y2)/2.0 );
433 : }
434 : }
435 :
436 0 : return OGRERR_NONE;
437 : }
438 :
|
