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