Changeset 6934 in josm for trunk/src/org/openstreetmap/josm/tools
 Timestamp:
 20140326T02:15:38+01:00 (10 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/org/openstreetmap/josm/tools/Geometry.java
r6920 r6934 688 688 * @param nodes Nodes for which the centroid is wanted 689 689 * @return the centroid of nodes 690 * @see Geometry#getCenter 690 691 */ 691 692 public static EastNorth getCentroid(List<Node> nodes) { … … 724 725 } 725 726 727 /** 728 * Compute center of the circle closest to different nodes. 729 * 730 * Ensure exact center computation in case nodes are already aligned in circle. 731 * This is done by least square method. 732 * Let be a_i x + b_i y + c_i = 0 equations of bisectors of each edges. 733 * Center must be intersection of all bisectors. 734 * <pre> 735 * [ a1 b1 ] [ c1 ] 736 * With A = [ ... ... ] and Y = [ ... ] 737 * [ an bn ] [ cn ] 738 * </pre> 739 * An approximation of center of circle is (At.A)^1.At.Y 740 * @param nodes Nodes parts of the circle (at least 3) 741 * @return An approximation of the center, of null if there is no solution. 742 * @see Geometry#getCentroid 743 * @since 6934 744 */ 745 public static EastNorth getCenter(List<Node> nodes) { 746 int nc = nodes.size(); 747 if(nc < 3) return null; 748 /** 749 * Equation of each bisector ax + by + c = 0 750 */ 751 double[] a = new double[nc]; 752 double[] b = new double[nc]; 753 double[] c = new double[nc]; 754 // Compute equation of bisector 755 for(int i = 0; i < nc; i++) { 756 EastNorth pt1 = nodes.get(i).getEastNorth(); 757 EastNorth pt2 = nodes.get((i+1) % nc).getEastNorth(); 758 a[i] = pt1.east()  pt2.east(); 759 b[i] = pt1.north()  pt2.north(); 760 double d = Math.sqrt(a[i]*a[i] + b[i]*b[i]); 761 if(d == 0) return null; 762 a[i] /= d; 763 b[i] /= d; 764 double xC = (pt1.east() + pt2.east()) / 2; 765 double yC = (pt1.north() + pt2.north()) / 2; 766 c[i] = (a[i]*xC + b[i]*yC); 767 } 768 // At.A = [aij] 769 double a11 = 0, a12 = 0, a22 = 0; 770 // At.Y = [bi] 771 double b1 = 0, b2 = 0; 772 for(int i = 0; i < nc; i++) { 773 a11 += a[i]*a[i]; 774 a12 += a[i]*b[i]; 775 a22 += b[i]*b[i]; 776 b1 = a[i]*c[i]; 777 b2 = b[i]*c[i]; 778 } 779 // (At.A)^1 = [invij] 780 double det = a11*a22  a12*a12; 781 if(Math.abs(det) < 1e5) return null; 782 double inv11 = a22/det; 783 double inv12 = a12/det; 784 double inv22 = a11/det; 785 // center (xC, yC) = (At.A)^1.At.y 786 double xC = inv11*b1 + inv12*b2; 787 double yC = inv12*b1 + inv22*b2; 788 return new EastNorth(xC, yC); 789 } 790 726 791 /** 727 792 * Returns the coordinate of intersection of segment sp1sp2 and an altitude
Note:
See TracChangeset
for help on using the changeset viewer.