001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.tools; 003 004import static org.openstreetmap.josm.data.projection.Ellipsoid.WGS84; 005 006import java.awt.geom.Area; 007import java.awt.geom.Line2D; 008import java.awt.geom.Path2D; 009import java.awt.geom.PathIterator; 010import java.awt.geom.Rectangle2D; 011import java.math.BigDecimal; 012import java.math.MathContext; 013import java.util.ArrayList; 014import java.util.Collection; 015import java.util.Collections; 016import java.util.Comparator; 017import java.util.Iterator; 018import java.util.LinkedHashSet; 019import java.util.List; 020import java.util.Set; 021import java.util.TreeSet; 022import java.util.function.Predicate; 023import java.util.stream.Collectors; 024 025import org.openstreetmap.josm.command.AddCommand; 026import org.openstreetmap.josm.command.ChangeNodesCommand; 027import org.openstreetmap.josm.command.Command; 028import org.openstreetmap.josm.data.coor.EastNorth; 029import org.openstreetmap.josm.data.coor.ILatLon; 030import org.openstreetmap.josm.data.coor.LatLon; 031import org.openstreetmap.josm.data.osm.BBox; 032import org.openstreetmap.josm.data.osm.DataSet; 033import org.openstreetmap.josm.data.osm.INode; 034import org.openstreetmap.josm.data.osm.IPrimitive; 035import org.openstreetmap.josm.data.osm.IWay; 036import org.openstreetmap.josm.data.osm.MultipolygonBuilder; 037import org.openstreetmap.josm.data.osm.MultipolygonBuilder.JoinedPolygon; 038import org.openstreetmap.josm.data.osm.Node; 039import org.openstreetmap.josm.data.osm.NodePositionComparator; 040import org.openstreetmap.josm.data.osm.OsmPrimitive; 041import org.openstreetmap.josm.data.osm.Relation; 042import org.openstreetmap.josm.data.osm.Way; 043import org.openstreetmap.josm.data.osm.WaySegment; 044import org.openstreetmap.josm.data.osm.visitor.paint.relations.Multipolygon; 045import org.openstreetmap.josm.data.osm.visitor.paint.relations.Multipolygon.PolyData; 046import org.openstreetmap.josm.data.osm.visitor.paint.relations.MultipolygonCache; 047import org.openstreetmap.josm.data.projection.Projection; 048import org.openstreetmap.josm.data.projection.ProjectionRegistry; 049import org.openstreetmap.josm.data.projection.Projections; 050 051/** 052 * Some tools for geometry related tasks. 053 * 054 * @author viesturs 055 */ 056public final class Geometry { 057 058 private Geometry() { 059 // Hide default constructor for utils classes 060 } 061 062 /** 063 * The result types for a {@link Geometry#polygonIntersection(Area, Area)} test 064 */ 065 public enum PolygonIntersection { 066 /** 067 * The first polygon is inside the second one 068 */ 069 FIRST_INSIDE_SECOND, 070 /** 071 * The second one is inside the first 072 */ 073 SECOND_INSIDE_FIRST, 074 /** 075 * The polygons do not overlap 076 */ 077 OUTSIDE, 078 /** 079 * The polygon borders cross each other 080 */ 081 CROSSING 082 } 083 084 /** threshold value for size of intersection area given in east/north space */ 085 public static final double INTERSECTION_EPS_EAST_NORTH = 1e-4; 086 087 /** 088 * Will find all intersection and add nodes there for list of given ways. 089 * Handles self-intersections too. 090 * And makes commands to add the intersection points to ways. 091 * 092 * Prerequisite: no two nodes have the same coordinates. 093 * 094 * @param ways a list of ways to test 095 * @param test if true, do not build list of Commands, just return nodes 096 * @param cmds list of commands, typically empty when handed to this method. 097 * Will be filled with commands that add intersection nodes to 098 * the ways. 099 * @return list of new nodes, if test is true the list might not contain all intersections 100 */ 101 public static Set<Node> addIntersections(List<Way> ways, boolean test, List<Command> cmds) { 102 103 int n = ways.size(); 104 @SuppressWarnings("unchecked") 105 List<Node>[] newNodes = new ArrayList[n]; 106 BBox[] wayBounds = new BBox[n]; 107 boolean[] changedWays = new boolean[n]; 108 109 Set<Node> intersectionNodes = new LinkedHashSet<>(); 110 111 //copy node arrays for local usage. 112 for (int pos = 0; pos < n; pos++) { 113 newNodes[pos] = new ArrayList<>(ways.get(pos).getNodes()); 114 wayBounds[pos] = ways.get(pos).getBBox(); 115 changedWays[pos] = false; 116 } 117 118 DataSet dataset = ways.get(0).getDataSet(); 119 120 //iterate over all way pairs and introduce the intersections 121 Comparator<Node> coordsComparator = new NodePositionComparator(); 122 for (int seg1Way = 0; seg1Way < n; seg1Way++) { 123 for (int seg2Way = seg1Way; seg2Way < n; seg2Way++) { 124 125 //do not waste time on bounds that do not intersect 126 if (!wayBounds[seg1Way].intersects(wayBounds[seg2Way])) { 127 continue; 128 } 129 130 List<Node> way1Nodes = newNodes[seg1Way]; 131 List<Node> way2Nodes = newNodes[seg2Way]; 132 133 //iterate over primary segmemt 134 for (int seg1Pos = 0; seg1Pos + 1 < way1Nodes.size(); seg1Pos++) { 135 136 //iterate over secondary segment 137 int seg2Start = seg1Way != seg2Way ? 0 : seg1Pos + 2; //skip the adjacent segment 138 139 for (int seg2Pos = seg2Start; seg2Pos + 1 < way2Nodes.size(); seg2Pos++) { 140 141 //need to get them again every time, because other segments may be changed 142 Node seg1Node1 = way1Nodes.get(seg1Pos); 143 Node seg1Node2 = way1Nodes.get(seg1Pos + 1); 144 Node seg2Node1 = way2Nodes.get(seg2Pos); 145 Node seg2Node2 = way2Nodes.get(seg2Pos + 1); 146 147 int commonCount = 0; 148 //test if we have common nodes to add. 149 if (seg1Node1 == seg2Node1 || seg1Node1 == seg2Node2) { 150 commonCount++; 151 152 if (seg1Way == seg2Way && 153 seg1Pos == 0 && 154 seg2Pos == way2Nodes.size() -2) { 155 //do not add - this is first and last segment of the same way. 156 } else { 157 intersectionNodes.add(seg1Node1); 158 } 159 } 160 161 if (seg1Node2 == seg2Node1 || seg1Node2 == seg2Node2) { 162 commonCount++; 163 164 intersectionNodes.add(seg1Node2); 165 } 166 167 //no common nodes - find intersection 168 if (commonCount == 0) { 169 EastNorth intersection = getSegmentSegmentIntersection( 170 seg1Node1.getEastNorth(), seg1Node2.getEastNorth(), 171 seg2Node1.getEastNorth(), seg2Node2.getEastNorth()); 172 173 if (intersection != null) { 174 Node newNode = new Node(ProjectionRegistry.getProjection().eastNorth2latlon(intersection)); 175 Node intNode = newNode; 176 boolean insertInSeg1 = false; 177 boolean insertInSeg2 = false; 178 //find if the intersection point is at end point of one of the segments, if so use that point 179 180 //segment 1 181 if (coordsComparator.compare(newNode, seg1Node1) == 0) { 182 intNode = seg1Node1; 183 } else if (coordsComparator.compare(newNode, seg1Node2) == 0) { 184 intNode = seg1Node2; 185 } else { 186 insertInSeg1 = true; 187 } 188 189 //segment 2 190 if (coordsComparator.compare(newNode, seg2Node1) == 0) { 191 intNode = seg2Node1; 192 } else if (coordsComparator.compare(newNode, seg2Node2) == 0) { 193 intNode = seg2Node2; 194 } else { 195 insertInSeg2 = true; 196 } 197 198 if (test) { 199 intersectionNodes.add(intNode); 200 return intersectionNodes; 201 } 202 203 if (insertInSeg1) { 204 way1Nodes.add(seg1Pos +1, intNode); 205 changedWays[seg1Way] = true; 206 207 //fix seg2 position, as indexes have changed, seg2Pos is always bigger than seg1Pos on the same segment. 208 if (seg2Way == seg1Way) { 209 seg2Pos++; 210 } 211 } 212 213 if (insertInSeg2) { 214 way2Nodes.add(seg2Pos +1, intNode); 215 changedWays[seg2Way] = true; 216 217 //Do not need to compare again to already split segment 218 seg2Pos++; 219 } 220 221 intersectionNodes.add(intNode); 222 223 if (intNode == newNode) { 224 cmds.add(new AddCommand(dataset, intNode)); 225 } 226 } 227 } else if (test && !intersectionNodes.isEmpty()) 228 return intersectionNodes; 229 } 230 } 231 } 232 } 233 234 235 for (int pos = 0; pos < ways.size(); pos++) { 236 if (changedWays[pos]) { 237 cmds.add(new ChangeNodesCommand(dataset, ways.get(pos), newNodes[pos])); 238 } 239 } 240 241 return intersectionNodes; 242 } 243 244 /** 245 * Tests if given point is to the right side of path consisting of 3 points. 246 * 247 * (Imagine the path is continued beyond the endpoints, so you get two rays 248 * starting from lineP2 and going through lineP1 and lineP3 respectively 249 * which divide the plane into two parts. The test returns true, if testPoint 250 * lies in the part that is to the right when traveling in the direction 251 * lineP1, lineP2, lineP3.) 252 * 253 * @param <N> type of node 254 * @param lineP1 first point in path 255 * @param lineP2 second point in path 256 * @param lineP3 third point in path 257 * @param testPoint point to test 258 * @return true if to the right side, false otherwise 259 */ 260 public static <N extends INode> boolean isToTheRightSideOfLine(N lineP1, N lineP2, N lineP3, N testPoint) { 261 boolean pathBendToRight = angleIsClockwise(lineP1, lineP2, lineP3); 262 boolean rightOfSeg1 = angleIsClockwise(lineP1, lineP2, testPoint); 263 boolean rightOfSeg2 = angleIsClockwise(lineP2, lineP3, testPoint); 264 265 if (pathBendToRight) 266 return rightOfSeg1 && rightOfSeg2; 267 else 268 return !(!rightOfSeg1 && !rightOfSeg2); 269 } 270 271 /** 272 * This method tests if secondNode is clockwise to first node. 273 * @param <N> type of node 274 * @param commonNode starting point for both vectors 275 * @param firstNode first vector end node 276 * @param secondNode second vector end node 277 * @return true if first vector is clockwise before second vector. 278 */ 279 public static <N extends INode> boolean angleIsClockwise(N commonNode, N firstNode, N secondNode) { 280 return angleIsClockwise(commonNode.getEastNorth(), firstNode.getEastNorth(), secondNode.getEastNorth()); 281 } 282 283 /** 284 * Finds the intersection of two line segments. 285 * @param p1 the coordinates of the start point of the first specified line segment 286 * @param p2 the coordinates of the end point of the first specified line segment 287 * @param p3 the coordinates of the start point of the second specified line segment 288 * @param p4 the coordinates of the end point of the second specified line segment 289 * @return EastNorth null if no intersection was found, the EastNorth coordinates of the intersection otherwise 290 */ 291 public static EastNorth getSegmentSegmentIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 292 293 CheckParameterUtil.ensureThat(p1.isValid(), () -> p1 + " invalid"); 294 CheckParameterUtil.ensureThat(p2.isValid(), () -> p2 + " invalid"); 295 CheckParameterUtil.ensureThat(p3.isValid(), () -> p3 + " invalid"); 296 CheckParameterUtil.ensureThat(p4.isValid(), () -> p4 + " invalid"); 297 298 double x1 = p1.getX(); 299 double y1 = p1.getY(); 300 double x2 = p2.getX(); 301 double y2 = p2.getY(); 302 double x3 = p3.getX(); 303 double y3 = p3.getY(); 304 double x4 = p4.getX(); 305 double y4 = p4.getY(); 306 307 //TODO: do this locally. 308 //TODO: remove this check after careful testing 309 if (!Line2D.linesIntersect(x1, y1, x2, y2, x3, y3, x4, y4)) return null; 310 311 // solve line-line intersection in parametric form: 312 // (x1,y1) + (x2-x1,y2-y1)* u = (x3,y3) + (x4-x3,y4-y3)* v 313 // (x2-x1,y2-y1)*u - (x4-x3,y4-y3)*v = (x3-x1,y3-y1) 314 // if 0<= u,v <=1, intersection exists at ( x1+ (x2-x1)*u, y1 + (y2-y1)*u ) 315 316 double a1 = x2 - x1; 317 double b1 = x3 - x4; 318 double c1 = x3 - x1; 319 320 double a2 = y2 - y1; 321 double b2 = y3 - y4; 322 double c2 = y3 - y1; 323 324 // Solve the equations 325 double det = a1*b2 - a2*b1; 326 327 double uu = b2*c1 - b1*c2; 328 double vv = a1*c2 - a2*c1; 329 double mag = Math.abs(uu)+Math.abs(vv); 330 331 if (Math.abs(det) > 1e-12 * mag) { 332 double u = uu/det, v = vv/det; 333 if (u > -1e-8 && u < 1+1e-8 && v > -1e-8 && v < 1+1e-8) { 334 if (u < 0) u = 0; 335 if (u > 1) u = 1.0; 336 return new EastNorth(x1+a1*u, y1+a2*u); 337 } else { 338 return null; 339 } 340 } else { 341 // parallel lines 342 return null; 343 } 344 } 345 346 /** 347 * Finds the intersection of two lines of infinite length. 348 * 349 * @param p1 first point on first line 350 * @param p2 second point on first line 351 * @param p3 first point on second line 352 * @param p4 second point on second line 353 * @return EastNorth null if no intersection was found, the coordinates of the intersection otherwise 354 * @throws IllegalArgumentException if a parameter is null or without valid coordinates 355 */ 356 public static EastNorth getLineLineIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 357 358 CheckParameterUtil.ensureThat(p1.isValid(), () -> p1 + " invalid"); 359 CheckParameterUtil.ensureThat(p2.isValid(), () -> p2 + " invalid"); 360 CheckParameterUtil.ensureThat(p3.isValid(), () -> p3 + " invalid"); 361 CheckParameterUtil.ensureThat(p4.isValid(), () -> p4 + " invalid"); 362 363 // Basically, the formula from wikipedia is used: 364 // https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection 365 // However, large numbers lead to rounding errors (see #10286). 366 // To avoid this, p1 is first subtracted from each of the points: 367 // p1' = 0 368 // p2' = p2 - p1 369 // p3' = p3 - p1 370 // p4' = p4 - p1 371 // In the end, p1 is added to the intersection point of segment p1'/p2' 372 // and segment p3'/p4'. 373 374 // Convert line from (point, point) form to ax+by=c 375 double a1 = p2.getY() - p1.getY(); 376 double b1 = p1.getX() - p2.getX(); 377 378 double a2 = p4.getY() - p3.getY(); 379 double b2 = p3.getX() - p4.getX(); 380 381 // Solve the equations 382 double det = a1 * b2 - a2 * b1; 383 if (det == 0) 384 return null; // Lines are parallel 385 386 double c2 = (p4.getX() - p1.getX()) * (p3.getY() - p1.getY()) - (p3.getX() - p1.getX()) * (p4.getY() - p1.getY()); 387 388 return new EastNorth(b1 * c2 / det + p1.getX(), -a1 * c2 / det + p1.getY()); 389 } 390 391 /** 392 * Check if the segment p1 - p2 is parallel to p3 - p4 393 * @param p1 First point for first segment 394 * @param p2 Second point for first segment 395 * @param p3 First point for second segment 396 * @param p4 Second point for second segment 397 * @return <code>true</code> if they are parallel or close to parallel 398 */ 399 public static boolean segmentsParallel(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 400 401 CheckParameterUtil.ensureThat(p1.isValid(), () -> p1 + " invalid"); 402 CheckParameterUtil.ensureThat(p2.isValid(), () -> p2 + " invalid"); 403 CheckParameterUtil.ensureThat(p3.isValid(), () -> p3 + " invalid"); 404 CheckParameterUtil.ensureThat(p4.isValid(), () -> p4 + " invalid"); 405 406 // Convert line from (point, point) form to ax+by=c 407 double a1 = p2.getY() - p1.getY(); 408 double b1 = p1.getX() - p2.getX(); 409 410 double a2 = p4.getY() - p3.getY(); 411 double b2 = p3.getX() - p4.getX(); 412 413 // Solve the equations 414 double det = a1 * b2 - a2 * b1; 415 // remove influence of of scaling factor 416 det /= Math.sqrt(a1*a1 + b1*b1) * Math.sqrt(a2*a2 + b2*b2); 417 return Math.abs(det) < 1e-3; 418 } 419 420 private static EastNorth closestPointTo(EastNorth p1, EastNorth p2, EastNorth point, boolean segmentOnly) { 421 CheckParameterUtil.ensureParameterNotNull(p1, "p1"); 422 CheckParameterUtil.ensureParameterNotNull(p2, "p2"); 423 CheckParameterUtil.ensureParameterNotNull(point, "point"); 424 425 double ldx = p2.getX() - p1.getX(); 426 double ldy = p2.getY() - p1.getY(); 427 428 //segment zero length 429 if (ldx == 0 && ldy == 0) 430 return p1; 431 432 double pdx = point.getX() - p1.getX(); 433 double pdy = point.getY() - p1.getY(); 434 435 double offset = (pdx * ldx + pdy * ldy) / (ldx * ldx + ldy * ldy); 436 437 if (segmentOnly && offset <= 0) 438 return p1; 439 else if (segmentOnly && offset >= 1) 440 return p2; 441 else 442 return p1.interpolate(p2, offset); 443 } 444 445 /** 446 * Calculates closest point to a line segment. 447 * @param segmentP1 First point determining line segment 448 * @param segmentP2 Second point determining line segment 449 * @param point Point for which a closest point is searched on line segment [P1,P2] 450 * @return segmentP1 if it is the closest point, segmentP2 if it is the closest point, 451 * a new point if closest point is between segmentP1 and segmentP2. 452 * @see #closestPointToLine 453 * @since 3650 454 */ 455 public static EastNorth closestPointToSegment(EastNorth segmentP1, EastNorth segmentP2, EastNorth point) { 456 return closestPointTo(segmentP1, segmentP2, point, true); 457 } 458 459 /** 460 * Calculates closest point to a line. 461 * @param lineP1 First point determining line 462 * @param lineP2 Second point determining line 463 * @param point Point for which a closest point is searched on line (P1,P2) 464 * @return The closest point found on line. It may be outside the segment [P1,P2]. 465 * @see #closestPointToSegment 466 * @since 4134 467 */ 468 public static EastNorth closestPointToLine(EastNorth lineP1, EastNorth lineP2, EastNorth point) { 469 return closestPointTo(lineP1, lineP2, point, false); 470 } 471 472 /** 473 * This method tests if secondNode is clockwise to first node. 474 * 475 * The line through the two points commonNode and firstNode divides the 476 * plane into two parts. The test returns true, if secondNode lies in 477 * the part that is to the right when traveling in the direction from 478 * commonNode to firstNode. 479 * 480 * @param commonNode starting point for both vectors 481 * @param firstNode first vector end node 482 * @param secondNode second vector end node 483 * @return true if first vector is clockwise before second vector. 484 */ 485 public static boolean angleIsClockwise(EastNorth commonNode, EastNorth firstNode, EastNorth secondNode) { 486 487 CheckParameterUtil.ensureThat(commonNode.isValid(), () -> commonNode + " invalid"); 488 CheckParameterUtil.ensureThat(firstNode.isValid(), () -> firstNode + " invalid"); 489 CheckParameterUtil.ensureThat(secondNode.isValid(), () -> secondNode + " invalid"); 490 491 double dy1 = firstNode.getY() - commonNode.getY(); 492 double dy2 = secondNode.getY() - commonNode.getY(); 493 double dx1 = firstNode.getX() - commonNode.getX(); 494 double dx2 = secondNode.getX() - commonNode.getX(); 495 496 return dy1 * dx2 - dx1 * dy2 > 0; 497 } 498 499 /** 500 * Returns the Area of a polygon, from its list of nodes. 501 * @param polygon List of nodes forming polygon 502 * @return Area for the given list of nodes (EastNorth coordinates) 503 * @since 6841 504 */ 505 public static Area getArea(List<? extends INode> polygon) { 506 Path2D path = new Path2D.Double(); 507 508 boolean begin = true; 509 for (INode n : polygon) { 510 EastNorth en = n.getEastNorth(); 511 if (en != null) { 512 if (begin) { 513 path.moveTo(en.getX(), en.getY()); 514 begin = false; 515 } else { 516 path.lineTo(en.getX(), en.getY()); 517 } 518 } 519 } 520 if (!begin) { 521 path.closePath(); 522 } 523 524 return new Area(path); 525 } 526 527 /** 528 * Builds a path from a list of nodes 529 * @param polygon Nodes, forming a closed polygon 530 * @param path2d path to add to; can be null, then a new path is created 531 * @return the path (LatLon coordinates) 532 * @since 13638 (signature) 533 */ 534 public static Path2D buildPath2DLatLon(List<? extends ILatLon> polygon, Path2D path2d) { 535 Path2D path = path2d != null ? path2d : new Path2D.Double(); 536 boolean begin = true; 537 for (ILatLon n : polygon) { 538 if (begin) { 539 path.moveTo(n.lon(), n.lat()); 540 begin = false; 541 } else { 542 path.lineTo(n.lon(), n.lat()); 543 } 544 } 545 if (!begin) { 546 path.closePath(); 547 } 548 return path; 549 } 550 551 /** 552 * Calculate area in east/north space for given primitive. Uses {@link MultipolygonCache} for multipolygon relations. 553 * @param p the primitive 554 * @return the area in east/north space, might be empty if the primitive is incomplete or not closed or a node 555 * since 15938 556 */ 557 public static Area getAreaEastNorth(IPrimitive p) { 558 if (p instanceof Way && ((Way) p).isClosed()) { 559 return Geometry.getArea(((Way) p).getNodes()); 560 } 561 if (p instanceof Relation && p.isMultipolygon() && !p.isIncomplete()) { 562 Multipolygon mp = MultipolygonCache.getInstance().get((Relation) p); 563 if (mp.getOpenEnds().isEmpty()) { 564 Path2D path = new Path2D.Double(); 565 path.setWindingRule(Path2D.WIND_EVEN_ODD); 566 for (PolyData pd : mp.getCombinedPolygons()) { 567 path.append(pd.get(), false); 568 } 569 return new Area(path); 570 } 571 } 572 return new Area(); 573 } 574 575 /** 576 * Returns the Area of a polygon, from the multipolygon relation. 577 * @param multipolygon the multipolygon relation 578 * @return Area for the multipolygon (LatLon coordinates) 579 */ 580 public static Area getAreaLatLon(Relation multipolygon) { 581 final Multipolygon mp = MultipolygonCache.getInstance().get(multipolygon); 582 Path2D path = new Path2D.Double(); 583 path.setWindingRule(Path2D.WIND_EVEN_ODD); 584 for (Multipolygon.PolyData pd : mp.getCombinedPolygons()) { 585 buildPath2DLatLon(pd.getNodes(), path); 586 for (Multipolygon.PolyData pdInner : pd.getInners()) { 587 buildPath2DLatLon(pdInner.getNodes(), path); 588 } 589 } 590 return new Area(path); 591 } 592 593 /** 594 * Tests if two polygons intersect. 595 * @param first List of nodes forming first polygon 596 * @param second List of nodes forming second polygon 597 * @return intersection kind 598 */ 599 public static PolygonIntersection polygonIntersection(List<? extends INode> first, List<? extends INode> second) { 600 Area a1 = getArea(first); 601 Area a2 = getArea(second); 602 return polygonIntersection(a1, a2, INTERSECTION_EPS_EAST_NORTH); 603 } 604 605 /** 606 * Tests if two polygons intersect. It is assumed that the area is given in East North points. 607 * @param a1 Area of first polygon 608 * @param a2 Area of second polygon 609 * @return intersection kind 610 * @since 6841 611 */ 612 public static PolygonIntersection polygonIntersection(Area a1, Area a2) { 613 return polygonIntersection(a1, a2, INTERSECTION_EPS_EAST_NORTH); 614 } 615 616 /** 617 * Tests if two polygons intersect. 618 * @param a1 Area of first polygon 619 * @param a2 Area of second polygon 620 * @param eps an area threshold, everything below is considered an empty intersection 621 * @return intersection kind 622 */ 623 public static PolygonIntersection polygonIntersection(Area a1, Area a2, double eps) { 624 return polygonIntersectionResult(a1, a2, eps).a; 625 } 626 627 /** 628 * Calculate intersection area and kind of intersection between two polygons. 629 * @param a1 Area of first polygon 630 * @param a2 Area of second polygon 631 * @param eps an area threshold, everything below is considered an empty intersection 632 * @return pair with intersection kind and intersection area (never null, but maybe empty) 633 * @since 15938 634 */ 635 public static Pair<PolygonIntersection, Area> polygonIntersectionResult(Area a1, Area a2, double eps) { 636 Area inter = new Area(a1); 637 inter.intersect(a2); 638 639 if (inter.isEmpty() || !checkIntersection(inter, eps)) { 640 return new Pair<>(PolygonIntersection.OUTSIDE, inter); 641 } else if (a2.getBounds2D().contains(a1.getBounds2D()) && inter.equals(a1)) { 642 return new Pair<>(PolygonIntersection.FIRST_INSIDE_SECOND, inter); 643 } else if (a1.getBounds2D().contains(a2.getBounds2D()) && inter.equals(a2)) { 644 return new Pair<>(PolygonIntersection.SECOND_INSIDE_FIRST, inter); 645 } else { 646 return new Pair<>(PolygonIntersection.CROSSING, inter); 647 } 648 } 649 650 /** 651 * Check an intersection area which might describe multiple small polygons. 652 * Return true if any of the polygons is bigger than the given threshold. 653 * @param inter the intersection area 654 * @param eps an area threshold, everything below is considered an empty intersection 655 * @return true if any of the polygons is bigger than the given threshold 656 */ 657 private static boolean checkIntersection(Area inter, double eps) { 658 PathIterator pit = inter.getPathIterator(null); 659 double[] res = new double[6]; 660 Rectangle2D r = new Rectangle2D.Double(); 661 while (!pit.isDone()) { 662 int type = pit.currentSegment(res); 663 switch (type) { 664 case PathIterator.SEG_MOVETO: 665 r = new Rectangle2D.Double(res[0], res[1], 0, 0); 666 break; 667 case PathIterator.SEG_LINETO: 668 r.add(res[0], res[1]); 669 break; 670 case PathIterator.SEG_CLOSE: 671 if (r.getWidth() > eps || r.getHeight() > eps) 672 return true; 673 break; 674 default: 675 break; 676 } 677 pit.next(); 678 } 679 return false; 680 } 681 682 /** 683 * Tests if point is inside a polygon. The polygon can be self-intersecting. In such case the contains function works in xor-like manner. 684 * @param polygonNodes list of nodes from polygon path. 685 * @param point the point to test 686 * @return true if the point is inside polygon. 687 */ 688 public static boolean nodeInsidePolygon(INode point, List<? extends INode> polygonNodes) { 689 if (polygonNodes.size() < 2) 690 return false; 691 692 //iterate each side of the polygon, start with the last segment 693 INode oldPoint = polygonNodes.get(polygonNodes.size() - 1); 694 695 if (!oldPoint.isLatLonKnown()) { 696 return false; 697 } 698 699 boolean inside = false; 700 INode p1, p2; 701 702 for (INode newPoint : polygonNodes) { 703 //skip duplicate points 704 if (newPoint.equals(oldPoint)) { 705 continue; 706 } 707 708 if (!newPoint.isLatLonKnown()) { 709 return false; 710 } 711 712 //order points so p1.lat <= p2.lat 713 if (newPoint.getEastNorth().getY() > oldPoint.getEastNorth().getY()) { 714 p1 = oldPoint; 715 p2 = newPoint; 716 } else { 717 p1 = newPoint; 718 p2 = oldPoint; 719 } 720 721 EastNorth pEN = point.getEastNorth(); 722 EastNorth opEN = oldPoint.getEastNorth(); 723 EastNorth npEN = newPoint.getEastNorth(); 724 EastNorth p1EN = p1.getEastNorth(); 725 EastNorth p2EN = p2.getEastNorth(); 726 727 if (pEN != null && opEN != null && npEN != null && p1EN != null && p2EN != null) { 728 //test if the line is crossed and if so invert the inside flag. 729 if ((npEN.getY() < pEN.getY()) == (pEN.getY() <= opEN.getY()) 730 && (pEN.getX() - p1EN.getX()) * (p2EN.getY() - p1EN.getY()) 731 < (p2EN.getX() - p1EN.getX()) * (pEN.getY() - p1EN.getY())) { 732 inside = !inside; 733 } 734 } 735 736 oldPoint = newPoint; 737 } 738 739 return inside; 740 } 741 742 /** 743 * Returns area of a closed way in square meters. 744 * 745 * @param way Way to measure, should be closed (first node is the same as last node) 746 * @return area of the closed way. 747 */ 748 public static double closedWayArea(Way way) { 749 return getAreaAndPerimeter(way.getNodes(), Projections.getProjectionByCode("EPSG:54008")).getArea(); 750 } 751 752 /** 753 * Returns area of a multipolygon in square meters. 754 * 755 * @param multipolygon the multipolygon to measure 756 * @return area of the multipolygon. 757 */ 758 public static double multipolygonArea(Relation multipolygon) { 759 final Multipolygon mp = MultipolygonCache.getInstance().get(multipolygon); 760 return mp.getCombinedPolygons().stream() 761 .mapToDouble(pd -> pd.getAreaAndPerimeter(Projections.getProjectionByCode("EPSG:54008")).getArea()) 762 .sum(); 763 } 764 765 /** 766 * Computes the area of a closed way and multipolygon in square meters, or {@code null} for other primitives 767 * 768 * @param osm the primitive to measure 769 * @return area of the primitive, or {@code null} 770 * @since 13638 (signature) 771 */ 772 public static Double computeArea(IPrimitive osm) { 773 if (osm instanceof Way && ((Way) osm).isClosed()) { 774 return closedWayArea((Way) osm); 775 } else if (osm instanceof Relation && ((Relation) osm).isMultipolygon() && !((Relation) osm).hasIncompleteMembers()) { 776 return multipolygonArea((Relation) osm); 777 } else { 778 return null; 779 } 780 } 781 782 /** 783 * Determines whether a way is oriented clockwise. 784 * 785 * Internals: Assuming a closed non-looping way, compute twice the area 786 * of the polygon using the formula {@code 2 * area = sum (X[n] * Y[n+1] - X[n+1] * Y[n])}. 787 * If the area is negative the way is ordered in a clockwise direction. 788 * 789 * See http://paulbourke.net/geometry/polyarea/ 790 * 791 * @param w the way to be checked. 792 * @return true if and only if way is oriented clockwise. 793 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}). 794 */ 795 public static boolean isClockwise(Way w) { 796 return isClockwise(w.getNodes()); 797 } 798 799 /** 800 * Determines whether path from nodes list is oriented clockwise. 801 * @param nodes Nodes list to be checked. 802 * @return true if and only if way is oriented clockwise. 803 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}). 804 * @see #isClockwise(Way) 805 */ 806 public static boolean isClockwise(List<? extends INode> nodes) { 807 int nodesCount = nodes.size(); 808 if (nodesCount < 3 || nodes.get(0) != nodes.get(nodesCount - 1)) { 809 throw new IllegalArgumentException("Way must be closed to check orientation."); 810 } 811 double area2 = 0.; 812 813 for (int node = 1; node <= /*sic! consider last-first as well*/ nodesCount; node++) { 814 INode coorPrev = nodes.get(node - 1); 815 INode coorCurr = nodes.get(node % nodesCount); 816 area2 += coorPrev.lon() * coorCurr.lat(); 817 area2 -= coorCurr.lon() * coorPrev.lat(); 818 } 819 return area2 < 0; 820 } 821 822 /** 823 * Returns angle of a segment defined with 2 point coordinates. 824 * 825 * @param p1 first point 826 * @param p2 second point 827 * @return Angle in radians (-pi, pi] 828 */ 829 public static double getSegmentAngle(EastNorth p1, EastNorth p2) { 830 831 CheckParameterUtil.ensureThat(p1.isValid(), () -> p1 + " invalid"); 832 CheckParameterUtil.ensureThat(p2.isValid(), () -> p2 + " invalid"); 833 834 return Math.atan2(p2.north() - p1.north(), p2.east() - p1.east()); 835 } 836 837 /** 838 * Returns angle of a corner defined with 3 point coordinates. 839 * 840 * @param p1 first point 841 * @param common Common end point 842 * @param p3 third point 843 * @return Angle in radians (-pi, pi] 844 */ 845 public static double getCornerAngle(EastNorth p1, EastNorth common, EastNorth p3) { 846 847 CheckParameterUtil.ensureThat(p1.isValid(), () -> p1 + " invalid"); 848 CheckParameterUtil.ensureThat(common.isValid(), () -> common + " invalid"); 849 CheckParameterUtil.ensureThat(p3.isValid(), () -> p3 + " invalid"); 850 851 double result = getSegmentAngle(common, p1) - getSegmentAngle(common, p3); 852 if (result <= -Math.PI) { 853 result += 2 * Math.PI; 854 } 855 856 if (result > Math.PI) { 857 result -= 2 * Math.PI; 858 } 859 860 return result; 861 } 862 863 /** 864 * Get angles in radians and return it's value in range [0, 180]. 865 * 866 * @param angle the angle in radians 867 * @return normalized angle in degrees 868 * @since 13670 869 */ 870 public static double getNormalizedAngleInDegrees(double angle) { 871 return Math.abs(180 * angle / Math.PI); 872 } 873 874 /** 875 * Compute the centroid/barycenter of nodes 876 * @param nodes Nodes for which the centroid is wanted 877 * @return the centroid of nodes 878 * @see Geometry#getCenter 879 */ 880 public static EastNorth getCentroid(List<? extends INode> nodes) { 881 return getCentroidEN(nodes.stream().map(INode::getEastNorth).collect(Collectors.toList())); 882 } 883 884 /** 885 * Compute the centroid/barycenter of nodes 886 * @param nodes Coordinates for which the centroid is wanted 887 * @return the centroid of nodes 888 * @since 13712 889 */ 890 public static EastNorth getCentroidEN(List<EastNorth> nodes) { 891 892 final int size = nodes.size(); 893 if (size == 1) { 894 return nodes.get(0); 895 } else if (size == 2) { 896 return nodes.get(0).getCenter(nodes.get(1)); 897 } 898 899 BigDecimal area = BigDecimal.ZERO; 900 BigDecimal north = BigDecimal.ZERO; 901 BigDecimal east = BigDecimal.ZERO; 902 903 // See https://en.wikipedia.org/wiki/Centroid#Of_a_polygon for the equation used here 904 for (int i = 0; i < size; i++) { 905 EastNorth n0 = nodes.get(i); 906 EastNorth n1 = nodes.get((i+1) % size); 907 908 if (n0 != null && n1 != null && n0.isValid() && n1.isValid()) { 909 BigDecimal x0 = BigDecimal.valueOf(n0.east()); 910 BigDecimal y0 = BigDecimal.valueOf(n0.north()); 911 BigDecimal x1 = BigDecimal.valueOf(n1.east()); 912 BigDecimal y1 = BigDecimal.valueOf(n1.north()); 913 914 BigDecimal k = x0.multiply(y1, MathContext.DECIMAL128).subtract(y0.multiply(x1, MathContext.DECIMAL128)); 915 916 area = area.add(k, MathContext.DECIMAL128); 917 east = east.add(k.multiply(x0.add(x1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 918 north = north.add(k.multiply(y0.add(y1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 919 } 920 } 921 922 BigDecimal d = new BigDecimal(3, MathContext.DECIMAL128); // 1/2 * 6 = 3 923 area = area.multiply(d, MathContext.DECIMAL128); 924 if (area.compareTo(BigDecimal.ZERO) != 0) { 925 north = north.divide(area, MathContext.DECIMAL128); 926 east = east.divide(area, MathContext.DECIMAL128); 927 } 928 929 return new EastNorth(east.doubleValue(), north.doubleValue()); 930 } 931 932 /** 933 * Compute center of the circle closest to different nodes. 934 * 935 * Ensure exact center computation in case nodes are already aligned in circle. 936 * This is done by least square method. 937 * Let be a_i x + b_i y + c_i = 0 equations of bisectors of each edges. 938 * Center must be intersection of all bisectors. 939 * <pre> 940 * [ a1 b1 ] [ -c1 ] 941 * With A = [ ... ... ] and Y = [ ... ] 942 * [ an bn ] [ -cn ] 943 * </pre> 944 * An approximation of center of circle is (At.A)^-1.At.Y 945 * @param nodes Nodes parts of the circle (at least 3) 946 * @return An approximation of the center, of null if there is no solution. 947 * @see Geometry#getCentroid 948 * @since 6934 949 */ 950 public static EastNorth getCenter(List<? extends INode> nodes) { 951 int nc = nodes.size(); 952 if (nc < 3) return null; 953 /** 954 * Equation of each bisector ax + by + c = 0 955 */ 956 double[] a = new double[nc]; 957 double[] b = new double[nc]; 958 double[] c = new double[nc]; 959 // Compute equation of bisector 960 for (int i = 0; i < nc; i++) { 961 EastNorth pt1 = nodes.get(i).getEastNorth(); 962 EastNorth pt2 = nodes.get((i+1) % nc).getEastNorth(); 963 a[i] = pt1.east() - pt2.east(); 964 b[i] = pt1.north() - pt2.north(); 965 double d = Math.sqrt(a[i]*a[i] + b[i]*b[i]); 966 if (d == 0) return null; 967 a[i] /= d; 968 b[i] /= d; 969 double xC = (pt1.east() + pt2.east()) / 2; 970 double yC = (pt1.north() + pt2.north()) / 2; 971 c[i] = -(a[i]*xC + b[i]*yC); 972 } 973 // At.A = [aij] 974 double a11 = 0, a12 = 0, a22 = 0; 975 // At.Y = [bi] 976 double b1 = 0, b2 = 0; 977 for (int i = 0; i < nc; i++) { 978 a11 += a[i]*a[i]; 979 a12 += a[i]*b[i]; 980 a22 += b[i]*b[i]; 981 b1 -= a[i]*c[i]; 982 b2 -= b[i]*c[i]; 983 } 984 // (At.A)^-1 = [invij] 985 double det = a11*a22 - a12*a12; 986 if (Math.abs(det) < 1e-5) return null; 987 double inv11 = a22/det; 988 double inv12 = -a12/det; 989 double inv22 = a11/det; 990 // center (xC, yC) = (At.A)^-1.At.y 991 double xC = inv11*b1 + inv12*b2; 992 double yC = inv12*b1 + inv22*b2; 993 return new EastNorth(xC, yC); 994 } 995 996 /** 997 * Tests if the {@code node} is inside the multipolygon {@code multiPolygon}. The nullable argument 998 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match. 999 * For repeated tests against {@code multiPolygon} better use {@link Geometry#filterInsideMultipolygon}. 1000 * @param node node 1001 * @param multiPolygon multipolygon 1002 * @param isOuterWayAMatch allows to decide if the immediate {@code outer} way of the multipolygon is a match 1003 * @return {@code true} if the node is inside the multipolygon 1004 */ 1005 public static boolean isNodeInsideMultiPolygon(INode node, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) { 1006 return isPolygonInsideMultiPolygon(Collections.singletonList(node), multiPolygon, isOuterWayAMatch); 1007 } 1008 1009 /** 1010 * Tests if the polygon formed by {@code nodes} is inside the multipolygon {@code multiPolygon}. The nullable argument 1011 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match. 1012 * For repeated tests against {@code multiPolygon} better use {@link Geometry#filterInsideMultipolygon}. 1013 * <p> 1014 * If {@code nodes} contains exactly one element, then it is checked whether that one node is inside the multipolygon. 1015 * @param nodes nodes forming the polygon 1016 * @param multiPolygon multipolygon 1017 * @param isOuterWayAMatch allows to decide if the immediate {@code outer} way of the multipolygon is a match 1018 * @return {@code true} if the multipolygon is valid and the polygon formed by nodes is inside the multipolygon 1019 */ 1020 public static boolean isPolygonInsideMultiPolygon(List<? extends INode> nodes, Relation multiPolygon, Predicate<Way> isOuterWayAMatch) { 1021 try { 1022 return isPolygonInsideMultiPolygon(nodes, MultipolygonBuilder.joinWays(multiPolygon), isOuterWayAMatch); 1023 } catch (MultipolygonBuilder.JoinedPolygonCreationException ex) { 1024 Logging.trace(ex); 1025 Logging.debug("Invalid multipolygon " + multiPolygon); 1026 return false; 1027 } 1028 } 1029 1030 /** 1031 * Tests if the polygon formed by {@code nodes} is inside the multipolygon {@code multiPolygon}. The nullable argument 1032 * {@code isOuterWayAMatch} allows to decide if the immediate {@code outer} way of the multipolygon is a match. 1033 * For repeated tests against {@code multiPolygon} better use {@link Geometry#filterInsideMultipolygon}. 1034 * <p> 1035 * If {@code nodes} contains exactly one element, then it is checked whether that one node is inside the multipolygon. 1036 * @param nodes nodes forming the polygon 1037 * @param outerInner result of {@link MultipolygonBuilder#joinWays(Relation)} 1038 * @param isOuterWayAMatch allows to decide if the immediate {@code outer} way of the multipolygon is a match 1039 * @return {@code true} if the multipolygon is valid and the polygon formed by nodes is inside the multipolygon 1040 * @since 15069 1041 */ 1042 public static boolean isPolygonInsideMultiPolygon(List<? extends INode> nodes, Pair<List<JoinedPolygon>, 1043 List<JoinedPolygon>> outerInner, Predicate<Way> isOuterWayAMatch) { 1044 Area a1 = nodes.size() == 1 ? null : getArea(nodes); 1045 // Test if object is inside an outer member 1046 for (JoinedPolygon out : outerInner.a) { 1047 if (a1 == null 1048 ? nodeInsidePolygon(nodes.get(0), out.nodes) 1049 : PolygonIntersection.FIRST_INSIDE_SECOND == polygonIntersection(a1, out.area)) { 1050 // If inside an outer, check it is not inside an inner 1051 boolean insideInner = outerInner.b.stream().anyMatch(in -> a1 == null 1052 ? nodeInsidePolygon(nodes.get(0), in.nodes) 1053 : in.area.getBounds2D().contains(a1.getBounds2D()) 1054 && polygonIntersection(a1, in.area) == PolygonIntersection.FIRST_INSIDE_SECOND 1055 && polygonIntersection(in.area, out.area) == PolygonIntersection.FIRST_INSIDE_SECOND); 1056 if (!insideInner) { 1057 // Final check using predicate 1058 if (isOuterWayAMatch == null || isOuterWayAMatch.test(out.ways.get(0) 1059 /* TODO give a better representation of the outer ring to the predicate */)) { 1060 return true; 1061 } 1062 } 1063 } 1064 } 1065 return false; 1066 } 1067 1068 /** 1069 * Find all primitives in the given collection which are inside the given polygon. 1070 * 1071 * @param primitives the primitives 1072 * @param polygon the closed way or multipolygon relation 1073 * @return a new list containing the found primitives, empty if polygon is invalid or nothing was found. 1074 * @see Geometry#filterInsidePolygon 1075 * @see Geometry#filterInsideMultipolygon 1076 * @since 15730 1077 */ 1078 public static List<IPrimitive> filterInsideAnyPolygon(Collection<IPrimitive> primitives, IPrimitive polygon) { 1079 if (polygon instanceof IWay<?>) { 1080 return filterInsidePolygon(primitives, (IWay<?>) polygon); 1081 } else if (polygon instanceof Relation && polygon.isMultipolygon()) { 1082 return filterInsideMultipolygon(primitives, (Relation) polygon); 1083 } 1084 return Collections.emptyList(); 1085 } 1086 1087 /** 1088 * Find all primitives in the given collection which are inside the given polygon. 1089 * Unclosed ways and multipolygon relations with unclosed outer rings are ignored. 1090 * 1091 * @param primitives the primitives 1092 * @param polygon the polygon 1093 * @return a new list containing the found primitives, empty if polygon is invalid or nothing was found. 1094 * @since 15069 (for {@link List} of {@code primitives}, 15730 for a {@link Collection} of {@code primitives}) 1095 */ 1096 public static List<IPrimitive> filterInsidePolygon(Collection<IPrimitive> primitives, IWay<?> polygon) { 1097 List<IPrimitive> res = new ArrayList<>(); 1098 if (!polygon.isClosed() || polygon.getNodesCount() <= 3) 1099 return res; 1100 /** polygon area in east north space, calculated only when really needed */ 1101 Area polygonArea = null; 1102 for (IPrimitive p : primitives) { 1103 if (p instanceof INode) { 1104 if (nodeInsidePolygon((INode) p, polygon.getNodes())) { 1105 res.add(p); 1106 } 1107 } else if (p instanceof IWay) { 1108 if (((IWay<?>) p).isClosed()) { 1109 if (polygonArea == null) { 1110 polygonArea = getArea(polygon.getNodes()); 1111 } 1112 if (PolygonIntersection.FIRST_INSIDE_SECOND == polygonIntersection(getArea(((IWay<?>) p).getNodes()), 1113 polygonArea)) { 1114 res.add(p); 1115 } 1116 } 1117 } else if (p.isMultipolygon()) { 1118 if (polygonArea == null) { 1119 polygonArea = getArea(polygon.getNodes()); 1120 } 1121 Multipolygon mp = new Multipolygon((Relation) p); 1122 boolean inside = true; 1123 // a (valid) multipolygon is inside the polygon if all outer rings are inside 1124 for (PolyData outer : mp.getOuterPolygons()) { 1125 if (!outer.isClosed() 1126 || PolygonIntersection.FIRST_INSIDE_SECOND != polygonIntersection(getArea(outer.getNodes()), 1127 polygonArea)) { 1128 inside = false; 1129 break; 1130 } 1131 } 1132 if (inside) { 1133 res.add(p); 1134 } 1135 } 1136 } 1137 return res; 1138 } 1139 1140 /** 1141 * Find all primitives in the given collection which are inside the given multipolygon. Members of the multipolygon are 1142 * ignored. Unclosed ways and multipolygon relations with unclosed outer rings are ignored. 1143 * @param primitives the primitives 1144 * @param multiPolygon the multipolygon relation 1145 * @return a new list containing the found primitives, empty if multipolygon is invalid or nothing was found. 1146 * @since 15069 1147 */ 1148 public static List<IPrimitive> filterInsideMultipolygon(Collection<IPrimitive> primitives, Relation multiPolygon) { 1149 List<IPrimitive> res = new ArrayList<>(); 1150 if (primitives.isEmpty()) 1151 return res; 1152 1153 final Pair<List<JoinedPolygon>, List<JoinedPolygon>> outerInner; 1154 try { 1155 outerInner = MultipolygonBuilder.joinWays(multiPolygon); 1156 } catch (MultipolygonBuilder.JoinedPolygonCreationException ex) { 1157 Logging.trace(ex); 1158 Logging.debug("Invalid multipolygon " + multiPolygon); 1159 return res; 1160 } 1161 1162 Set<OsmPrimitive> members = multiPolygon.getMemberPrimitives(); 1163 for (IPrimitive p : primitives) { 1164 if (members.contains(p)) 1165 continue; 1166 if (p instanceof Node) { 1167 if (isPolygonInsideMultiPolygon(Collections.singletonList((Node) p), outerInner, null)) { 1168 res.add(p); 1169 } 1170 } else if (p instanceof Way) { 1171 if (((IWay<?>) p).isClosed() && isPolygonInsideMultiPolygon(((Way) p).getNodes(), outerInner, null)) { 1172 res.add(p); 1173 } 1174 } else if (p.isMultipolygon()) { 1175 Multipolygon mp = new Multipolygon((Relation) p); 1176 // a (valid) multipolygon is inside multiPolygon if all outer rings are inside 1177 boolean inside = mp.getOuterPolygons().stream() 1178 .allMatch(outer -> outer.isClosed() && isPolygonInsideMultiPolygon(outer.getNodes(), outerInner, null)); 1179 if (inside) { 1180 res.add(p); 1181 } 1182 } 1183 } 1184 return res; 1185 } 1186 1187 /** 1188 * Data class to hold two double values (area and perimeter of a polygon). 1189 */ 1190 public static class AreaAndPerimeter { 1191 private final double area; 1192 private final double perimeter; 1193 1194 /** 1195 * Create a new {@link AreaAndPerimeter} 1196 * @param area The area 1197 * @param perimeter The perimeter 1198 */ 1199 public AreaAndPerimeter(double area, double perimeter) { 1200 this.area = area; 1201 this.perimeter = perimeter; 1202 } 1203 1204 /** 1205 * Gets the area 1206 * @return The area size 1207 */ 1208 public double getArea() { 1209 return area; 1210 } 1211 1212 /** 1213 * Gets the perimeter 1214 * @return The perimeter length 1215 */ 1216 public double getPerimeter() { 1217 return perimeter; 1218 } 1219 } 1220 1221 /** 1222 * Calculate area and perimeter length of a polygon. 1223 * 1224 * Uses current projection; units are that of the projected coordinates. 1225 * 1226 * @param nodes the list of nodes representing the polygon 1227 * @return area and perimeter 1228 */ 1229 public static AreaAndPerimeter getAreaAndPerimeter(List<? extends ILatLon> nodes) { 1230 return getAreaAndPerimeter(nodes, null); 1231 } 1232 1233 /** 1234 * Calculate area and perimeter length of a polygon in the given projection. 1235 * 1236 * @param nodes the list of nodes representing the polygon 1237 * @param projection the projection to use for the calculation, {@code null} defaults to {@link ProjectionRegistry#getProjection()} 1238 * @return area and perimeter 1239 * @since 13638 (signature) 1240 */ 1241 public static AreaAndPerimeter getAreaAndPerimeter(List<? extends ILatLon> nodes, Projection projection) { 1242 CheckParameterUtil.ensureParameterNotNull(nodes, "nodes"); 1243 double area = 0; 1244 double perimeter = 0; 1245 Projection useProjection = projection == null ? ProjectionRegistry.getProjection() : projection; 1246 1247 if (!nodes.isEmpty()) { 1248 boolean closed = nodes.get(0) == nodes.get(nodes.size() - 1); 1249 int numSegments = closed ? nodes.size() - 1 : nodes.size(); 1250 EastNorth p1 = nodes.get(0).getEastNorth(useProjection); 1251 for (int i = 1; i <= numSegments; i++) { 1252 final ILatLon node = nodes.get(i == numSegments ? 0 : i); 1253 final EastNorth p2 = node.getEastNorth(useProjection); 1254 if (p1 != null && p2 != null) { 1255 area += p1.east() * p2.north() - p2.east() * p1.north(); 1256 perimeter += p1.distance(p2); 1257 } 1258 p1 = p2; 1259 } 1260 } 1261 return new AreaAndPerimeter(Math.abs(area) / 2, perimeter); 1262 } 1263 1264 /** 1265 * Get the closest primitive to {@code osm} from the collection of 1266 * OsmPrimitive {@code primitives} 1267 * 1268 * The {@code primitives} should be fully downloaded to ensure accuracy. 1269 * 1270 * Note: The complexity of this method is O(n*m), where n is the number of 1271 * children {@code osm} has plus 1, m is the number of children the 1272 * collection of primitives have plus the number of primitives in the 1273 * collection. 1274 * 1275 * @param <T> The return type of the primitive 1276 * @param osm The primitive to get the distances from 1277 * @param primitives The collection of primitives to get the distance to 1278 * @return The closest {@link OsmPrimitive}. This is not determinative. 1279 * To get all primitives that share the same distance, use 1280 * {@link Geometry#getClosestPrimitives}. 1281 * @since 15035 1282 */ 1283 public static <T extends OsmPrimitive> T getClosestPrimitive(OsmPrimitive osm, Collection<T> primitives) { 1284 Collection<T> collection = getClosestPrimitives(osm, primitives); 1285 return collection.iterator().next(); 1286 } 1287 1288 /** 1289 * Get the closest primitives to {@code osm} from the collection of 1290 * OsmPrimitive {@code primitives} 1291 * 1292 * The {@code primitives} should be fully downloaded to ensure accuracy. 1293 * 1294 * Note: The complexity of this method is O(n*m), where n is the number of 1295 * children {@code osm} has plus 1, m is the number of children the 1296 * collection of primitives have plus the number of primitives in the 1297 * collection. 1298 * 1299 * @param <T> The return type of the primitive 1300 * @param osm The primitive to get the distances from 1301 * @param primitives The collection of primitives to get the distance to 1302 * @return The closest {@link OsmPrimitive}s. May be empty. 1303 * @since 15035 1304 */ 1305 public static <T extends OsmPrimitive> Collection<T> getClosestPrimitives(OsmPrimitive osm, Collection<T> primitives) { 1306 double lowestDistance = Double.MAX_VALUE; 1307 TreeSet<T> closest = new TreeSet<>(); 1308 for (T primitive : primitives) { 1309 double distance = getDistance(osm, primitive); 1310 if (Double.isNaN(distance)) continue; 1311 int comp = Double.compare(distance, lowestDistance); 1312 if (comp < 0) { 1313 closest.clear(); 1314 lowestDistance = distance; 1315 closest.add(primitive); 1316 } else if (comp == 0) { 1317 closest.add(primitive); 1318 } 1319 } 1320 return closest; 1321 } 1322 1323 /** 1324 * Get the furthest primitive to {@code osm} from the collection of 1325 * OsmPrimitive {@code primitives} 1326 * 1327 * The {@code primitives} should be fully downloaded to ensure accuracy. 1328 * 1329 * It does NOT give the furthest primitive based off of the furthest 1330 * part of that primitive 1331 * 1332 * Note: The complexity of this method is O(n*m), where n is the number of 1333 * children {@code osm} has plus 1, m is the number of children the 1334 * collection of primitives have plus the number of primitives in the 1335 * collection. 1336 * 1337 * @param <T> The return type of the primitive 1338 * @param osm The primitive to get the distances from 1339 * @param primitives The collection of primitives to get the distance to 1340 * @return The furthest {@link OsmPrimitive}. This is not determinative. 1341 * To get all primitives that share the same distance, use 1342 * {@link Geometry#getFurthestPrimitives} 1343 * @since 15035 1344 */ 1345 public static <T extends OsmPrimitive> T getFurthestPrimitive(OsmPrimitive osm, Collection<T> primitives) { 1346 return getFurthestPrimitives(osm, primitives).iterator().next(); 1347 } 1348 1349 /** 1350 * Get the furthest primitives to {@code osm} from the collection of 1351 * OsmPrimitive {@code primitives} 1352 * 1353 * The {@code primitives} should be fully downloaded to ensure accuracy. 1354 * 1355 * It does NOT give the furthest primitive based off of the furthest 1356 * part of that primitive 1357 * 1358 * Note: The complexity of this method is O(n*m), where n is the number of 1359 * children {@code osm} has plus 1, m is the number of children the 1360 * collection of primitives have plus the number of primitives in the 1361 * collection. 1362 * 1363 * @param <T> The return type of the primitive 1364 * @param osm The primitive to get the distances from 1365 * @param primitives The collection of primitives to get the distance to 1366 * @return The furthest {@link OsmPrimitive}s. It may return an empty collection. 1367 * @since 15035 1368 */ 1369 public static <T extends OsmPrimitive> Collection<T> getFurthestPrimitives(OsmPrimitive osm, Collection<T> primitives) { 1370 double furthestDistance = Double.NEGATIVE_INFINITY; 1371 TreeSet<T> furthest = new TreeSet<>(); 1372 for (T primitive : primitives) { 1373 double distance = getDistance(osm, primitive); 1374 if (Double.isNaN(distance)) continue; 1375 int comp = Double.compare(distance, furthestDistance); 1376 if (comp > 0) { 1377 furthest.clear(); 1378 furthestDistance = distance; 1379 furthest.add(primitive); 1380 } else if (comp == 0) { 1381 furthest.add(primitive); 1382 } 1383 } 1384 return furthest; 1385 } 1386 1387 /** 1388 * Get the distance between different {@link OsmPrimitive}s 1389 * @param one The primitive to get the distance from 1390 * @param two The primitive to get the distance to 1391 * @return The distance between the primitives in meters 1392 * (or the unit of the current projection, see {@link Projection}). 1393 * May return {@link Double#NaN} if one of the primitives is incomplete. 1394 * 1395 * Note: The complexity is O(n*m), where (n,m) are the number of child 1396 * objects the {@link OsmPrimitive}s have. 1397 * @since 15035 1398 */ 1399 public static double getDistance(OsmPrimitive one, OsmPrimitive two) { 1400 double rValue = Double.MAX_VALUE; 1401 if (one == null || two == null || one.isIncomplete() 1402 || two.isIncomplete()) return Double.NaN; 1403 if (one instanceof Node && two instanceof Node) { 1404 rValue = ((Node) one).getCoor().greatCircleDistance(((Node) two).getCoor()); 1405 } else if (one instanceof Node && two instanceof Way) { 1406 rValue = getDistanceWayNode((Way) two, (Node) one); 1407 } else if (one instanceof Way && two instanceof Node) { 1408 rValue = getDistanceWayNode((Way) one, (Node) two); 1409 } else if (one instanceof Way && two instanceof Way) { 1410 rValue = getDistanceWayWay((Way) one, (Way) two); 1411 } else if (one instanceof Relation && !(two instanceof Relation)) { 1412 for (OsmPrimitive osmPrimitive: ((Relation) one).getMemberPrimitives()) { 1413 double currentDistance = getDistance(osmPrimitive, two); 1414 if (currentDistance < rValue) rValue = currentDistance; 1415 } 1416 } else if (!(one instanceof Relation) && two instanceof Relation) { 1417 rValue = getDistance(two, one); 1418 } else if (one instanceof Relation && two instanceof Relation) { 1419 for (OsmPrimitive osmPrimitive1 : ((Relation) one).getMemberPrimitives()) { 1420 for (OsmPrimitive osmPrimitive2 : ((Relation) two).getMemberPrimitives()) { 1421 double currentDistance = getDistance(osmPrimitive1, osmPrimitive2); 1422 if (currentDistance < rValue) rValue = currentDistance; 1423 } 1424 } 1425 } 1426 return rValue != Double.MAX_VALUE ? rValue : Double.NaN; 1427 } 1428 1429 /** 1430 * Get the distance between a way and a node 1431 * @param way The way to get the distance from 1432 * @param node The node to get the distance to 1433 * @return The distance between the {@code way} and the {@code node} in 1434 * meters (or the unit of the current projection, see {@link Projection}). 1435 * May return {@link Double#NaN} if the primitives are incomplete. 1436 * @since 15035 1437 */ 1438 public static double getDistanceWayNode(Way way, Node node) { 1439 if (way == null || node == null || way.getNodesCount() < 2 || !node.isLatLonKnown()) 1440 return Double.NaN; 1441 1442 double smallest = Double.MAX_VALUE; 1443 EastNorth en0 = node.getEastNorth(); 1444 // go through the nodes as if they were paired 1445 Iterator<Node> iter = way.getNodes().iterator(); 1446 EastNorth en1 = iter.next().getEastNorth(); 1447 while (iter.hasNext()) { 1448 EastNorth en2 = iter.next().getEastNorth(); 1449 double distance = getSegmentNodeDistSq(en1, en2, en0); 1450 if (distance < smallest) 1451 smallest = distance; 1452 en1 = en2; 1453 } 1454 return smallest != Double.MAX_VALUE ? Math.sqrt(smallest) : Double.NaN; 1455 } 1456 1457 /** 1458 * Get the closest {@link WaySegment} from a way to a primitive. 1459 * @param way The {@link Way} to get the distance from and the {@link WaySegment} 1460 * @param primitive The {@link OsmPrimitive} to get the distance to 1461 * @return The {@link WaySegment} that is closest to {@code primitive} from {@code way}. 1462 * If there are multiple {@link WaySegment}s with the same distance, the last 1463 * {@link WaySegment} with the same distance will be returned. 1464 * May return {@code null} if the way has fewer than two nodes or one 1465 * of the primitives is incomplete. 1466 * @since 15035 1467 */ 1468 public static WaySegment getClosestWaySegment(Way way, OsmPrimitive primitive) { 1469 if (way == null || primitive == null || way.isIncomplete() 1470 || primitive.isIncomplete()) return null; 1471 double lowestDistance = Double.MAX_VALUE; 1472 Pair<Node, Node> closestNodes = null; 1473 for (Pair<Node, Node> nodes : way.getNodePairs(false)) { 1474 Way tWay = new Way(); 1475 tWay.addNode(nodes.a); 1476 tWay.addNode(nodes.b); 1477 double distance = getDistance(tWay, primitive); 1478 if (distance < lowestDistance) { 1479 lowestDistance = distance; 1480 closestNodes = nodes; 1481 } 1482 } 1483 if (closestNodes == null) return null; 1484 return lowestDistance != Double.MAX_VALUE ? WaySegment.forNodePair(way, closestNodes.a, closestNodes.b) : null; 1485 } 1486 1487 /** 1488 * Get the distance between different ways. Iterates over the nodes of the ways, complexity is O(n*m) 1489 * (n,m giving the number of nodes) 1490 * @param w1 The first {@link Way} 1491 * @param w2 The second {@link Way} 1492 * @return The shortest distance between the ways in meters 1493 * (or the unit of the current projection, see {@link Projection}). 1494 * May return {@link Double#NaN}. 1495 * @since 15035 1496 */ 1497 public static double getDistanceWayWay(Way w1, Way w2) { 1498 if (w1 == null || w2 == null || w1.getNodesCount() < 2 || w2.getNodesCount() < 2) 1499 return Double.NaN; 1500 double rValue = Double.MAX_VALUE; 1501 Iterator<Node> iter1 = w1.getNodes().iterator(); 1502 Node w1N1 = iter1.next(); 1503 while (iter1.hasNext()) { 1504 Node w1N2 = iter1.next(); 1505 Iterator<Node> iter2 = w2.getNodes().iterator(); 1506 Node w2N1 = iter2.next(); 1507 while (iter2.hasNext()) { 1508 Node w2N2 = iter2.next(); 1509 double distance = getDistanceSegmentSegment(w1N1, w1N2, w2N1, w2N2); 1510 if (distance < rValue) 1511 rValue = distance; 1512 w2N1 = w2N2; 1513 } 1514 w1N1 = w1N2; 1515 } 1516 return rValue != Double.MAX_VALUE ? rValue : Double.NaN; 1517 } 1518 1519 /** 1520 * Get the distance between different {@link WaySegment}s 1521 * @param ws1 A {@link WaySegment} 1522 * @param ws2 A {@link WaySegment} 1523 * @return The distance between the two {@link WaySegment}s in meters 1524 * (or the unit of the current projection, see {@link Projection}). 1525 * May return {@link Double#NaN}. 1526 * @since 15035 1527 */ 1528 public static double getDistanceSegmentSegment(WaySegment ws1, WaySegment ws2) { 1529 return getDistanceSegmentSegment(ws1.getFirstNode(), ws1.getSecondNode(), ws2.getFirstNode(), ws2.getSecondNode()); 1530 } 1531 1532 /** 1533 * Get the distance between different {@link WaySegment}s 1534 * @param ws1Node1 The first node of the first WaySegment 1535 * @param ws1Node2 The second node of the second WaySegment 1536 * @param ws2Node1 The first node of the second WaySegment 1537 * @param ws2Node2 The second node of the second WaySegment 1538 * @return The distance between the two {@link WaySegment}s in meters 1539 * (or the unit of the current projection, see {@link Projection}). 1540 * May return {@link Double#NaN}. 1541 * @since 15035 1542 */ 1543 public static double getDistanceSegmentSegment(Node ws1Node1, Node ws1Node2, Node ws2Node1, Node ws2Node2) { 1544 EastNorth enWs1Node1 = ws1Node1.getEastNorth(); 1545 EastNorth enWs1Node2 = ws1Node2.getEastNorth(); 1546 EastNorth enWs2Node1 = ws2Node1.getEastNorth(); 1547 EastNorth enWs2Node2 = ws2Node2.getEastNorth(); 1548 if (enWs1Node1 == null || enWs1Node2 == null || enWs2Node1 == null || enWs2Node2 == null) 1549 return Double.NaN; 1550 if (getSegmentSegmentIntersection(enWs1Node1, enWs1Node2, enWs2Node1, enWs2Node2) != null) 1551 return 0; 1552 1553 double dist1sq = getSegmentNodeDistSq(enWs1Node1, enWs1Node2, enWs2Node1); 1554 double dist2sq = getSegmentNodeDistSq(enWs1Node1, enWs1Node2, enWs2Node2); 1555 double dist3sq = getSegmentNodeDistSq(enWs2Node1, enWs2Node2, enWs1Node1); 1556 double dist4sq = getSegmentNodeDistSq(enWs2Node1, enWs2Node2, enWs1Node2); 1557 double smallest = Math.min(Math.min(dist1sq, dist2sq), Math.min(dist3sq, dist4sq)); 1558 return smallest != Double.MAX_VALUE ? Math.sqrt(smallest) : Double.NaN; 1559 } 1560 1561 /** 1562 * Create a new LatLon at a specified distance. Currently uses WGS84, but may change randomly in the future. 1563 * This does not currently attempt to be hugely accurate. The actual location may be off 1564 * depending upon the distance and the elevation, but should be within 0.0002 meters. 1565 * 1566 * @param original The originating point 1567 * @param angle The angle (from true north) in radians 1568 * @param offset The distance to the new point in the current projection's units 1569 * @return The location at the specified angle and distance from the originating point 1570 * @since 18109 1571 */ 1572 public static ILatLon getLatLonFrom(final ILatLon original, final double angle, final double offset) { 1573 final double meterOffset = ProjectionRegistry.getProjection().getMetersPerUnit() * offset; 1574 final double radianLat = Math.toRadians(original.lat()); 1575 final double radianLon = Math.toRadians(original.lon()); 1576 final double angularDistance = meterOffset / WGS84.a; 1577 final double lat = Math.asin(Math.sin(radianLat) * Math.cos(angularDistance) 1578 + Math.cos(radianLat) * Math.sin(angularDistance) * Math.cos(angle)); 1579 final double lon = radianLon + Math.atan2(Math.sin(angle) * Math.sin(angularDistance) * Math.cos(radianLat), 1580 Math.cos(angularDistance) - Math.sin(radianLat) * Math.sin(lat)); 1581 return new LatLon(Math.toDegrees(lat), Math.toDegrees(lon)); 1582 } 1583 1584 /** 1585 * Calculate closest distance between a line segment s1-s2 and a point p 1586 * @param s1 start of segment 1587 * @param s2 end of segment 1588 * @param p the point 1589 * @return the square of the euclidean distance from p to the closest point on the segment 1590 */ 1591 private static double getSegmentNodeDistSq(EastNorth s1, EastNorth s2, EastNorth p) { 1592 EastNorth c1 = closestPointTo(s1, s2, p, true); 1593 return c1.distanceSq(p); 1594 } 1595}