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