001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.tools; 003 004import java.awt.Rectangle; 005import java.awt.geom.Area; 006import java.awt.geom.Line2D; 007import java.awt.geom.Path2D; 008import java.math.BigDecimal; 009import java.math.MathContext; 010import java.util.ArrayList; 011import java.util.Collections; 012import java.util.Comparator; 013import java.util.EnumSet; 014import java.util.LinkedHashSet; 015import java.util.List; 016import java.util.Set; 017import java.util.function.Predicate; 018 019import org.openstreetmap.josm.Main; 020import org.openstreetmap.josm.command.AddCommand; 021import org.openstreetmap.josm.command.ChangeCommand; 022import org.openstreetmap.josm.command.Command; 023import org.openstreetmap.josm.data.coor.EastNorth; 024import org.openstreetmap.josm.data.osm.BBox; 025import org.openstreetmap.josm.data.osm.DataSet; 026import org.openstreetmap.josm.data.osm.MultipolygonBuilder; 027import org.openstreetmap.josm.data.osm.MultipolygonBuilder.JoinedPolygon; 028import org.openstreetmap.josm.data.osm.Node; 029import org.openstreetmap.josm.data.osm.NodePositionComparator; 030import org.openstreetmap.josm.data.osm.OsmPrimitive; 031import org.openstreetmap.josm.data.osm.Relation; 032import org.openstreetmap.josm.data.osm.Way; 033import org.openstreetmap.josm.data.osm.visitor.paint.relations.Multipolygon; 034import org.openstreetmap.josm.data.osm.visitor.paint.relations.MultipolygonCache; 035import org.openstreetmap.josm.data.projection.Projection; 036import org.openstreetmap.josm.data.projection.Projections; 037 038/** 039 * Some tools for geometry related tasks. 040 * 041 * @author viesturs 042 */ 043public final class Geometry { 044 045 private Geometry() { 046 // Hide default constructor for utils classes 047 } 048 049 /** 050 * The result types for a {@link Geometry#polygonIntersection(Area, Area)} test 051 */ 052 public enum PolygonIntersection { 053 /** 054 * The first polygon is inside the second one 055 */ 056 FIRST_INSIDE_SECOND, 057 /** 058 * The second one is inside the first 059 */ 060 SECOND_INSIDE_FIRST, 061 /** 062 * The polygons do not overlap 063 */ 064 OUTSIDE, 065 /** 066 * The polygon borders cross each other 067 */ 068 CROSSING 069 } 070 071 /** 072 * Will find all intersection and add nodes there for list of given ways. 073 * Handles self-intersections too. 074 * And makes commands to add the intersection points to ways. 075 * 076 * Prerequisite: no two nodes have the same coordinates. 077 * 078 * @param ways a list of ways to test 079 * @param test if false, do not build list of Commands, just return nodes 080 * @param cmds list of commands, typically empty when handed to this method. 081 * Will be filled with commands that add intersection nodes to 082 * the ways. 083 * @return list of new nodes 084 */ 085 public static Set<Node> addIntersections(List<Way> ways, boolean test, List<Command> cmds) { 086 087 int n = ways.size(); 088 @SuppressWarnings("unchecked") 089 List<Node>[] newNodes = new ArrayList[n]; 090 BBox[] wayBounds = new BBox[n]; 091 boolean[] changedWays = new boolean[n]; 092 093 Set<Node> intersectionNodes = new LinkedHashSet<>(); 094 095 //copy node arrays for local usage. 096 for (int pos = 0; pos < n; pos++) { 097 newNodes[pos] = new ArrayList<>(ways.get(pos).getNodes()); 098 wayBounds[pos] = getNodesBounds(newNodes[pos]); 099 changedWays[pos] = false; 100 } 101 102 DataSet dataset = ways.get(0).getDataSet(); 103 104 //iterate over all way pairs and introduce the intersections 105 Comparator<Node> coordsComparator = new NodePositionComparator(); 106 for (int seg1Way = 0; seg1Way < n; seg1Way++) { 107 for (int seg2Way = seg1Way; seg2Way < n; seg2Way++) { 108 109 //do not waste time on bounds that do not intersect 110 if (!wayBounds[seg1Way].intersects(wayBounds[seg2Way])) { 111 continue; 112 } 113 114 List<Node> way1Nodes = newNodes[seg1Way]; 115 List<Node> way2Nodes = newNodes[seg2Way]; 116 117 //iterate over primary segmemt 118 for (int seg1Pos = 0; seg1Pos + 1 < way1Nodes.size(); seg1Pos++) { 119 120 //iterate over secondary segment 121 int seg2Start = seg1Way != seg2Way ? 0 : seg1Pos + 2; //skip the adjacent segment 122 123 for (int seg2Pos = seg2Start; seg2Pos + 1 < way2Nodes.size(); seg2Pos++) { 124 125 //need to get them again every time, because other segments may be changed 126 Node seg1Node1 = way1Nodes.get(seg1Pos); 127 Node seg1Node2 = way1Nodes.get(seg1Pos + 1); 128 Node seg2Node1 = way2Nodes.get(seg2Pos); 129 Node seg2Node2 = way2Nodes.get(seg2Pos + 1); 130 131 int commonCount = 0; 132 //test if we have common nodes to add. 133 if (seg1Node1 == seg2Node1 || seg1Node1 == seg2Node2) { 134 commonCount++; 135 136 if (seg1Way == seg2Way && 137 seg1Pos == 0 && 138 seg2Pos == way2Nodes.size() -2) { 139 //do not add - this is first and last segment of the same way. 140 } else { 141 intersectionNodes.add(seg1Node1); 142 } 143 } 144 145 if (seg1Node2 == seg2Node1 || seg1Node2 == seg2Node2) { 146 commonCount++; 147 148 intersectionNodes.add(seg1Node2); 149 } 150 151 //no common nodes - find intersection 152 if (commonCount == 0) { 153 EastNorth intersection = getSegmentSegmentIntersection( 154 seg1Node1.getEastNorth(), seg1Node2.getEastNorth(), 155 seg2Node1.getEastNorth(), seg2Node2.getEastNorth()); 156 157 if (intersection != null) { 158 if (test) { 159 intersectionNodes.add(seg2Node1); 160 return intersectionNodes; 161 } 162 163 Node newNode = new Node(Main.getProjection().eastNorth2latlon(intersection)); 164 Node intNode = newNode; 165 boolean insertInSeg1 = false; 166 boolean insertInSeg2 = false; 167 //find if the intersection point is at end point of one of the segments, if so use that point 168 169 //segment 1 170 if (coordsComparator.compare(newNode, seg1Node1) == 0) { 171 intNode = seg1Node1; 172 } else if (coordsComparator.compare(newNode, seg1Node2) == 0) { 173 intNode = seg1Node2; 174 } else { 175 insertInSeg1 = true; 176 } 177 178 //segment 2 179 if (coordsComparator.compare(newNode, seg2Node1) == 0) { 180 intNode = seg2Node1; 181 } else if (coordsComparator.compare(newNode, seg2Node2) == 0) { 182 intNode = seg2Node2; 183 } else { 184 insertInSeg2 = true; 185 } 186 187 if (insertInSeg1) { 188 way1Nodes.add(seg1Pos +1, intNode); 189 changedWays[seg1Way] = true; 190 191 //fix seg2 position, as indexes have changed, seg2Pos is always bigger than seg1Pos on the same segment. 192 if (seg2Way == seg1Way) { 193 seg2Pos++; 194 } 195 } 196 197 if (insertInSeg2) { 198 way2Nodes.add(seg2Pos +1, intNode); 199 changedWays[seg2Way] = true; 200 201 //Do not need to compare again to already split segment 202 seg2Pos++; 203 } 204 205 intersectionNodes.add(intNode); 206 207 if (intNode == newNode) { 208 cmds.add(new AddCommand(dataset, intNode)); 209 } 210 } 211 } else if (test && !intersectionNodes.isEmpty()) 212 return intersectionNodes; 213 } 214 } 215 } 216 } 217 218 219 for (int pos = 0; pos < ways.size(); pos++) { 220 if (!changedWays[pos]) { 221 continue; 222 } 223 224 Way way = ways.get(pos); 225 Way newWay = new Way(way); 226 newWay.setNodes(newNodes[pos]); 227 228 cmds.add(new ChangeCommand(dataset, way, newWay)); 229 } 230 231 return intersectionNodes; 232 } 233 234 private static BBox getNodesBounds(List<Node> nodes) { 235 236 BBox bounds = new BBox(nodes.get(0)); 237 for (Node n: nodes) { 238 bounds.add(n); 239 } 240 return bounds; 241 } 242 243 /** 244 * Tests if given point is to the right side of path consisting of 3 points. 245 * 246 * (Imagine the path is continued beyond the endpoints, so you get two rays 247 * starting from lineP2 and going through lineP1 and lineP3 respectively 248 * which divide the plane into two parts. The test returns true, if testPoint 249 * lies in the part that is to the right when traveling in the direction 250 * lineP1, lineP2, lineP3.) 251 * 252 * @param lineP1 first point in path 253 * @param lineP2 second point in path 254 * @param lineP3 third point in path 255 * @param testPoint point to test 256 * @return true if to the right side, false otherwise 257 */ 258 public static boolean isToTheRightSideOfLine(Node lineP1, Node lineP2, Node lineP3, Node testPoint) { 259 boolean pathBendToRight = angleIsClockwise(lineP1, lineP2, lineP3); 260 boolean rightOfSeg1 = angleIsClockwise(lineP1, lineP2, testPoint); 261 boolean rightOfSeg2 = angleIsClockwise(lineP2, lineP3, testPoint); 262 263 if (pathBendToRight) 264 return rightOfSeg1 && rightOfSeg2; 265 else 266 return !(!rightOfSeg1 && !rightOfSeg2); 267 } 268 269 /** 270 * This method tests if secondNode is clockwise to first node. 271 * @param commonNode starting point for both vectors 272 * @param firstNode first vector end node 273 * @param secondNode second vector end node 274 * @return true if first vector is clockwise before second vector. 275 */ 276 public static boolean angleIsClockwise(Node commonNode, Node firstNode, Node secondNode) { 277 return angleIsClockwise(commonNode.getEastNorth(), firstNode.getEastNorth(), secondNode.getEastNorth()); 278 } 279 280 /** 281 * Finds the intersection of two line segments. 282 * @param p1 the coordinates of the start point of the first specified line segment 283 * @param p2 the coordinates of the end point of the first specified line segment 284 * @param p3 the coordinates of the start point of the second specified line segment 285 * @param p4 the coordinates of the end point of the second specified line segment 286 * @return EastNorth null if no intersection was found, the EastNorth coordinates of the intersection otherwise 287 */ 288 public static EastNorth getSegmentSegmentIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 289 290 CheckParameterUtil.ensure(p1, "p1", EastNorth::isValid); 291 CheckParameterUtil.ensure(p2, "p2", EastNorth::isValid); 292 CheckParameterUtil.ensure(p3, "p3", EastNorth::isValid); 293 CheckParameterUtil.ensure(p4, "p4", EastNorth::isValid); 294 295 double x1 = p1.getX(); 296 double y1 = p1.getY(); 297 double x2 = p2.getX(); 298 double y2 = p2.getY(); 299 double x3 = p3.getX(); 300 double y3 = p3.getY(); 301 double x4 = p4.getX(); 302 double y4 = p4.getY(); 303 304 //TODO: do this locally. 305 //TODO: remove this check after careful testing 306 if (!Line2D.linesIntersect(x1, y1, x2, y2, x3, y3, x4, y4)) return null; 307 308 // solve line-line intersection in parametric form: 309 // (x1,y1) + (x2-x1,y2-y1)* u = (x3,y3) + (x4-x3,y4-y3)* v 310 // (x2-x1,y2-y1)*u - (x4-x3,y4-y3)*v = (x3-x1,y3-y1) 311 // if 0<= u,v <=1, intersection exists at ( x1+ (x2-x1)*u, y1 + (y2-y1)*u ) 312 313 double a1 = x2 - x1; 314 double b1 = x3 - x4; 315 double c1 = x3 - x1; 316 317 double a2 = y2 - y1; 318 double b2 = y3 - y4; 319 double c2 = y3 - y1; 320 321 // Solve the equations 322 double det = a1*b2 - a2*b1; 323 324 double uu = b2*c1 - b1*c2; 325 double vv = a1*c2 - a2*c1; 326 double mag = Math.abs(uu)+Math.abs(vv); 327 328 if (Math.abs(det) > 1e-12 * mag) { 329 double u = uu/det, v = vv/det; 330 if (u > -1e-8 && u < 1+1e-8 && v > -1e-8 && v < 1+1e-8) { 331 if (u < 0) u = 0; 332 if (u > 1) u = 1.0; 333 return new EastNorth(x1+a1*u, y1+a2*u); 334 } else { 335 return null; 336 } 337 } else { 338 // parallel lines 339 return null; 340 } 341 } 342 343 /** 344 * Finds the intersection of two lines of infinite length. 345 * 346 * @param p1 first point on first line 347 * @param p2 second point on first line 348 * @param p3 first point on second line 349 * @param p4 second point on second line 350 * @return EastNorth null if no intersection was found, the coordinates of the intersection otherwise 351 * @throws IllegalArgumentException if a parameter is null or without valid coordinates 352 */ 353 public static EastNorth getLineLineIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 354 355 CheckParameterUtil.ensure(p1, "p1", EastNorth::isValid); 356 CheckParameterUtil.ensure(p2, "p2", EastNorth::isValid); 357 CheckParameterUtil.ensure(p3, "p3", EastNorth::isValid); 358 CheckParameterUtil.ensure(p4, "p4", EastNorth::isValid); 359 360 // Basically, the formula from wikipedia is used: 361 // https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection 362 // However, large numbers lead to rounding errors (see #10286). 363 // To avoid this, p1 is first substracted from each of the points: 364 // p1' = 0 365 // p2' = p2 - p1 366 // p3' = p3 - p1 367 // p4' = p4 - p1 368 // In the end, p1 is added to the intersection point of segment p1'/p2' 369 // and segment p3'/p4'. 370 371 // Convert line from (point, point) form to ax+by=c 372 double a1 = p2.getY() - p1.getY(); 373 double b1 = p1.getX() - p2.getX(); 374 375 double a2 = p4.getY() - p3.getY(); 376 double b2 = p3.getX() - p4.getX(); 377 378 // Solve the equations 379 double det = a1 * b2 - a2 * b1; 380 if (det == 0) 381 return null; // Lines are parallel 382 383 double c2 = (p4.getX() - p1.getX()) * (p3.getY() - p1.getY()) - (p3.getX() - p1.getX()) * (p4.getY() - p1.getY()); 384 385 return new EastNorth(b1 * c2 / det + p1.getX(), -a1 * c2 / det + p1.getY()); 386 } 387 388 /** 389 * Check if the segment p1 - p2 is parallel to p3 - p4 390 * @param p1 First point for first segment 391 * @param p2 Second point for first segment 392 * @param p3 First point for second segment 393 * @param p4 Second point for second segment 394 * @return <code>true</code> if they are parallel or close to parallel 395 */ 396 public static boolean segmentsParallel(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 397 398 CheckParameterUtil.ensure(p1, "p1", EastNorth::isValid); 399 CheckParameterUtil.ensure(p2, "p2", EastNorth::isValid); 400 CheckParameterUtil.ensure(p3, "p3", EastNorth::isValid); 401 CheckParameterUtil.ensure(p4, "p4", EastNorth::isValid); 402 403 // Convert line from (point, point) form to ax+by=c 404 double a1 = p2.getY() - p1.getY(); 405 double b1 = p1.getX() - p2.getX(); 406 407 double a2 = p4.getY() - p3.getY(); 408 double b2 = p3.getX() - p4.getX(); 409 410 // Solve the equations 411 double det = a1 * b2 - a2 * b1; 412 // remove influence of of scaling factor 413 det /= Math.sqrt(a1*a1 + b1*b1) * Math.sqrt(a2*a2 + b2*b2); 414 return Math.abs(det) < 1e-3; 415 } 416 417 private static EastNorth closestPointTo(EastNorth p1, EastNorth p2, EastNorth point, boolean segmentOnly) { 418 CheckParameterUtil.ensureParameterNotNull(p1, "p1"); 419 CheckParameterUtil.ensureParameterNotNull(p2, "p2"); 420 CheckParameterUtil.ensureParameterNotNull(point, "point"); 421 422 double ldx = p2.getX() - p1.getX(); 423 double ldy = p2.getY() - p1.getY(); 424 425 //segment zero length 426 if (ldx == 0 && ldy == 0) 427 return p1; 428 429 double pdx = point.getX() - p1.getX(); 430 double pdy = point.getY() - p1.getY(); 431 432 double offset = (pdx * ldx + pdy * ldy) / (ldx * ldx + ldy * ldy); 433 434 if (segmentOnly && offset <= 0) 435 return p1; 436 else if (segmentOnly && offset >= 1) 437 return p2; 438 else 439 return p1.interpolate(p2, offset); 440 } 441 442 /** 443 * Calculates closest point to a line segment. 444 * @param segmentP1 First point determining line segment 445 * @param segmentP2 Second point determining line segment 446 * @param point Point for which a closest point is searched on line segment [P1,P2] 447 * @return segmentP1 if it is the closest point, segmentP2 if it is the closest point, 448 * a new point if closest point is between segmentP1 and segmentP2. 449 * @see #closestPointToLine 450 * @since 3650 451 */ 452 public static EastNorth closestPointToSegment(EastNorth segmentP1, EastNorth segmentP2, EastNorth point) { 453 return closestPointTo(segmentP1, segmentP2, point, true); 454 } 455 456 /** 457 * Calculates closest point to a line. 458 * @param lineP1 First point determining line 459 * @param lineP2 Second point determining line 460 * @param point Point for which a closest point is searched on line (P1,P2) 461 * @return The closest point found on line. It may be outside the segment [P1,P2]. 462 * @see #closestPointToSegment 463 * @since 4134 464 */ 465 public static EastNorth closestPointToLine(EastNorth lineP1, EastNorth lineP2, EastNorth point) { 466 return closestPointTo(lineP1, lineP2, point, false); 467 } 468 469 /** 470 * This method tests if secondNode is clockwise to first node. 471 * 472 * The line through the two points commonNode and firstNode divides the 473 * plane into two parts. The test returns true, if secondNode lies in 474 * the part that is to the right when traveling in the direction from 475 * commonNode to firstNode. 476 * 477 * @param commonNode starting point for both vectors 478 * @param firstNode first vector end node 479 * @param secondNode second vector end node 480 * @return true if first vector is clockwise before second vector. 481 */ 482 public static boolean angleIsClockwise(EastNorth commonNode, EastNorth firstNode, EastNorth secondNode) { 483 484 CheckParameterUtil.ensure(commonNode, "commonNode", EastNorth::isValid); 485 CheckParameterUtil.ensure(firstNode, "firstNode", EastNorth::isValid); 486 CheckParameterUtil.ensure(secondNode, "secondNode", EastNorth::isValid); 487 488 double dy1 = firstNode.getY() - commonNode.getY(); 489 double dy2 = secondNode.getY() - commonNode.getY(); 490 double dx1 = firstNode.getX() - commonNode.getX(); 491 double dx2 = secondNode.getX() - commonNode.getX(); 492 493 return dy1 * dx2 - dx1 * dy2 > 0; 494 } 495 496 /** 497 * Returns the Area of a polygon, from its list of nodes. 498 * @param polygon List of nodes forming polygon 499 * @return Area for the given list of nodes (EastNorth coordinates) 500 * @since 6841 501 */ 502 public static Area getArea(List<Node> polygon) { 503 Path2D path = new Path2D.Double(); 504 505 boolean begin = true; 506 for (Node n : polygon) { 507 EastNorth en = n.getEastNorth(); 508 if (en != null) { 509 if (begin) { 510 path.moveTo(en.getX(), en.getY()); 511 begin = false; 512 } else { 513 path.lineTo(en.getX(), en.getY()); 514 } 515 } 516 } 517 if (!begin) { 518 path.closePath(); 519 } 520 521 return new Area(path); 522 } 523 524 /** 525 * Builds a path from a list of nodes 526 * @param polygon Nodes, forming a closed polygon 527 * @param path2d path to add to; can be null, then a new path is created 528 * @return the path (LatLon coordinates) 529 */ 530 public static Path2D buildPath2DLatLon(List<Node> polygon, Path2D path2d) { 531 Path2D path = path2d != null ? path2d : new Path2D.Double(); 532 boolean begin = true; 533 for (Node n : polygon) { 534 if (begin) { 535 path.moveTo(n.lon(), n.lat()); 536 begin = false; 537 } else { 538 path.lineTo(n.lon(), n.lat()); 539 } 540 } 541 if (!begin) { 542 path.closePath(); 543 } 544 return path; 545 } 546 547 /** 548 * Returns the Area of a polygon, from the multipolygon relation. 549 * @param multipolygon the multipolygon relation 550 * @return Area for the multipolygon (LatLon coordinates) 551 */ 552 public static Area getAreaLatLon(Relation multipolygon) { 553 final Multipolygon mp = MultipolygonCache.getInstance().get(multipolygon); 554 Path2D path = new Path2D.Double(); 555 path.setWindingRule(Path2D.WIND_EVEN_ODD); 556 for (Multipolygon.PolyData pd : mp.getCombinedPolygons()) { 557 buildPath2DLatLon(pd.getNodes(), path); 558 for (Multipolygon.PolyData pdInner : pd.getInners()) { 559 buildPath2DLatLon(pdInner.getNodes(), path); 560 } 561 } 562 return new Area(path); 563 } 564 565 /** 566 * Tests if two polygons intersect. 567 * @param first List of nodes forming first polygon 568 * @param second List of nodes forming second polygon 569 * @return intersection kind 570 */ 571 public static PolygonIntersection polygonIntersection(List<Node> first, List<Node> second) { 572 Area a1 = getArea(first); 573 Area a2 = getArea(second); 574 return polygonIntersection(a1, a2); 575 } 576 577 /** 578 * Tests if two polygons intersect. 579 * @param a1 Area of first polygon 580 * @param a2 Area of second polygon 581 * @return intersection kind 582 * @since 6841 583 */ 584 public static PolygonIntersection polygonIntersection(Area a1, Area a2) { 585 return polygonIntersection(a1, a2, 1.0); 586 } 587 588 /** 589 * Tests if two polygons intersect. 590 * @param a1 Area of first polygon 591 * @param a2 Area of second polygon 592 * @param eps an area threshold, everything below is considered an empty intersection 593 * @return intersection kind 594 */ 595 public static PolygonIntersection polygonIntersection(Area a1, Area a2, double eps) { 596 597 Area inter = new Area(a1); 598 inter.intersect(a2); 599 600 Rectangle bounds = inter.getBounds(); 601 602 if (inter.isEmpty() || bounds.getHeight()*bounds.getWidth() <= eps) { 603 return PolygonIntersection.OUTSIDE; 604 } else if (a2.getBounds2D().contains(a1.getBounds2D()) && inter.equals(a1)) { 605 return PolygonIntersection.FIRST_INSIDE_SECOND; 606 } else if (a1.getBounds2D().contains(a2.getBounds2D()) && inter.equals(a2)) { 607 return PolygonIntersection.SECOND_INSIDE_FIRST; 608 } else { 609 return PolygonIntersection.CROSSING; 610 } 611 } 612 613 /** 614 * Tests if point is inside a polygon. The polygon can be self-intersecting. In such case the contains function works in xor-like manner. 615 * @param polygonNodes list of nodes from polygon path. 616 * @param point the point to test 617 * @return true if the point is inside polygon. 618 */ 619 public static boolean nodeInsidePolygon(Node point, List<Node> polygonNodes) { 620 if (polygonNodes.size() < 2) 621 return false; 622 623 //iterate each side of the polygon, start with the last segment 624 Node oldPoint = polygonNodes.get(polygonNodes.size() - 1); 625 626 if (!oldPoint.isLatLonKnown()) { 627 return false; 628 } 629 630 boolean inside = false; 631 Node p1, p2; 632 633 for (Node newPoint : polygonNodes) { 634 //skip duplicate points 635 if (newPoint.equals(oldPoint)) { 636 continue; 637 } 638 639 if (!newPoint.isLatLonKnown()) { 640 return false; 641 } 642 643 //order points so p1.lat <= p2.lat 644 if (newPoint.getEastNorth().getY() > oldPoint.getEastNorth().getY()) { 645 p1 = oldPoint; 646 p2 = newPoint; 647 } else { 648 p1 = newPoint; 649 p2 = oldPoint; 650 } 651 652 EastNorth pEN = point.getEastNorth(); 653 EastNorth opEN = oldPoint.getEastNorth(); 654 EastNorth npEN = newPoint.getEastNorth(); 655 EastNorth p1EN = p1.getEastNorth(); 656 EastNorth p2EN = p2.getEastNorth(); 657 658 if (pEN != null && opEN != null && npEN != null && p1EN != null && p2EN != null) { 659 //test if the line is crossed and if so invert the inside flag. 660 if ((npEN.getY() < pEN.getY()) == (pEN.getY() <= opEN.getY()) 661 && (pEN.getX() - p1EN.getX()) * (p2EN.getY() - p1EN.getY()) 662 < (p2EN.getX() - p1EN.getX()) * (pEN.getY() - p1EN.getY())) { 663 inside = !inside; 664 } 665 } 666 667 oldPoint = newPoint; 668 } 669 670 return inside; 671 } 672 673 /** 674 * Returns area of a closed way in square meters. 675 * 676 * @param way Way to measure, should be closed (first node is the same as last node) 677 * @return area of the closed way. 678 */ 679 public static double closedWayArea(Way way) { 680 return getAreaAndPerimeter(way.getNodes(), Projections.getProjectionByCode("EPSG:54008")).getArea(); 681 } 682 683 /** 684 * Returns area of a multipolygon in square meters. 685 * 686 * @param multipolygon the multipolygon to measure 687 * @return area of the multipolygon. 688 */ 689 public static double multipolygonArea(Relation multipolygon) { 690 double area = 0.0; 691 final Multipolygon mp = MultipolygonCache.getInstance().get(multipolygon); 692 for (Multipolygon.PolyData pd : mp.getCombinedPolygons()) { 693 area += pd.getAreaAndPerimeter(Projections.getProjectionByCode("EPSG:54008")).getArea(); 694 } 695 return area; 696 } 697 698 /** 699 * Computes the area of a closed way and multipolygon in square meters, or {@code null} for other primitives 700 * 701 * @param osm the primitive to measure 702 * @return area of the primitive, or {@code null} 703 */ 704 public static Double computeArea(OsmPrimitive osm) { 705 if (osm instanceof Way && ((Way) osm).isClosed()) { 706 return closedWayArea((Way) osm); 707 } else if (osm instanceof Relation && ((Relation) osm).isMultipolygon() && !((Relation) osm).hasIncompleteMembers()) { 708 return multipolygonArea((Relation) osm); 709 } else { 710 return null; 711 } 712 } 713 714 /** 715 * Determines whether a way is oriented clockwise. 716 * 717 * Internals: Assuming a closed non-looping way, compute twice the area 718 * of the polygon using the formula {@code 2 * area = sum (X[n] * Y[n+1] - X[n+1] * Y[n])}. 719 * If the area is negative the way is ordered in a clockwise direction. 720 * 721 * See http://paulbourke.net/geometry/polyarea/ 722 * 723 * @param w the way to be checked. 724 * @return true if and only if way is oriented clockwise. 725 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}). 726 */ 727 public static boolean isClockwise(Way w) { 728 return isClockwise(w.getNodes()); 729 } 730 731 /** 732 * Determines whether path from nodes list is oriented clockwise. 733 * @param nodes Nodes list to be checked. 734 * @return true if and only if way is oriented clockwise. 735 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}). 736 * @see #isClockwise(Way) 737 */ 738 public static boolean isClockwise(List<Node> nodes) { 739 int nodesCount = nodes.size(); 740 if (nodesCount < 3 || nodes.get(0) != nodes.get(nodesCount - 1)) { 741 throw new IllegalArgumentException("Way must be closed to check orientation."); 742 } 743 double area2 = 0.; 744 745 for (int node = 1; node <= /*sic! consider last-first as well*/ nodesCount; node++) { 746 Node coorPrev = nodes.get(node - 1); 747 Node coorCurr = nodes.get(node % nodesCount); 748 area2 += coorPrev.lon() * coorCurr.lat(); 749 area2 -= coorCurr.lon() * coorPrev.lat(); 750 } 751 return area2 < 0; 752 } 753 754 /** 755 * Returns angle of a segment defined with 2 point coordinates. 756 * 757 * @param p1 first point 758 * @param p2 second point 759 * @return Angle in radians (-pi, pi] 760 */ 761 public static double getSegmentAngle(EastNorth p1, EastNorth p2) { 762 763 CheckParameterUtil.ensure(p1, "p1", EastNorth::isValid); 764 CheckParameterUtil.ensure(p2, "p2", EastNorth::isValid); 765 766 return Math.atan2(p2.north() - p1.north(), p2.east() - p1.east()); 767 } 768 769 /** 770 * Returns angle of a corner defined with 3 point coordinates. 771 * 772 * @param p1 first point 773 * @param p2 Common endpoint 774 * @param p3 third point 775 * @return Angle in radians (-pi, pi] 776 */ 777 public static double getCornerAngle(EastNorth p1, EastNorth p2, EastNorth p3) { 778 779 CheckParameterUtil.ensure(p1, "p1", EastNorth::isValid); 780 CheckParameterUtil.ensure(p2, "p2", EastNorth::isValid); 781 CheckParameterUtil.ensure(p3, "p3", EastNorth::isValid); 782 783 Double result = getSegmentAngle(p2, p1) - getSegmentAngle(p2, p3); 784 if (result <= -Math.PI) { 785 result += 2 * Math.PI; 786 } 787 788 if (result > Math.PI) { 789 result -= 2 * Math.PI; 790 } 791 792 return result; 793 } 794 795 /** 796 * Compute the centroid/barycenter of nodes 797 * @param nodes Nodes for which the centroid is wanted 798 * @return the centroid of nodes 799 * @see Geometry#getCenter 800 */ 801 public static EastNorth getCentroid(List<Node> nodes) { 802 803 BigDecimal area = BigDecimal.ZERO; 804 BigDecimal north = BigDecimal.ZERO; 805 BigDecimal east = BigDecimal.ZERO; 806 807 // See https://en.wikipedia.org/wiki/Centroid#Centroid_of_polygon for the equation used here 808 for (int i = 0; i < nodes.size(); i++) { 809 EastNorth n0 = nodes.get(i).getEastNorth(); 810 EastNorth n1 = nodes.get((i+1) % nodes.size()).getEastNorth(); 811 812 if (n0 != null && n1 != null && n0.isValid() && n1.isValid()) { 813 BigDecimal x0 = BigDecimal.valueOf(n0.east()); 814 BigDecimal y0 = BigDecimal.valueOf(n0.north()); 815 BigDecimal x1 = BigDecimal.valueOf(n1.east()); 816 BigDecimal y1 = BigDecimal.valueOf(n1.north()); 817 818 BigDecimal k = x0.multiply(y1, MathContext.DECIMAL128).subtract(y0.multiply(x1, MathContext.DECIMAL128)); 819 820 area = area.add(k, MathContext.DECIMAL128); 821 east = east.add(k.multiply(x0.add(x1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 822 north = north.add(k.multiply(y0.add(y1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 823 } 824 } 825 826 BigDecimal d = new BigDecimal(3, MathContext.DECIMAL128); // 1/2 * 6 = 3 827 area = area.multiply(d, MathContext.DECIMAL128); 828 if (area.compareTo(BigDecimal.ZERO) != 0) { 829 north = north.divide(area, MathContext.DECIMAL128); 830 east = east.divide(area, MathContext.DECIMAL128); 831 } 832 833 return new EastNorth(east.doubleValue(), north.doubleValue()); 834 } 835 836 /** 837 * Compute center of the circle closest to different nodes. 838 * 839 * Ensure exact center computation in case nodes are already aligned in circle. 840 * This is done by least square method. 841 * Let be a_i x + b_i y + c_i = 0 equations of bisectors of each edges. 842 * Center must be intersection of all bisectors. 843 * <pre> 844 * [ a1 b1 ] [ -c1 ] 845 * With A = [ ... ... ] and Y = [ ... ] 846 * [ an bn ] [ -cn ] 847 * </pre> 848 * An approximation of center of circle is (At.A)^-1.At.Y 849 * @param nodes Nodes parts of the circle (at least 3) 850 * @return An approximation of the center, of null if there is no solution. 851 * @see Geometry#getCentroid 852 * @since 6934 853 */ 854 public static EastNorth getCenter(List<Node> nodes) { 855 int nc = nodes.size(); 856 if (nc < 3) return null; 857 /** 858 * Equation of each bisector ax + by + c = 0 859 */ 860 double[] a = new double[nc]; 861 double[] b = new double[nc]; 862 double[] c = new double[nc]; 863 // Compute equation of bisector 864 for (int i = 0; i < nc; i++) { 865 EastNorth pt1 = nodes.get(i).getEastNorth(); 866 EastNorth pt2 = nodes.get((i+1) % nc).getEastNorth(); 867 a[i] = pt1.east() - pt2.east(); 868 b[i] = pt1.north() - pt2.north(); 869 double d = Math.sqrt(a[i]*a[i] + b[i]*b[i]); 870 if (d == 0) return null; 871 a[i] /= d; 872 b[i] /= d; 873 double xC = (pt1.east() + pt2.east()) / 2; 874 double yC = (pt1.north() + pt2.north()) / 2; 875 c[i] = -(a[i]*xC + b[i]*yC); 876 } 877 // At.A = [aij] 878 double a11 = 0, a12 = 0, a22 = 0; 879 // At.Y = [bi] 880 double b1 = 0, b2 = 0; 881 for (int i = 0; i < nc; i++) { 882 a11 += a[i]*a[i]; 883 a12 += a[i]*b[i]; 884 a22 += b[i]*b[i]; 885 b1 -= a[i]*c[i]; 886 b2 -= b[i]*c[i]; 887 } 888 // (At.A)^-1 = [invij] 889 double det = a11*a22 - a12*a12; 890 if (Math.abs(det) < 1e-5) return null; 891 double inv11 = a22/det; 892 double inv12 = -a12/det; 893 double inv22 = a11/det; 894 // center (xC, yC) = (At.A)^-1.At.y 895 double xC = inv11*b1 + inv12*b2; 896 double yC = inv12*b1 + inv22*b2; 897 return new EastNorth(xC, yC); 898 } 899 900 /** 901 * Tests if the {@code node} is inside the multipolygon {@code multiPolygon}. The nullable argument 902 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match. 903 * @param node node 904 * @param multiPolygon multipolygon 905 * @param isOuterWayAMatch allows to decide if the immediate {@code outer} way of the multipolygon is a match 906 * @return {@code true} if the node is inside the multipolygon 907 */ 908 public static boolean isNodeInsideMultiPolygon(Node node, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) { 909 return isPolygonInsideMultiPolygon(Collections.singletonList(node), multiPolygon, isOuterWayAMatch); 910 } 911 912 /** 913 * Tests if the polygon formed by {@code nodes} is inside the multipolygon {@code multiPolygon}. The nullable argument 914 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match. 915 * <p> 916 * If {@code nodes} contains exactly one element, then it is checked whether that one node is inside the multipolygon. 917 * @param nodes nodes forming the polygon 918 * @param multiPolygon multipolygon 919 * @param isOuterWayAMatch allows to decide if the immediate {@code outer} way of the multipolygon is a match 920 * @return {@code true} if the polygon formed by nodes is inside the multipolygon 921 */ 922 public static boolean isPolygonInsideMultiPolygon(List<Node> nodes, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) { 923 // Extract outer/inner members from multipolygon 924 final Pair<List<JoinedPolygon>, List<JoinedPolygon>> outerInner; 925 try { 926 outerInner = MultipolygonBuilder.joinWays(multiPolygon); 927 } catch (MultipolygonBuilder.JoinedPolygonCreationException ex) { 928 Logging.trace(ex); 929 Logging.debug("Invalid multipolygon " + multiPolygon); 930 return false; 931 } 932 // Test if object is inside an outer member 933 for (JoinedPolygon out : outerInner.a) { 934 if (nodes.size() == 1 935 ? nodeInsidePolygon(nodes.get(0), out.getNodes()) 936 : EnumSet.of(PolygonIntersection.FIRST_INSIDE_SECOND, PolygonIntersection.CROSSING).contains( 937 polygonIntersection(nodes, out.getNodes()))) { 938 boolean insideInner = false; 939 // If inside an outer, check it is not inside an inner 940 for (JoinedPolygon in : outerInner.b) { 941 if (polygonIntersection(in.getNodes(), out.getNodes()) == PolygonIntersection.FIRST_INSIDE_SECOND 942 && (nodes.size() == 1 943 ? nodeInsidePolygon(nodes.get(0), in.getNodes()) 944 : polygonIntersection(nodes, in.getNodes()) == PolygonIntersection.FIRST_INSIDE_SECOND)) { 945 insideInner = true; 946 break; 947 } 948 } 949 // Inside outer but not inside inner -> the polygon appears to be inside a the multipolygon 950 if (!insideInner) { 951 // Final check using predicate 952 if (isOuterWayAMatch == null || isOuterWayAMatch.test(out.ways.get(0) 953 /* TODO give a better representation of the outer ring to the predicate */)) { 954 return true; 955 } 956 } 957 } 958 } 959 return false; 960 } 961 962 /** 963 * Data class to hold two double values (area and perimeter of a polygon). 964 */ 965 public static class AreaAndPerimeter { 966 private final double area; 967 private final double perimeter; 968 969 /** 970 * Create a new {@link AreaAndPerimeter} 971 * @param area The area 972 * @param perimeter The perimeter 973 */ 974 public AreaAndPerimeter(double area, double perimeter) { 975 this.area = area; 976 this.perimeter = perimeter; 977 } 978 979 /** 980 * Gets the area 981 * @return The area size 982 */ 983 public double getArea() { 984 return area; 985 } 986 987 /** 988 * Gets the perimeter 989 * @return The perimeter length 990 */ 991 public double getPerimeter() { 992 return perimeter; 993 } 994 } 995 996 /** 997 * Calculate area and perimeter length of a polygon. 998 * 999 * Uses current projection; units are that of the projected coordinates. 1000 * 1001 * @param nodes the list of nodes representing the polygon 1002 * @return area and perimeter 1003 */ 1004 public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes) { 1005 return getAreaAndPerimeter(nodes, null); 1006 } 1007 1008 /** 1009 * Calculate area and perimeter length of a polygon in the given projection. 1010 * 1011 * @param nodes the list of nodes representing the polygon 1012 * @param projection the projection to use for the calculation, {@code null} defaults to {@link Main#getProjection()} 1013 * @return area and perimeter 1014 */ 1015 public static AreaAndPerimeter getAreaAndPerimeter(List<Node> nodes, Projection projection) { 1016 CheckParameterUtil.ensureParameterNotNull(nodes, "nodes"); 1017 double area = 0; 1018 double perimeter = 0; 1019 Projection useProjection = projection == null ? Main.getProjection() : projection; 1020 1021 if (!nodes.isEmpty()) { 1022 boolean closed = nodes.get(0) == nodes.get(nodes.size() - 1); 1023 int numSegments = closed ? nodes.size() - 1 : nodes.size(); 1024 EastNorth p1 = nodes.get(0).getEastNorth(useProjection); 1025 for (int i = 1; i <= numSegments; i++) { 1026 final Node node = nodes.get(i == numSegments ? 0 : i); 1027 final EastNorth p2 = node.getEastNorth(useProjection); 1028 if (p1 != null && p2 != null) { 1029 area += p1.east() * p2.north() - p2.east() * p1.north(); 1030 perimeter += p1.distance(p2); 1031 } 1032 p1 = p2; 1033 } 1034 } 1035 return new AreaAndPerimeter(Math.abs(area) / 2, perimeter); 1036 } 1037}