source: josm/trunk/src/org/openstreetmap/josm/tools/Geometry.java@ 6566

Last change on this file since 6566 was 6566, checked in by simon04, 10 years ago

fix #9494 - Advanced object info: add "Center of bounding box", and for ways "Centroid"
  • Property svn:eol-style set to native
File size: 27.7 KB
Line 
1// License: GPL. For details, see LICENSE file.
2package org.openstreetmap.josm.tools;
3
4import java.awt.Rectangle;
5import java.awt.geom.Area;
6import java.awt.geom.Line2D;
7import java.awt.geom.Path2D;
8import java.math.BigDecimal;
9import java.math.MathContext;
10import java.util.ArrayList;
11import java.util.Comparator;
12import java.util.LinkedHashSet;
13import java.util.List;
14import java.util.Set;
15
16import org.openstreetmap.josm.Main;
17import org.openstreetmap.josm.command.AddCommand;
18import org.openstreetmap.josm.command.ChangeCommand;
19import org.openstreetmap.josm.command.Command;
20import org.openstreetmap.josm.data.coor.EastNorth;
21import org.openstreetmap.josm.data.coor.LatLon;
22import org.openstreetmap.josm.data.osm.BBox;
23import org.openstreetmap.josm.data.osm.Node;
24import org.openstreetmap.josm.data.osm.NodePositionComparator;
25import org.openstreetmap.josm.data.osm.Way;
26
27/**
28 * Some tools for geometry related tasks.
29 *
30 * @author viesturs
31 */
32public final class Geometry {
33
34 private Geometry() {
35 // Hide default constructor for utils classes
36 }
37
38 public enum PolygonIntersection {FIRST_INSIDE_SECOND, SECOND_INSIDE_FIRST, OUTSIDE, CROSSING}
39
40 /**
41 * Will find all intersection and add nodes there for list of given ways.
42 * Handles self-intersections too.
43 * And makes commands to add the intersection points to ways.
44 *
45 * Prerequisite: no two nodes have the same coordinates.
46 *
47 * @param ways a list of ways to test
48 * @param test if false, do not build list of Commands, just return nodes
49 * @param cmds list of commands, typically empty when handed to this method.
50 * Will be filled with commands that add intersection nodes to
51 * the ways.
52 * @return list of new nodes
53 */
54 public static Set<Node> addIntersections(List<Way> ways, boolean test, List<Command> cmds) {
55
56 int n = ways.size();
57 @SuppressWarnings("unchecked")
58 List<Node>[] newNodes = new ArrayList[n];
59 BBox[] wayBounds = new BBox[n];
60 boolean[] changedWays = new boolean[n];
61
62 Set<Node> intersectionNodes = new LinkedHashSet<Node>();
63
64 //copy node arrays for local usage.
65 for (int pos = 0; pos < n; pos ++) {
66 newNodes[pos] = new ArrayList<Node>(ways.get(pos).getNodes());
67 wayBounds[pos] = getNodesBounds(newNodes[pos]);
68 changedWays[pos] = false;
69 }
70
71 //iterate over all way pairs and introduce the intersections
72 Comparator<Node> coordsComparator = new NodePositionComparator();
73 for (int seg1Way = 0; seg1Way < n; seg1Way ++) {
74 for (int seg2Way = seg1Way; seg2Way < n; seg2Way ++) {
75
76 //do not waste time on bounds that do not intersect
77 if (!wayBounds[seg1Way].intersects(wayBounds[seg2Way])) {
78 continue;
79 }
80
81 List<Node> way1Nodes = newNodes[seg1Way];
82 List<Node> way2Nodes = newNodes[seg2Way];
83
84 //iterate over primary segmemt
85 for (int seg1Pos = 0; seg1Pos + 1 < way1Nodes.size(); seg1Pos ++) {
86
87 //iterate over secondary segment
88 int seg2Start = seg1Way != seg2Way ? 0: seg1Pos + 2;//skip the adjacent segment
89
90 for (int seg2Pos = seg2Start; seg2Pos + 1< way2Nodes.size(); seg2Pos ++) {
91
92 //need to get them again every time, because other segments may be changed
93 Node seg1Node1 = way1Nodes.get(seg1Pos);
94 Node seg1Node2 = way1Nodes.get(seg1Pos + 1);
95 Node seg2Node1 = way2Nodes.get(seg2Pos);
96 Node seg2Node2 = way2Nodes.get(seg2Pos + 1);
97
98 int commonCount = 0;
99 //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/barycenter 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}
Note: See TracBrowser for help on using the repository browser.