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