001// License: GPL. For details, see LICENSE file. 002package org.openstreetmap.josm.tools; 003 004import java.awt.Rectangle; 005import java.awt.geom.Area; 006import java.awt.geom.Line2D; 007import java.awt.geom.Path2D; 008import java.math.BigDecimal; 009import java.math.MathContext; 010import java.util.ArrayList; 011import java.util.Comparator; 012import java.util.LinkedHashSet; 013import java.util.List; 014import java.util.Set; 015 016import org.openstreetmap.josm.Main; 017import org.openstreetmap.josm.command.AddCommand; 018import org.openstreetmap.josm.command.ChangeCommand; 019import org.openstreetmap.josm.command.Command; 020import org.openstreetmap.josm.data.coor.EastNorth; 021import org.openstreetmap.josm.data.coor.LatLon; 022import org.openstreetmap.josm.data.osm.BBox; 023import org.openstreetmap.josm.data.osm.Node; 024import org.openstreetmap.josm.data.osm.NodePositionComparator; 025import org.openstreetmap.josm.data.osm.Way; 026 027/** 028 * Some tools for geometry related tasks. 029 * 030 * @author viesturs 031 */ 032public final class Geometry { 033 034 private Geometry() { 035 // Hide default constructor for utils classes 036 } 037 038 public enum PolygonIntersection {FIRST_INSIDE_SECOND, SECOND_INSIDE_FIRST, OUTSIDE, CROSSING} 039 040 /** 041 * Will find all intersection and add nodes there for list of given ways. 042 * Handles self-intersections too. 043 * And makes commands to add the intersection points to ways. 044 * 045 * Prerequisite: no two nodes have the same coordinates. 046 * 047 * @param ways a list of ways to test 048 * @param test if false, do not build list of Commands, just return nodes 049 * @param cmds list of commands, typically empty when handed to this method. 050 * Will be filled with commands that add intersection nodes to 051 * the ways. 052 * @return list of new nodes 053 */ 054 public static Set<Node> addIntersections(List<Way> ways, boolean test, List<Command> cmds) { 055 056 int n = ways.size(); 057 @SuppressWarnings("unchecked") 058 List<Node>[] newNodes = new ArrayList[n]; 059 BBox[] wayBounds = new BBox[n]; 060 boolean[] changedWays = new boolean[n]; 061 062 Set<Node> intersectionNodes = new LinkedHashSet<Node>(); 063 064 //copy node arrays for local usage. 065 for (int pos = 0; pos < n; pos ++) { 066 newNodes[pos] = new ArrayList<Node>(ways.get(pos).getNodes()); 067 wayBounds[pos] = getNodesBounds(newNodes[pos]); 068 changedWays[pos] = false; 069 } 070 071 //iterate over all way pairs and introduce the intersections 072 Comparator<Node> coordsComparator = new NodePositionComparator(); 073 for (int seg1Way = 0; seg1Way < n; seg1Way ++) { 074 for (int seg2Way = seg1Way; seg2Way < n; seg2Way ++) { 075 076 //do not waste time on bounds that do not intersect 077 if (!wayBounds[seg1Way].intersects(wayBounds[seg2Way])) { 078 continue; 079 } 080 081 List<Node> way1Nodes = newNodes[seg1Way]; 082 List<Node> way2Nodes = newNodes[seg2Way]; 083 084 //iterate over primary segmemt 085 for (int seg1Pos = 0; seg1Pos + 1 < way1Nodes.size(); seg1Pos ++) { 086 087 //iterate over secondary segment 088 int seg2Start = seg1Way != seg2Way ? 0: seg1Pos + 2;//skip the adjacent segment 089 090 for (int seg2Pos = seg2Start; seg2Pos + 1< way2Nodes.size(); seg2Pos ++) { 091 092 //need to get them again every time, because other segments may be changed 093 Node seg1Node1 = way1Nodes.get(seg1Pos); 094 Node seg1Node2 = way1Nodes.get(seg1Pos + 1); 095 Node seg2Node1 = way2Nodes.get(seg2Pos); 096 Node seg2Node2 = way2Nodes.get(seg2Pos + 1); 097 098 int commonCount = 0; 099 //test if we have common nodes to add. 100 if (seg1Node1 == seg2Node1 || seg1Node1 == seg2Node2) { 101 commonCount ++; 102 103 if (seg1Way == seg2Way && 104 seg1Pos == 0 && 105 seg2Pos == way2Nodes.size() -2) { 106 //do not add - this is first and last segment of the same way. 107 } else { 108 intersectionNodes.add(seg1Node1); 109 } 110 } 111 112 if (seg1Node2 == seg2Node1 || seg1Node2 == seg2Node2) { 113 commonCount ++; 114 115 intersectionNodes.add(seg1Node2); 116 } 117 118 //no common nodes - find intersection 119 if (commonCount == 0) { 120 EastNorth intersection = getSegmentSegmentIntersection( 121 seg1Node1.getEastNorth(), seg1Node2.getEastNorth(), 122 seg2Node1.getEastNorth(), seg2Node2.getEastNorth()); 123 124 if (intersection != null) { 125 if (test) { 126 intersectionNodes.add(seg2Node1); 127 return intersectionNodes; 128 } 129 130 Node newNode = new Node(Main.getProjection().eastNorth2latlon(intersection)); 131 Node intNode = newNode; 132 boolean insertInSeg1 = false; 133 boolean insertInSeg2 = false; 134 //find if the intersection point is at end point of one of the segments, if so use that point 135 136 //segment 1 137 if (coordsComparator.compare(newNode, seg1Node1) == 0) { 138 intNode = seg1Node1; 139 } else if (coordsComparator.compare(newNode, seg1Node2) == 0) { 140 intNode = seg1Node2; 141 } else { 142 insertInSeg1 = true; 143 } 144 145 //segment 2 146 if (coordsComparator.compare(newNode, seg2Node1) == 0) { 147 intNode = seg2Node1; 148 } else if (coordsComparator.compare(newNode, seg2Node2) == 0) { 149 intNode = seg2Node2; 150 } else { 151 insertInSeg2 = true; 152 } 153 154 if (insertInSeg1) { 155 way1Nodes.add(seg1Pos +1, intNode); 156 changedWays[seg1Way] = true; 157 158 //fix seg2 position, as indexes have changed, seg2Pos is always bigger than seg1Pos on the same segment. 159 if (seg2Way == seg1Way) { 160 seg2Pos ++; 161 } 162 } 163 164 if (insertInSeg2) { 165 way2Nodes.add(seg2Pos +1, intNode); 166 changedWays[seg2Way] = true; 167 168 //Do not need to compare again to already split segment 169 seg2Pos ++; 170 } 171 172 intersectionNodes.add(intNode); 173 174 if (intNode == newNode) { 175 cmds.add(new AddCommand(intNode)); 176 } 177 } 178 } 179 else if (test && !intersectionNodes.isEmpty()) 180 return intersectionNodes; 181 } 182 } 183 } 184 } 185 186 187 for (int pos = 0; pos < ways.size(); pos ++) { 188 if (!changedWays[pos]) { 189 continue; 190 } 191 192 Way way = ways.get(pos); 193 Way newWay = new Way(way); 194 newWay.setNodes(newNodes[pos]); 195 196 cmds.add(new ChangeCommand(way, newWay)); 197 } 198 199 return intersectionNodes; 200 } 201 202 private static BBox getNodesBounds(List<Node> nodes) { 203 204 BBox bounds = new BBox(nodes.get(0)); 205 for (Node n: nodes) { 206 bounds.add(n.getCoor()); 207 } 208 return bounds; 209 } 210 211 /** 212 * Tests if given point is to the right side of path consisting of 3 points. 213 * 214 * (Imagine the path is continued beyond the endpoints, so you get two rays 215 * starting from lineP2 and going through lineP1 and lineP3 respectively 216 * which divide the plane into two parts. The test returns true, if testPoint 217 * lies in the part that is to the right when traveling in the direction 218 * lineP1, lineP2, lineP3.) 219 * 220 * @param lineP1 first point in path 221 * @param lineP2 second point in path 222 * @param lineP3 third point in path 223 * @param testPoint 224 * @return true if to the right side, false otherwise 225 */ 226 public static boolean isToTheRightSideOfLine(Node lineP1, Node lineP2, Node lineP3, Node testPoint) { 227 boolean pathBendToRight = angleIsClockwise(lineP1, lineP2, lineP3); 228 boolean rightOfSeg1 = angleIsClockwise(lineP1, lineP2, testPoint); 229 boolean rightOfSeg2 = angleIsClockwise(lineP2, lineP3, testPoint); 230 231 if (pathBendToRight) 232 return rightOfSeg1 && rightOfSeg2; 233 else 234 return !(!rightOfSeg1 && !rightOfSeg2); 235 } 236 237 /** 238 * This method tests if secondNode is clockwise to first node. 239 * @param commonNode starting point for both vectors 240 * @param firstNode first vector end node 241 * @param secondNode second vector end node 242 * @return true if first vector is clockwise before second vector. 243 */ 244 public static boolean angleIsClockwise(Node commonNode, Node firstNode, Node secondNode) { 245 return angleIsClockwise(commonNode.getEastNorth(), firstNode.getEastNorth(), secondNode.getEastNorth()); 246 } 247 248 /** 249 * Finds the intersection of two line segments 250 * @return EastNorth null if no intersection was found, the EastNorth coordinates of the intersection otherwise 251 */ 252 public static EastNorth getSegmentSegmentIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 253 254 CheckParameterUtil.ensureValidCoordinates(p1, "p1"); 255 CheckParameterUtil.ensureValidCoordinates(p2, "p2"); 256 CheckParameterUtil.ensureValidCoordinates(p3, "p3"); 257 CheckParameterUtil.ensureValidCoordinates(p4, "p4"); 258 259 double x1 = p1.getX(); 260 double y1 = p1.getY(); 261 double x2 = p2.getX(); 262 double y2 = p2.getY(); 263 double x3 = p3.getX(); 264 double y3 = p3.getY(); 265 double x4 = p4.getX(); 266 double y4 = p4.getY(); 267 268 //TODO: do this locally. 269 //TODO: remove this check after careful testing 270 if (!Line2D.linesIntersect(x1, y1, x2, y2, x3, y3, x4, y4)) return null; 271 272 // solve line-line intersection in parametric form: 273 // (x1,y1) + (x2-x1,y2-y1)* u = (x3,y3) + (x4-x3,y4-y3)* v 274 // (x2-x1,y2-y1)*u - (x4-x3,y4-y3)*v = (x3-x1,y3-y1) 275 // if 0<= u,v <=1, intersection exists at ( x1+ (x2-x1)*u, y1 + (y2-y1)*u ) 276 277 double a1 = x2 - x1; 278 double b1 = x3 - x4; 279 double c1 = x3 - x1; 280 281 double a2 = y2 - y1; 282 double b2 = y3 - y4; 283 double c2 = y3 - y1; 284 285 // Solve the equations 286 double det = a1*b2 - a2*b1; 287 288 double uu = b2*c1 - b1*c2 ; 289 double vv = a1*c2 - a2*c1; 290 double mag = Math.abs(uu)+Math.abs(vv); 291 292 if (Math.abs(det) > 1e-12 * mag) { 293 double u = uu/det, v = vv/det; 294 if (u>-1e-8 && u < 1+1e-8 && v>-1e-8 && v < 1+1e-8 ) { 295 if (u<0) u=0; 296 if (u>1) u=1.0; 297 return new EastNorth(x1+a1*u, y1+a2*u); 298 } else { 299 return null; 300 } 301 } else { 302 // parallel lines 303 return null; 304 } 305 } 306 307 /** 308 * Finds the intersection of two lines of infinite length. 309 * @return EastNorth null if no intersection was found, the coordinates of the intersection otherwise 310 * @throws IllegalArgumentException if a parameter is null or without valid coordinates 311 */ 312 public static EastNorth getLineLineIntersection(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 313 314 CheckParameterUtil.ensureValidCoordinates(p1, "p1"); 315 CheckParameterUtil.ensureValidCoordinates(p2, "p2"); 316 CheckParameterUtil.ensureValidCoordinates(p3, "p3"); 317 CheckParameterUtil.ensureValidCoordinates(p4, "p4"); 318 319 if (!p1.isValid()) throw new IllegalArgumentException(); 320 321 // Convert line from (point, point) form to ax+by=c 322 double a1 = p2.getY() - p1.getY(); 323 double b1 = p1.getX() - p2.getX(); 324 double c1 = p2.getX() * p1.getY() - p1.getX() * p2.getY(); 325 326 double a2 = p4.getY() - p3.getY(); 327 double b2 = p3.getX() - p4.getX(); 328 double c2 = p4.getX() * p3.getY() - p3.getX() * p4.getY(); 329 330 // Solve the equations 331 double det = a1 * b2 - a2 * b1; 332 if (det == 0) 333 return null; // Lines are parallel 334 335 return new EastNorth((b1 * c2 - b2 * c1) / det, (a2 * c1 - a1 * c2) / det); 336 } 337 338 public static boolean segmentsParallel(EastNorth p1, EastNorth p2, EastNorth p3, EastNorth p4) { 339 340 CheckParameterUtil.ensureValidCoordinates(p1, "p1"); 341 CheckParameterUtil.ensureValidCoordinates(p2, "p2"); 342 CheckParameterUtil.ensureValidCoordinates(p3, "p3"); 343 CheckParameterUtil.ensureValidCoordinates(p4, "p4"); 344 345 // Convert line from (point, point) form to ax+by=c 346 double a1 = p2.getY() - p1.getY(); 347 double b1 = p1.getX() - p2.getX(); 348 349 double a2 = p4.getY() - p3.getY(); 350 double b2 = p3.getX() - p4.getX(); 351 352 // Solve the equations 353 double det = a1 * b2 - a2 * b1; 354 // remove influence of of scaling factor 355 det /= Math.sqrt(a1*a1 + b1*b1) * Math.sqrt(a2*a2 + b2*b2); 356 return Math.abs(det) < 1e-3; 357 } 358 359 private static EastNorth closestPointTo(EastNorth p1, EastNorth p2, EastNorth point, boolean segmentOnly) { 360 CheckParameterUtil.ensureParameterNotNull(p1, "p1"); 361 CheckParameterUtil.ensureParameterNotNull(p2, "p2"); 362 CheckParameterUtil.ensureParameterNotNull(point, "point"); 363 364 double ldx = p2.getX() - p1.getX(); 365 double ldy = p2.getY() - p1.getY(); 366 367 if (ldx == 0 && ldy == 0) //segment zero length 368 return p1; 369 370 double pdx = point.getX() - p1.getX(); 371 double pdy = point.getY() - p1.getY(); 372 373 double offset = (pdx * ldx + pdy * ldy) / (ldx * ldx + ldy * ldy); 374 375 if (segmentOnly && offset <= 0) 376 return p1; 377 else if (segmentOnly && offset >= 1) 378 return p2; 379 else 380 return new EastNorth(p1.getX() + ldx * offset, p1.getY() + ldy * offset); 381 } 382 383 /** 384 * Calculates closest point to a line segment. 385 * @param segmentP1 First point determining line segment 386 * @param segmentP2 Second point determining line segment 387 * @param point Point for which a closest point is searched on line segment [P1,P2] 388 * @return segmentP1 if it is the closest point, segmentP2 if it is the closest point, 389 * a new point if closest point is between segmentP1 and segmentP2. 390 * @since 3650 391 * @see #closestPointToLine 392 */ 393 public static EastNorth closestPointToSegment(EastNorth segmentP1, EastNorth segmentP2, EastNorth point) { 394 return closestPointTo(segmentP1, segmentP2, point, true); 395 } 396 397 /** 398 * Calculates closest point to a line. 399 * @param lineP1 First point determining line 400 * @param lineP2 Second point determining line 401 * @param point Point for which a closest point is searched on line (P1,P2) 402 * @return The closest point found on line. It may be outside the segment [P1,P2]. 403 * @since 4134 404 * @see #closestPointToSegment 405 */ 406 public static EastNorth closestPointToLine(EastNorth lineP1, EastNorth lineP2, EastNorth point) { 407 return closestPointTo(lineP1, lineP2, point, false); 408 } 409 410 /** 411 * This method tests if secondNode is clockwise to first node. 412 * 413 * The line through the two points commonNode and firstNode divides the 414 * plane into two parts. The test returns true, if secondNode lies in 415 * the part that is to the right when traveling in the direction from 416 * commonNode to firstNode. 417 * 418 * @param commonNode starting point for both vectors 419 * @param firstNode first vector end node 420 * @param secondNode second vector end node 421 * @return true if first vector is clockwise before second vector. 422 */ 423 public static boolean angleIsClockwise(EastNorth commonNode, EastNorth firstNode, EastNorth secondNode) { 424 425 CheckParameterUtil.ensureValidCoordinates(commonNode, "commonNode"); 426 CheckParameterUtil.ensureValidCoordinates(firstNode, "firstNode"); 427 CheckParameterUtil.ensureValidCoordinates(secondNode, "secondNode"); 428 429 double dy1 = (firstNode.getY() - commonNode.getY()); 430 double dy2 = (secondNode.getY() - commonNode.getY()); 431 double dx1 = (firstNode.getX() - commonNode.getX()); 432 double dx2 = (secondNode.getX() - commonNode.getX()); 433 434 return dy1 * dx2 - dx1 * dy2 > 0; 435 } 436 437 private static Area getArea(List<Node> polygon) { 438 Path2D path = new Path2D.Double(); 439 440 boolean begin = true; 441 for (Node n : polygon) { 442 if (begin) { 443 path.moveTo(n.getEastNorth().getX(), n.getEastNorth().getY()); 444 begin = false; 445 } else { 446 path.lineTo(n.getEastNorth().getX(), n.getEastNorth().getY()); 447 } 448 } 449 if (!begin) { 450 path.closePath(); 451 } 452 453 return new Area(path); 454 } 455 456 /** 457 * Tests if two polygons intersect. 458 * @param first 459 * @param second 460 * @return intersection kind 461 */ 462 public static PolygonIntersection polygonIntersection(List<Node> first, List<Node> second) { 463 464 Area a1 = getArea(first); 465 Area a2 = getArea(second); 466 467 Area inter = new Area(a1); 468 inter.intersect(a2); 469 470 Rectangle bounds = inter.getBounds(); 471 472 if (inter.isEmpty() || bounds.getHeight()*bounds.getWidth() <= 1.0) { 473 return PolygonIntersection.OUTSIDE; 474 } else if (inter.equals(a1)) { 475 return PolygonIntersection.FIRST_INSIDE_SECOND; 476 } else if (inter.equals(a2)) { 477 return PolygonIntersection.SECOND_INSIDE_FIRST; 478 } else { 479 return PolygonIntersection.CROSSING; 480 } 481 } 482 483 /** 484 * Tests if point is inside a polygon. The polygon can be self-intersecting. In such case the contains function works in xor-like manner. 485 * @param polygonNodes list of nodes from polygon path. 486 * @param point the point to test 487 * @return true if the point is inside polygon. 488 */ 489 public static boolean nodeInsidePolygon(Node point, List<Node> polygonNodes) { 490 if (polygonNodes.size() < 2) 491 return false; 492 493 boolean inside = false; 494 Node p1, p2; 495 496 //iterate each side of the polygon, start with the last segment 497 Node oldPoint = polygonNodes.get(polygonNodes.size() - 1); 498 499 for (Node newPoint : polygonNodes) { 500 //skip duplicate points 501 if (newPoint.equals(oldPoint)) { 502 continue; 503 } 504 505 //order points so p1.lat <= p2.lat 506 if (newPoint.getEastNorth().getY() > oldPoint.getEastNorth().getY()) { 507 p1 = oldPoint; 508 p2 = newPoint; 509 } else { 510 p1 = newPoint; 511 p2 = oldPoint; 512 } 513 514 //test if the line is crossed and if so invert the inside flag. 515 if ((newPoint.getEastNorth().getY() < point.getEastNorth().getY()) == (point.getEastNorth().getY() <= oldPoint.getEastNorth().getY()) 516 && (point.getEastNorth().getX() - p1.getEastNorth().getX()) * (p2.getEastNorth().getY() - p1.getEastNorth().getY()) 517 < (p2.getEastNorth().getX() - p1.getEastNorth().getX()) * (point.getEastNorth().getY() - p1.getEastNorth().getY())) 518 { 519 inside = !inside; 520 } 521 522 oldPoint = newPoint; 523 } 524 525 return inside; 526 } 527 528 /** 529 * Returns area of a closed way in square meters. 530 * (approximate(?), but should be OK for small areas) 531 * 532 * Relies on the current projection: Works correctly, when 533 * one unit in projected coordinates corresponds to one meter. 534 * This is true for most projections, but not for WGS84 and 535 * Mercator (EPSG:3857). 536 * 537 * @param way Way to measure, should be closed (first node is the same as last node) 538 * @return area of the closed way. 539 */ 540 public static double closedWayArea(Way way) { 541 542 //http://local.wasp.uwa.edu.au/~pbourke/geometry/polyarea/ 543 double area = 0; 544 Node lastN = null; 545 for (Node n : way.getNodes()) { 546 if (lastN != null) { 547 n.getEastNorth().getX(); 548 549 area += (calcX(n) * calcY(lastN)) - (calcY(n) * calcX(lastN)); 550 } 551 lastN = n; 552 } 553 return Math.abs(area/2); 554 } 555 556 protected static double calcX(Node p1){ 557 double lat1, lon1, lat2, lon2; 558 double dlon, dlat; 559 560 lat1 = p1.getCoor().lat() * Math.PI / 180.0; 561 lon1 = p1.getCoor().lon() * Math.PI / 180.0; 562 lat2 = lat1; 563 lon2 = 0; 564 565 dlon = lon2 - lon1; 566 dlat = lat2 - lat1; 567 568 double a = (Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2)); 569 double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 570 return 6367000 * c; 571 } 572 573 protected static double calcY(Node p1){ 574 double lat1, lon1, lat2, lon2; 575 double dlon, dlat; 576 577 lat1 = p1.getCoor().lat() * Math.PI / 180.0; 578 lon1 = p1.getCoor().lon() * Math.PI / 180.0; 579 lat2 = 0; 580 lon2 = lon1; 581 582 dlon = lon2 - lon1; 583 dlat = lat2 - lat1; 584 585 double a = (Math.pow(Math.sin(dlat/2), 2) + Math.cos(lat1) * Math.cos(lat2) * Math.pow(Math.sin(dlon/2), 2)); 586 double c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a)); 587 return 6367000 * c; 588 } 589 590 /** 591 * Determines whether a way is oriented clockwise. 592 * 593 * Internals: Assuming a closed non-looping way, compute twice the area 594 * of the polygon using the formula {@code 2 * area = sum (X[n] * Y[n+1] - X[n+1] * Y[n])}. 595 * If the area is negative the way is ordered in a clockwise direction. 596 * 597 * See http://paulbourke.net/geometry/polyarea/ 598 * 599 * @param w the way to be checked. 600 * @return true if and only if way is oriented clockwise. 601 * @throws IllegalArgumentException if way is not closed (see {@link Way#isClosed}). 602 */ 603 public static boolean isClockwise(Way w) { 604 if (!w.isClosed()) { 605 throw new IllegalArgumentException("Way must be closed to check orientation."); 606 } 607 608 double area2 = 0.; 609 int nodesCount = w.getNodesCount(); 610 611 for (int node = 1; node <= /*sic! consider last-first as well*/ nodesCount; node++) { 612 LatLon coorPrev = w.getNode(node - 1).getCoor(); 613 LatLon coorCurr = w.getNode(node % nodesCount).getCoor(); 614 area2 += coorPrev.lon() * coorCurr.lat(); 615 area2 -= coorCurr.lon() * coorPrev.lat(); 616 } 617 return area2 < 0; 618 } 619 620 /** 621 * Returns angle of a segment defined with 2 point coordinates. 622 * 623 * @param p1 624 * @param p2 625 * @return Angle in radians (-pi, pi] 626 */ 627 public static double getSegmentAngle(EastNorth p1, EastNorth p2) { 628 629 CheckParameterUtil.ensureValidCoordinates(p1, "p1"); 630 CheckParameterUtil.ensureValidCoordinates(p2, "p2"); 631 632 return Math.atan2(p2.north() - p1.north(), p2.east() - p1.east()); 633 } 634 635 /** 636 * Returns angle of a corner defined with 3 point coordinates. 637 * 638 * @param p1 639 * @param p2 Common endpoint 640 * @param p3 641 * @return Angle in radians (-pi, pi] 642 */ 643 public static double getCornerAngle(EastNorth p1, EastNorth p2, EastNorth p3) { 644 645 CheckParameterUtil.ensureValidCoordinates(p1, "p1"); 646 CheckParameterUtil.ensureValidCoordinates(p2, "p2"); 647 CheckParameterUtil.ensureValidCoordinates(p3, "p3"); 648 649 Double result = getSegmentAngle(p2, p1) - getSegmentAngle(p2, p3); 650 if (result <= -Math.PI) { 651 result += 2 * Math.PI; 652 } 653 654 if (result > Math.PI) { 655 result -= 2 * Math.PI; 656 } 657 658 return result; 659 } 660 661 /** 662 * Compute the centroid of nodes 663 * @param nodes Nodes for which the centroid is wanted 664 * @return the centroid of nodes 665 */ 666 public static EastNorth getCentroid(List<Node> nodes) { 667 668 BigDecimal area = BigDecimal.ZERO; 669 BigDecimal north = BigDecimal.ZERO; 670 BigDecimal east = BigDecimal.ZERO; 671 672 // See http://en.wikipedia.org/w/index.php?title=Centroid&oldid=294224857#Centroid_of_polygon for the equation used here 673 for (int i = 0; i < nodes.size(); i++) { 674 EastNorth n0 = nodes.get(i).getEastNorth(); 675 EastNorth n1 = nodes.get((i+1) % nodes.size()).getEastNorth(); 676 677 if (n0.isValid() && n1.isValid()) { 678 BigDecimal x0 = new BigDecimal(n0.east()); 679 BigDecimal y0 = new BigDecimal(n0.north()); 680 BigDecimal x1 = new BigDecimal(n1.east()); 681 BigDecimal y1 = new BigDecimal(n1.north()); 682 683 BigDecimal k = x0.multiply(y1, MathContext.DECIMAL128).subtract(y0.multiply(x1, MathContext.DECIMAL128)); 684 685 area = area.add(k, MathContext.DECIMAL128); 686 east = east.add(k.multiply(x0.add(x1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 687 north = north.add(k.multiply(y0.add(y1, MathContext.DECIMAL128), MathContext.DECIMAL128)); 688 } 689 } 690 691 BigDecimal d = new BigDecimal(3, MathContext.DECIMAL128); // 1/2 * 6 = 3 692 area = area.multiply(d, MathContext.DECIMAL128); 693 if (area.compareTo(BigDecimal.ZERO) != 0) { 694 north = north.divide(area, MathContext.DECIMAL128); 695 east = east.divide(area, MathContext.DECIMAL128); 696 } 697 698 return new EastNorth(east.doubleValue(), north.doubleValue()); 699 } 700 701 /** 702 * Returns the coordinate of intersection of segment sp1-sp2 and an altitude 703 * to it starting at point ap. If the line defined with sp1-sp2 intersects 704 * its altitude out of sp1-sp2, null is returned. 705 * 706 * @param sp1 707 * @param sp2 708 * @param ap 709 * @return Intersection coordinate or null 710 */ 711 public static EastNorth getSegmentAltituteIntersection(EastNorth sp1, EastNorth sp2, EastNorth ap) { 712 713 CheckParameterUtil.ensureValidCoordinates(sp1, "sp1"); 714 CheckParameterUtil.ensureValidCoordinates(sp2, "sp2"); 715 CheckParameterUtil.ensureValidCoordinates(ap, "ap"); 716 717 Double segmentLenght = sp1.distance(sp2); 718 Double altitudeAngle = getSegmentAngle(sp1, sp2) + Math.PI / 2; 719 720 // Taking a random point on the altitude line (angle is known). 721 EastNorth ap2 = new EastNorth(ap.east() + 1000 722 * Math.cos(altitudeAngle), ap.north() + 1000 723 * Math.sin(altitudeAngle)); 724 725 // Finding the intersection of two lines 726 EastNorth resultCandidate = Geometry.getLineLineIntersection(sp1, sp2, 727 ap, ap2); 728 729 // Filtering result 730 if (resultCandidate != null 731 && resultCandidate.distance(sp1) * .999 < segmentLenght 732 && resultCandidate.distance(sp2) * .999 < segmentLenght) { 733 return resultCandidate; 734 } else { 735 return null; 736 } 737 } 738}