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