| 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 | */ |
| 744 | public static EastNorth getCenter(List<Node> nodes) { |
| 745 | int nc = nodes.size(); |
| 746 | if(nc < 3) return null; |
| 747 | /** |
| 748 | * Equation of each bisector ax + by + c = 0 |
| 749 | */ |
| 750 | double[] a = new double[nc]; |
| 751 | double[] b = new double[nc]; |
| 752 | double[] c = new double[nc]; |
| 753 | // Compute equation of bisector |
| 754 | for(int i = 0; i < nc; i++) { |
| 755 | EastNorth pt1 = nodes.get(i).getEastNorth(); |
| 756 | EastNorth pt2 = nodes.get((i+1) % nc).getEastNorth(); |
| 757 | a[i] = pt1.east() - pt2.east(); |
| 758 | b[i] = pt1.north() - pt2.north(); |
| 759 | double d = Math.sqrt(a[i]*a[i] + b[i]*b[i]); |
| 760 | if(d == 0) return null; |
| 761 | a[i] /= d; |
| 762 | b[i] /= d; |
| 763 | double xC = (pt1.east() + pt2.east()) / 2; |
| 764 | double yC = (pt1.north() + pt2.north()) / 2; |
| 765 | c[i] = -(a[i]*xC + b[i]*yC); |
| 766 | } |
| 767 | // At.A = [aij] |
| 768 | double a11 = 0, a12 = 0, a22 = 0; |
| 769 | // At.Y = [bi] |
| 770 | double b1 = 0, b2 = 0; |
| 771 | for(int i = 0; i < nc; i++) { |
| 772 | a11 += a[i]*a[i]; |
| 773 | a12 += a[i]*b[i]; |
| 774 | a22 += b[i]*b[i]; |
| 775 | b1 -= a[i]*c[i]; |
| 776 | b2 -= b[i]*c[i]; |
| 777 | } |
| 778 | // (At.A)^-1 = [invij] |
| 779 | double det = a11*a22 - a12*a12; |
| 780 | if(Math.abs(det) < 1e-5) return null; |
| 781 | double inv11 = a22/det; |
| 782 | double inv12 = -a12/det; |
| 783 | double inv22 = a11/det; |
| 784 | // center (xC, yC) = (At.A)^-1.At.y |
| 785 | double xC = inv11*b1 + inv12*b2; |
| 786 | double yC = inv12*b1 + inv22*b2; |
| 787 | return new EastNorth(xC, yC); |
| 788 | } |
| 789 | |
| 790 | /** |