source: josm/trunk/src/org/openstreetmap/josm/tools/Geometry.java@ 7318

Last change on this file since 7318 was 7193, checked in by bastiK, 10 years ago

added right- and left-hand traffic database
(new mapcss pseudo-class-condition :righthandtraffic)

Data file (data/left-right-hand-traffic.osm) license: CC0.
Loosely based on boundary lines data from naturalearthdata.com.
Traffic law info from wikipedia.
Boundaries generalized, especially in areas where there seem to be no roads.
Should be mostly correct, but some remote islands and atolls may be assigned wrongly.
Alway check before you start driving. :)

Added fast lookup index similar to QuadBuckets.
  • Property svn:eol-style set to native
File size: 35.8 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.tools;
3
4import java.awt.Rectangle;
5import java.awt.geom.Area;
6import java.awt.geom.Line2D;
7import java.awt.geom.Path2D;
8import java.math.BigDecimal;
9import java.math.MathContext;
10import java.util.ArrayList;
11import java.util.Collections;
12import java.util.Comparator;
13import java.util.EnumSet;
14import java.util.HashSet;
15import java.util.LinkedHashSet;
16import java.util.List;
17import java.util.Set;
18
19import org.openstreetmap.josm.Main;
20import org.openstreetmap.josm.command.AddCommand;
21import org.openstreetmap.josm.command.ChangeCommand;
22import org.openstreetmap.josm.command.Command;
23import org.openstreetmap.josm.data.coor.EastNorth;
24import org.openstreetmap.josm.data.coor.LatLon;
25import org.openstreetmap.josm.data.osm.BBox;
26import org.openstreetmap.josm.data.osm.MultipolygonCreate;
27import org.openstreetmap.josm.data.osm.Node;
28import org.openstreetmap.josm.data.osm.NodePositionComparator;
29import org.openstreetmap.josm.data.osm.OsmPrimitiveType;
30import org.openstreetmap.josm.data.osm.Relation;
31import org.openstreetmap.josm.data.osm.RelationMember;
32import org.openstreetmap.josm.data.osm.Way;
33
34/**
35 * Some tools for geometry related tasks.
36 *
37 * @author viesturs
38 */
39public final class Geometry {
40
41 private Geometry() {
42 // Hide default constructor for utils classes
43 }
44
45 public enum PolygonIntersection {FIRST_INSIDE_SECOND, SECOND_INSIDE_FIRST, OUTSIDE, CROSSING}
46
47 /**
48 * Will find all intersection and add nodes there for list of given ways.
49 * Handles self-intersections too.
50 * And makes commands to add the intersection points to ways.
51 *
52 * Prerequisite: no two nodes have the same coordinates.
53 *
54 * @param ways a list of ways to test
55 * @param test if false, do not build list of Commands, just return nodes
56 * @param cmds list of commands, typically empty when handed to this method.
57 * Will be filled with commands that add intersection nodes to
58 * the ways.
59 * @return list of new nodes
60 */
61 public static Set<Node> addIntersections(List<Way> ways, boolean test, List<Command> cmds) {
62
63 int n = ways.size();
64 @SuppressWarnings("unchecked")
65 List<Node>[] newNodes = new ArrayList[n];
66 BBox[] wayBounds = new BBox[n];
67 boolean[] changedWays = new boolean[n];
68
69 Set<Node> intersectionNodes = new LinkedHashSet<>();
70
71 //copy node arrays for local usage.
72 for (int pos = 0; pos < n; pos ++) {
73 newNodes[pos] = new ArrayList<>(ways.get(pos).getNodes());
74 wayBounds[pos] = getNodesBounds(newNodes[pos]);
75 changedWays[pos] = false;
76 }
77
78 //iterate over all way pairs and introduce the intersections
79 Comparator<Node> coordsComparator = new NodePositionComparator();
80 for (int seg1Way = 0; seg1Way < n; seg1Way ++) {
81 for (int seg2Way = seg1Way; seg2Way < n; seg2Way ++) {
82
83 //do not waste time on bounds that do not intersect
84 if (!wayBounds[seg1Way].intersects(wayBounds[seg2Way])) {
85 continue;
86 }
87
88 List<Node> way1Nodes = newNodes[seg1Way];
89 List<Node> way2Nodes = newNodes[seg2Way];
90
91 //iterate over primary segmemt
92 for (int seg1Pos = 0; seg1Pos + 1 < way1Nodes.size(); seg1Pos ++) {
93
94 //iterate over secondary segment
95 int seg2Start = seg1Way != seg2Way ? 0: seg1Pos + 2;//skip the adjacent segment
96
97 for (int seg2Pos = seg2Start; seg2Pos + 1< way2Nodes.size(); seg2Pos ++) {
98
99 //need to get them again every time, because other segments may be changed
100 Node seg1Node1 = way1Nodes.get(seg1Pos);
101 Node seg1Node2 = way1Nodes.get(seg1Pos + 1);
102 Node seg2Node1 = way2Nodes.get(seg2Pos);
103 Node seg2Node2 = way2Nodes.get(seg2Pos + 1);
104
105 int commonCount = 0;
106 //test if we have common nodes to add.
107 if (seg1Node1 == seg2Node1 || seg1Node1 == seg2Node2) {
108 commonCount ++;
109
110 if (seg1Way == seg2Way &&
111 seg1Pos == 0 &&
112 seg2Pos == way2Nodes.size() -2) {
113 //do not add - this is first and last segment of the same way.
114 } else {
115 intersectionNodes.add(seg1Node1);
116 }
117 }
118
119 if (seg1Node2 == seg2Node1 || seg1Node2 == seg2Node2) {
120 commonCount ++;
121
122 intersectionNodes.add(seg1Node2);
123 }
124
125 //no common nodes - find intersection
126 if (commonCount == 0) {
127 EastNorth intersection = getSegmentSegmentIntersection(
128 seg1Node1.getEastNorth(), seg1Node2.getEastNorth(),
129 seg2Node1.getEastNorth(), seg2Node2.getEastNorth());
130
131 if (intersection != null) {
132 if (test) {
133 intersectionNodes.add(seg2Node1);
134 return intersectionNodes;
135 }
136
137 Node newNode = new Node(Main.getProjection().eastNorth2latlon(intersection));
138 Node intNode = newNode;
139 boolean insertInSeg1 = false;
140 boolean insertInSeg2 = false;
141 //find if the intersection point is at end point of one of the segments, if so use that point
142
143 //segment 1
144 if (coordsComparator.compare(newNode, seg1Node1) == 0) {
145 intNode = seg1Node1;
146 } else if (coordsComparator.compare(newNode, seg1Node2) == 0) {
147 intNode = seg1Node2;
148 } else {
149 insertInSeg1 = true;
150 }
151
152 //segment 2
153 if (coordsComparator.compare(newNode, seg2Node1) == 0) {
154 intNode = seg2Node1;
155 } else if (coordsComparator.compare(newNode, seg2Node2) == 0) {
156 intNode = seg2Node2;
157 } else {
158 insertInSeg2 = true;
159 }
160
161 if (insertInSeg1) {
162 way1Nodes.add(seg1Pos +1, intNode);
163 changedWays[seg1Way] = true;
164
165 //fix seg2 position, as indexes have changed, seg2Pos is always bigger than seg1Pos on the same segment.
166 if (seg2Way == seg1Way) {
167 seg2Pos ++;
168 }
169 }
170
171 if (insertInSeg2) {
172 way2Nodes.add(seg2Pos +1, intNode);
173 changedWays[seg2Way] = true;
174
175 //Do not need to compare again to already split segment
176 seg2Pos ++;
177 }
178
179 intersectionNodes.add(intNode);
180
181 if (intNode == newNode) {
182 cmds.add(new AddCommand(intNode));
183 }
184 }
185 }
186 else if (test && !intersectionNodes.isEmpty())
187 return intersectionNodes;
188 }
189 }
190 }
191 }
192
193
194 for (int pos = 0; pos < ways.size(); pos ++) {
195 if (!changedWays[pos]) {
196 continue;
197 }
198
199 Way way = ways.get(pos);
200 Way newWay = new Way(way);
201 newWay.setNodes(newNodes[pos]);
202
203 cmds.add(new ChangeCommand(way, newWay));
204 }
205
206 return intersectionNodes;
207 }
208
209 private static BBox getNodesBounds(List<Node> nodes) {
210
211 BBox bounds = new BBox(nodes.get(0));
212 for (Node n: nodes) {
213 bounds.add(n.getCoor());
214 }
215 return bounds;
216 }
217
218 /**
219 * Tests if given point is to the right side of path consisting of 3 points.
220 *
221 * (Imagine the path is continued beyond the endpoints, so you get two rays
222 * starting from lineP2 and going through lineP1 and lineP3 respectively
223 * which divide the plane into two parts. The test returns true, if testPoint
224 * lies in the part that is to the right when traveling in the direction
225 * lineP1, lineP2, lineP3.)
226 *
227 * @param lineP1 first point in path
228 * @param lineP2 second point in path
229 * @param lineP3 third point in path
230 * @param testPoint
231 * @return true if to the right side, false otherwise
232 */
233 public static boolean isToTheRightSideOfLine(Node lineP1, Node lineP2, Node lineP3, Node testPoint) {
234 boolean pathBendToRight = angleIsClockwise(lineP1, lineP2, lineP3);
235 boolean rightOfSeg1 = angleIsClockwise(lineP1, lineP2, testPoint);
236 boolean rightOfSeg2 = angleIsClockwise(lineP2, lineP3, testPoint);
237
238 if (pathBendToRight)
239 return rightOfSeg1 && rightOfSeg2;
240 else
241 return !(!rightOfSeg1 && !rightOfSeg2);
242 }
243
244 /**
245 * This method tests if secondNode is clockwise to first node.
246 * @param commonNode starting point for both vectors
247 * @param firstNode first vector end node
248 * @param secondNode second vector end node
249 * @return true if first vector is clockwise before second vector.
250 */
251 public static boolean angleIsClockwise(Node commonNode, Node firstNode, Node secondNode) {
252 return angleIsClockwise(commonNode.getEastNorth(), firstNode.getEastNorth(), secondNode.getEastNorth());
253 }
254
255 /**
256 * Finds the intersection of two line segments
257 * @return EastNorth null if no intersection was found, the EastNorth coordinates of the intersection otherwise
258 */
259 public static EastNorth getSegmentSegmentIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) {
260
261 CheckParameterUtil.ensureValidCoordinates(p1, "p1");
262 CheckParameterUtil.ensureValidCoordinates(p2, "p2");
263 CheckParameterUtil.ensureValidCoordinates(p3, "p3");
264 CheckParameterUtil.ensureValidCoordinates(p4, "p4");
265
266 double x1 = p1.getX();
267 double y1 = p1.getY();
268 double x2 = p2.getX();
269 double y2 = p2.getY();
270 double x3 = p3.getX();
271 double y3 = p3.getY();
272 double x4 = p4.getX();
273 double y4 = p4.getY();
274
275 //TODO: do this locally.
276 //TODO: remove this check after careful testing
277 if (!Line2D.linesIntersect(x1, y1, x2, y2, x3, y3, x4, y4)) return null;
278
279 // solve line-line intersection in parametric form:
280 // (x1,y1) + (x2-x1,y2-y1)* u = (x3,y3) + (x4-x3,y4-y3)* v
281 // (x2-x1,y2-y1)*u - (x4-x3,y4-y3)*v = (x3-x1,y3-y1)
282 // if 0<= u,v <=1, intersection exists at ( x1+ (x2-x1)*u, y1 + (y2-y1)*u )
283
284 double a1 = x2 - x1;
285 double b1 = x3 - x4;
286 double c1 = x3 - x1;
287
288 double a2 = y2 - y1;
289 double b2 = y3 - y4;
290 double c2 = y3 - y1;
291
292 // Solve the equations
293 double det = a1*b2 - a2*b1;
294
295 double uu = b2*c1 - b1*c2 ;
296 double vv = a1*c2 - a2*c1;
297 double mag = Math.abs(uu)+Math.abs(vv);
298
299 if (Math.abs(det) > 1e-12 * mag) {
300 double u = uu/det, v = vv/det;
301 if (u>-1e-8 && u < 1+1e-8 && v>-1e-8 && v < 1+1e-8 ) {
302 if (u<0) u=0;
303 if (u>1) u=1.0;
304 return new EastNorth(x1+a1*u, y1+a2*u);
305 } else {
306 return null;
307 }
308 } else {
309 // parallel lines
310 return null;
311 }
312 }
313
314 /**
315 * Finds the intersection of two lines of infinite length.
316 * @return EastNorth null if no intersection was found, the coordinates of the intersection otherwise
317 * @throws IllegalArgumentException if a parameter is null or without valid coordinates
318 */
319 public static EastNorth getLineLineIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) {
320
321 CheckParameterUtil.ensureValidCoordinates(p1, "p1");
322 CheckParameterUtil.ensureValidCoordinates(p2, "p2");
323 CheckParameterUtil.ensureValidCoordinates(p3, "p3");
324 CheckParameterUtil.ensureValidCoordinates(p4, "p4");
325
326 if (!p1.isValid()) throw new IllegalArgumentException();
327
328 // Convert line from (point, point) form to ax+by=c
329 double a1 = p2.getY() - p1.getY();
330 double b1 = p1.getX() - p2.getX();
331 double c1 = p2.getX() * p1.getY() - p1.getX() * p2.getY();
332
333 double a2 = p4.getY() - p3.getY();
334 double b2 = p3.getX() - p4.getX();
335 double c2 = p4.getX() * p3.getY() - p3.getX() * p4.getY();
336
337 // Solve the equations
338 double det = a1 * b2 - a2 * b1;
339 if (det == 0)
340 return null; // Lines are parallel
341
342 return new EastNorth((b1 * c2 - b2 * c1) / det, (a2 * c1 - a1 * c2) / det);
343 }
344
345 public static boolean segmentsParallel(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) {
346
347 CheckParameterUtil.ensureValidCoordinates(p1, "p1");
348 CheckParameterUtil.ensureValidCoordinates(p2, "p2");
349 CheckParameterUtil.ensureValidCoordinates(p3, "p3");
350 CheckParameterUtil.ensureValidCoordinates(p4, "p4");
351
352 // Convert line from (point, point) form to ax+by=c
353 double a1 = p2.getY() - p1.getY();
354 double b1 = p1.getX() - p2.getX();
355
356 double a2 = p4.getY() - p3.getY();
357 double b2 = p3.getX() - p4.getX();
358
359 // Solve the equations
360 double det = a1 * b2 - a2 * b1;
361 // remove influence of of scaling factor
362 det /= Math.sqrt(a1*a1 + b1*b1) * Math.sqrt(a2*a2 + b2*b2);
363 return Math.abs(det) < 1e-3;
364 }
365
366 private static EastNorth closestPointTo(EastNorth p1, EastNorth p2, EastNorth point, boolean segmentOnly) {
367 CheckParameterUtil.ensureParameterNotNull(p1, "p1");
368 CheckParameterUtil.ensureParameterNotNull(p2, "p2");
369 CheckParameterUtil.ensureParameterNotNull(point, "point");
370
371 double ldx = p2.getX() - p1.getX();
372 double ldy = p2.getY() - p1.getY();
373
374 if (ldx == 0 && ldy == 0) //segment zero length
375 return p1;
376
377 double pdx = point.getX() - p1.getX();
378 double pdy = point.getY() - p1.getY();
379
380 double offset = (pdx * ldx + pdy * ldy) / (ldx * ldx + ldy * ldy);
381
382 if (segmentOnly && offset <= 0)
383 return p1;
384 else if (segmentOnly && offset >= 1)
385 return p2;
386 else
387 return new EastNorth(p1.getX() + ldx * offset, p1.getY() + ldy * offset);
388 }
389
390 /**
391 * Calculates closest point to a line segment.
392 * @param segmentP1 First point determining line segment
393 * @param segmentP2 Second point determining line segment
394 * @param point Point for which a closest point is searched on line segment [P1,P2]
395 * @return segmentP1 if it is the closest point, segmentP2 if it is the closest point,
396 * a new point if closest point is between segmentP1 and segmentP2.
397 * @since 3650
398 * @see #closestPointToLine
399 */
400 public static EastNorth closestPointToSegment(EastNorth segmentP1, EastNorth segmentP2, EastNorth point) {
401 return closestPointTo(segmentP1, segmentP2, point, true);
402 }
403
404 /**
405 * Calculates closest point to a line.
406 * @param lineP1 First point determining line
407 * @param lineP2 Second point determining line
408 * @param point Point for which a closest point is searched on line (P1,P2)
409 * @return The closest point found on line. It may be outside the segment [P1,P2].
410 * @since 4134
411 * @see #closestPointToSegment
412 */
413 public static EastNorth closestPointToLine(EastNorth lineP1, EastNorth lineP2, EastNorth point) {
414 return closestPointTo(lineP1, lineP2, point, false);
415 }
416
417 /**
418 * This method tests if secondNode is clockwise to first node.
419 *
420 * The line through the two points commonNode and firstNode divides the
421 * plane into two parts. The test returns true, if secondNode lies in
422 * the part that is to the right when traveling in the direction from
423 * commonNode to firstNode.
424 *
425 * @param commonNode starting point for both vectors
426 * @param firstNode first vector end node
427 * @param secondNode second vector end node
428 * @return true if first vector is clockwise before second vector.
429 */
430 public static boolean angleIsClockwise(EastNorth commonNode, EastNorth firstNode, EastNorth secondNode) {
431
432 CheckParameterUtil.ensureValidCoordinates(commonNode, "commonNode");
433 CheckParameterUtil.ensureValidCoordinates(firstNode, "firstNode");
434 CheckParameterUtil.ensureValidCoordinates(secondNode, "secondNode");
435
436 double dy1 = (firstNode.getY() - commonNode.getY());
437 double dy2 = (secondNode.getY() - commonNode.getY());
438 double dx1 = (firstNode.getX() - commonNode.getX());
439 double dx2 = (secondNode.getX() - commonNode.getX());
440
441 return dy1 * dx2 - dx1 * dy2 > 0;
442 }
443
444 /**
445 * Returns the Area of a polygon, from its list of nodes.
446 * @param polygon List of nodes forming polygon (EastNorth coordinates)
447 * @return Area for the given list of nodes
448 * @since 6841
449 */
450 public static Area getArea(List<Node> polygon) {
451 Path2D path = new Path2D.Double();
452
453 boolean begin = true;
454 for (Node n : polygon) {
455 EastNorth en = n.getEastNorth();
456 if (en != null) {
457 if (begin) {
458 path.moveTo(en.getX(), en.getY());
459 begin = false;
460 } else {
461 path.lineTo(en.getX(), en.getY());
462 }
463 }
464 }
465 if (!begin) {
466 path.closePath();
467 }
468
469 return new Area(path);
470 }
471
472 /**
473 * Returns the Area of a polygon, from its list of nodes.
474 * @param polygon List of nodes forming polygon (LatLon coordinates)
475 * @return Area for the given list of nodes
476 * @since 6841
477 */
478 public static Area getAreaLatLon(List<Node> polygon) {
479 Path2D path = new Path2D.Double();
480
481 boolean begin = true;
482 for (Node n : polygon) {
483 if (begin) {
484 path.moveTo(n.getCoor().lon(), n.getCoor().lat());
485 begin = false;
486 } else {
487 path.lineTo(n.getCoor().lon(), n.getCoor().lat());
488 }
489 }
490 if (!begin) {
491 path.closePath();
492 }
493
494 return new Area(path);
495 }
496
497 /**
498 * Tests if two polygons intersect.
499 * @param first List of nodes forming first polygon
500 * @param second List of nodes forming second polygon
501 * @return intersection kind
502 */
503 public static PolygonIntersection polygonIntersection(List<Node> first, List<Node> second) {
504 Area a1 = getArea(first);
505 Area a2 = getArea(second);
506 return polygonIntersection(a1, a2);
507 }
508
509 /**
510 * Tests if two polygons intersect.
511 * @param a1 Area of first polygon
512 * @param a2 Area of second polygon
513 * @return intersection kind
514 * @since 6841
515 */
516 public static PolygonIntersection polygonIntersection(Area a1, Area a2) {
517 return polygonIntersection(a1, a2, 1.0);
518 }
519
520 /**
521 * Tests if two polygons intersect.
522 * @param a1 Area of first polygon
523 * @param a2 Area of second polygon
524 * @param eps an area threshold, everything below is considered an empty intersection
525 * @return intersection kind
526 */
527 public static PolygonIntersection polygonIntersection(Area a1, Area a2, double eps) {
528
529 Area inter = new Area(a1);
530 inter.intersect(a2);
531
532 Rectangle bounds = inter.getBounds();
533
534 if (inter.isEmpty() || bounds.getHeight()*bounds.getWidth() <= eps) {
535 return PolygonIntersection.OUTSIDE;
536 } else if (inter.equals(a1)) {
537 return PolygonIntersection.FIRST_INSIDE_SECOND;
538 } else if (inter.equals(a2)) {
539 return PolygonIntersection.SECOND_INSIDE_FIRST;
540 } else {
541 return PolygonIntersection.CROSSING;
542 }
543 }
544
545 /**
546 * Tests if point is inside a polygon. The polygon can be self-intersecting. In such case the contains function works in xor-like manner.
547 * @param polygonNodes list of nodes from polygon path.
548 * @param point the point to test
549 * @return true if the point is inside polygon.
550 */
551 public static boolean nodeInsidePolygon(Node point, List<Node> polygonNodes) {
552 if (polygonNodes.size() < 2)
553 return false;
554
555 boolean inside = false;
556 Node p1, p2;
557
558 //iterate each side of the polygon, start with the last segment
559 Node oldPoint = polygonNodes.get(polygonNodes.size() - 1);
560
561 for (Node newPoint : polygonNodes) {
562 //skip duplicate points
563 if (newPoint.equals(oldPoint)) {
564 continue;
565 }
566
567 //order points so p1.lat <= p2.lat
568 if (newPoint.getEastNorth().getY() > oldPoint.getEastNorth().getY()) {
569 p1 = oldPoint;
570 p2 = newPoint;
571 } else {
572 p1 = newPoint;
573 p2 = oldPoint;
574 }
575
576 //test if the line is crossed and if so invert the inside flag.
577 if ((newPoint.getEastNorth().getY() < point.getEastNorth().getY()) == (point.getEastNorth().getY() <= oldPoint.getEastNorth().getY())
578 && (point.getEastNorth().getX() - p1.getEastNorth().getX()) * (p2.getEastNorth().getY() - p1.getEastNorth().getY())
579 < (p2.getEastNorth().getX() - p1.getEastNorth().getX()) * (point.getEastNorth().getY() - p1.getEastNorth().getY()))
580 {
581 inside = !inside;
582 }
583
584 oldPoint = newPoint;
585 }
586
587 return inside;
588 }
589
590 /**
591 * Returns area of a closed way in square meters.
592 * (approximate(?), but should be OK for small areas)
593 *
594 * Relies on the current projection: Works correctly, when
595 * one unit in projected coordinates corresponds to one meter.
596 * This is true for most projections, but not for WGS84 and
597 * Mercator (EPSG:3857).
598 *
599 * @param way Way to measure, should be closed (first node is the same as last node)
600 * @return area of the closed way.
601 */
602 public static double closedWayArea(Way way) {
603
604 //http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/
605 double area = 0;
606 Node lastN = null;
607 for (Node n : way.getNodes()) {
608 if (lastN != null) {
609 n.getEastNorth().getX();
610
611 area += (calcX(n) * calcY(lastN)) - (calcY(n) * calcX(lastN));
612 }
613 lastN = n;
614 }
615 return Math.abs(area/2);
616 }
617
618 protected static double calcX(Node p1){
619 double lat1, lon1, lat2, lon2;
620 double dlon, dlat;
621
622 lat1 = p1.getCoor().lat() * Math.PI / 180.0;
623 lon1 = p1.getCoor().lon() * Math.PI / 180.0;
624 lat2 = lat1;
625 lon2 = 0;
626
627 dlon = lon2 - lon1;
628 dlat = lat2 - lat1;
629
630 double a = (Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2));
631 double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
632 return 6367000 * c;
633 }
634
635 protected static double calcY(Node p1){
636 double lat1, lon1, lat2, lon2;
637 double dlon, dlat;
638
639 lat1 = p1.getCoor().lat() * Math.PI / 180.0;
640 lon1 = p1.getCoor().lon() * Math.PI / 180.0;
641 lat2 = 0;
642 lon2 = lon1;
643
644 dlon = lon2 - lon1;
645 dlat = lat2 - lat1;
646
647 double a = (Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2));
648 double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
649 return 6367000 * c;
650 }
651
652 /**
653 * Determines whether a way is oriented clockwise.
654 *
655 * Internals: Assuming a closed non-looping way, compute twice the area
656 * of the polygon using the formula {@code 2 * area = sum (X[n] * Y[n+1] - X[n+1] * Y[n])}.
657 * If the area is negative the way is ordered in a clockwise direction.
658 *
659 * See http://paulbourke.net/geometry/polyarea/
660 *
661 * @param w the way to be checked.
662 * @return true if and only if way is oriented clockwise.
663 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}).
664 */
665 public static boolean isClockwise(Way w) {
666 if (!w.isClosed()) {
667 throw new IllegalArgumentException("Way must be closed to check orientation.");
668 }
669
670 double area2 = 0.;
671 int nodesCount = w.getNodesCount();
672
673 for (int node = 1; node <= /*sic! consider last-first as well*/ nodesCount; node++) {
674 LatLon coorPrev = w.getNode(node - 1).getCoor();
675 LatLon coorCurr = w.getNode(node % nodesCount).getCoor();
676 area2 += coorPrev.lon() * coorCurr.lat();
677 area2 -= coorCurr.lon() * coorPrev.lat();
678 }
679 return area2 < 0;
680 }
681
682 /**
683 * Returns angle of a segment defined with 2 point coordinates.
684 *
685 * @param p1
686 * @param p2
687 * @return Angle in radians (-pi, pi]
688 */
689 public static double getSegmentAngle(EastNorth p1, EastNorth p2) {
690
691 CheckParameterUtil.ensureValidCoordinates(p1, "p1");
692 CheckParameterUtil.ensureValidCoordinates(p2, "p2");
693
694 return Math.atan2(p2.north() - p1.north(), p2.east() - p1.east());
695 }
696
697 /**
698 * Returns angle of a corner defined with 3 point coordinates.
699 *
700 * @param p1
701 * @param p2 Common endpoint
702 * @param p3
703 * @return Angle in radians (-pi, pi]
704 */
705 public static double getCornerAngle(EastNorth p1, EastNorth p2, EastNorth p3) {
706
707 CheckParameterUtil.ensureValidCoordinates(p1, "p1");
708 CheckParameterUtil.ensureValidCoordinates(p2, "p2");
709 CheckParameterUtil.ensureValidCoordinates(p3, "p3");
710
711 Double result = getSegmentAngle(p2, p1) - getSegmentAngle(p2, p3);
712 if (result <= -Math.PI) {
713 result += 2 * Math.PI;
714 }
715
716 if (result > Math.PI) {
717 result -= 2 * Math.PI;
718 }
719
720 return result;
721 }
722
723 /**
724 * Compute the centroid/barycenter of nodes
725 * @param nodes Nodes for which the centroid is wanted
726 * @return the centroid of nodes
727 * @see Geometry#getCenter
728 */
729 public static EastNorth getCentroid(List<Node> nodes) {
730
731 BigDecimal area = BigDecimal.ZERO;
732 BigDecimal north = BigDecimal.ZERO;
733 BigDecimal east = BigDecimal.ZERO;
734
735 // See https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon for the equation used here
736 for (int i = 0; i < nodes.size(); i++) {
737 EastNorth n0 = nodes.get(i).getEastNorth();
738 EastNorth n1 = nodes.get((i+1) % nodes.size()).getEastNorth();
739
740 if (n0 != null && n1 != null && n0.isValid() && n1.isValid()) {
741 BigDecimal x0 = new BigDecimal(n0.east());
742 BigDecimal y0 = new BigDecimal(n0.north());
743 BigDecimal x1 = new BigDecimal(n1.east());
744 BigDecimal y1 = new BigDecimal(n1.north());
745
746 BigDecimal k = x0.multiply(y1, MathContext.DECIMAL128).subtract(y0.multiply(x1, MathContext.DECIMAL128));
747
748 area = area.add(k, MathContext.DECIMAL128);
749 east = east.add(k.multiply(x0.add(x1, MathContext.DECIMAL128), MathContext.DECIMAL128));
750 north = north.add(k.multiply(y0.add(y1, MathContext.DECIMAL128), MathContext.DECIMAL128));
751 }
752 }
753
754 BigDecimal d = new BigDecimal(3, MathContext.DECIMAL128); // 1/2 * 6 = 3
755 area = area.multiply(d, MathContext.DECIMAL128);
756 if (area.compareTo(BigDecimal.ZERO) != 0) {
757 north = north.divide(area, MathContext.DECIMAL128);
758 east = east.divide(area, MathContext.DECIMAL128);
759 }
760
761 return new EastNorth(east.doubleValue(), north.doubleValue());
762 }
763
764 /**
765 * Compute center of the circle closest to different nodes.
766 *
767 * Ensure exact center computation in case nodes are already aligned in circle.
768 * This is done by least square method.
769 * Let be a_i x + b_i y + c_i = 0 equations of bisectors of each edges.
770 * Center must be intersection of all bisectors.
771 * <pre>
772 * [ a1 b1 ] [ -c1 ]
773 * With A = [ ... ... ] and Y = [ ... ]
774 * [ an bn ] [ -cn ]
775 * </pre>
776 * An approximation of center of circle is (At.A)^-1.At.Y
777 * @param nodes Nodes parts of the circle (at least 3)
778 * @return An approximation of the center, of null if there is no solution.
779 * @see Geometry#getCentroid
780 * @since 6934
781 */
782 public static EastNorth getCenter(List<Node> nodes) {
783 int nc = nodes.size();
784 if(nc < 3) return null;
785 /**
786 * Equation of each bisector ax + by + c = 0
787 */
788 double[] a = new double[nc];
789 double[] b = new double[nc];
790 double[] c = new double[nc];
791 // Compute equation of bisector
792 for(int i = 0; i < nc; i++) {
793 EastNorth pt1 = nodes.get(i).getEastNorth();
794 EastNorth pt2 = nodes.get((i+1) % nc).getEastNorth();
795 a[i] = pt1.east() - pt2.east();
796 b[i] = pt1.north() - pt2.north();
797 double d = Math.sqrt(a[i]*a[i] + b[i]*b[i]);
798 if(d == 0) return null;
799 a[i] /= d;
800 b[i] /= d;
801 double xC = (pt1.east() + pt2.east()) / 2;
802 double yC = (pt1.north() + pt2.north()) / 2;
803 c[i] = -(a[i]*xC + b[i]*yC);
804 }
805 // At.A = [aij]
806 double a11 = 0, a12 = 0, a22 = 0;
807 // At.Y = [bi]
808 double b1 = 0, b2 = 0;
809 for(int i = 0; i < nc; i++) {
810 a11 += a[i]*a[i];
811 a12 += a[i]*b[i];
812 a22 += b[i]*b[i];
813 b1 -= a[i]*c[i];
814 b2 -= b[i]*c[i];
815 }
816 // (At.A)^-1 = [invij]
817 double det = a11*a22 - a12*a12;
818 if(Math.abs(det) < 1e-5) return null;
819 double inv11 = a22/det;
820 double inv12 = -a12/det;
821 double inv22 = a11/det;
822 // center (xC, yC) = (At.A)^-1.At.y
823 double xC = inv11*b1 + inv12*b2;
824 double yC = inv12*b1 + inv22*b2;
825 return new EastNorth(xC, yC);
826 }
827
828 /**
829 * Returns the coordinate of intersection of segment sp1-sp2 and an altitude
830 * to it starting at point ap. If the line defined with sp1-sp2 intersects
831 * its altitude out of sp1-sp2, null is returned.
832 *
833 * @param sp1
834 * @param sp2
835 * @param ap
836 * @return Intersection coordinate or null
837 */
838 public static EastNorth getSegmentAltituteIntersection(EastNorth sp1, EastNorth sp2, EastNorth ap) {
839
840 CheckParameterUtil.ensureValidCoordinates(sp1, "sp1");
841 CheckParameterUtil.ensureValidCoordinates(sp2, "sp2");
842 CheckParameterUtil.ensureValidCoordinates(ap, "ap");
843
844 Double segmentLenght = sp1.distance(sp2);
845 Double altitudeAngle = getSegmentAngle(sp1, sp2) + Math.PI / 2;
846
847 // Taking a random point on the altitude line (angle is known).
848 EastNorth ap2 = new EastNorth(ap.east() + 1000
849 * Math.cos(altitudeAngle), ap.north() + 1000
850 * Math.sin(altitudeAngle));
851
852 // Finding the intersection of two lines
853 EastNorth resultCandidate = Geometry.getLineLineIntersection(sp1, sp2,
854 ap, ap2);
855
856 // Filtering result
857 if (resultCandidate != null
858 && resultCandidate.distance(sp1) * .999 < segmentLenght
859 && resultCandidate.distance(sp2) * .999 < segmentLenght) {
860 return resultCandidate;
861 } else {
862 return null;
863 }
864 }
865
866 public static class MultiPolygonMembers {
867 public final Set<Way> outers = new HashSet<>();
868 public final Set<Way> inners = new HashSet<>();
869
870 public MultiPolygonMembers(Relation multiPolygon) {
871 for (RelationMember m : multiPolygon.getMembers()) {
872 if (m.getType().equals(OsmPrimitiveType.WAY)) {
873 if ("outer".equals(m.getRole())) {
874 outers.add(m.getWay());
875 } else if ("inner".equals(m.getRole())) {
876 inners.add(m.getWay());
877 }
878 }
879 }
880 }
881 }
882
883 /**
884 * Tests if the {@code node} is inside the multipolygon {@code multiPolygon}. The nullable argument
885 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match.
886 */
887 public static boolean isNodeInsideMultiPolygon(Node node, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) {
888 return isPolygonInsideMultiPolygon(Collections.singletonList(node), multiPolygon, isOuterWayAMatch);
889 }
890
891 /**
892 * Tests if the polygon formed by {@code nodes} is inside the multipolygon {@code multiPolygon}. The nullable argument
893 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match.
894 * <p>
895 * If {@code nodes} contains exactly one element, then it is checked whether that one node is inside the multipolygon.
896 */
897 public static boolean isPolygonInsideMultiPolygon(List<Node> nodes, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) {
898 // Extract outer/inner members from multipolygon
899 final MultiPolygonMembers mpm = new MultiPolygonMembers(multiPolygon);
900 // Construct complete rings for the inner/outer members
901 final List<MultipolygonCreate.JoinedPolygon> outerRings;
902 final List<MultipolygonCreate.JoinedPolygon> innerRings;
903 try {
904 outerRings = MultipolygonCreate.joinWays(mpm.outers);
905 innerRings = MultipolygonCreate.joinWays(mpm.inners);
906 } catch (MultipolygonCreate.JoinedPolygonCreationException ex) {
907 Main.debug("Invalid multipolygon " + multiPolygon);
908 return false;
909 }
910 // Test if object is inside an outer member
911 for (MultipolygonCreate.JoinedPolygon out : outerRings) {
912 if (nodes.size() == 1
913 ? nodeInsidePolygon(nodes.get(0), out.getNodes())
914 : EnumSet.of(PolygonIntersection.FIRST_INSIDE_SECOND, PolygonIntersection.CROSSING).contains(polygonIntersection(nodes, out.getNodes()))) {
915 boolean insideInner = false;
916 // If inside an outer, check it is not inside an inner
917 for (MultipolygonCreate.JoinedPolygon in : innerRings) {
918 if (polygonIntersection(in.getNodes(), out.getNodes()) == PolygonIntersection.FIRST_INSIDE_SECOND
919 && (nodes.size() == 1
920 ? nodeInsidePolygon(nodes.get(0), in.getNodes())
921 : polygonIntersection(nodes, in.getNodes()) == PolygonIntersection.FIRST_INSIDE_SECOND)) {
922 insideInner = true;
923 break;
924 }
925 }
926 // Inside outer but not inside inner -> the polygon appears to be inside a the multipolygon
927 if (!insideInner) {
928 // Final check using predicate
929 if (isOuterWayAMatch == null || isOuterWayAMatch.evaluate(out.ways.get(0) /* TODO give a better representation of the outer ring to the predicate */)) {
930 return true;
931 }
932 }
933 }
934 }
935 return false;
936 }
937}
Note: See TracBrowser for help on using the repository browser.