#include<cstdio>#include<cstdlib>#include<cmath>#define eps 1e-8#define max 1000000 struct Point { float x ; float y ;};Point point[max];using namespace std ; int main(){ int n = 0 ; state: printf("请输入点的个数(上限不要超过1000,000):\n"); scanf("%d",&n) ; //input the number of points while(n<3) { printf("输入点个数小于三个,无法确定唯一圆,请重新输入!\n"); scanf("%d",&n); } double sum_x = 0 ; double sum_y = 0 ; double sum_xy = 0 ; double sum_x_2 = 0 ; double sum_y_2 = 0 ; double sum_x_y_2 = 0 ; double sum_y_x_2 = 0 ; double sum_x_3 = 0 ; double sum_y_3 = 0 ; double temp_x_2 = 0 ; double temp_y_2 = 0 ; double P , Q ,R ,T , W ; printf("请输入点的坐标,中间以空格间隔:\n"); for(int i=0;i<n;++i) //input n point location { scanf("%f %f",&point[i].x,&point[i].y); sum_x += point[i].x; sum_y += point[i].y; sum_xy +=point[i].x*point[i].y ; temp_x_2 = point[i].x*point[i].x ; temp_y_2 = point[i].y*point[i].y ; sum_x_2 += temp_x_2 ; sum_y_2 += temp_y_2 ; sum_x_y_2 += temp_y_2*point[i].x ; sum_y_x_2 += temp_x_2*point[i].y ; sum_x_3 += temp_x_2*point[i].x ; sum_y_3 += temp_y_2*point[i].y ; } P = n*sum_x_2 - sum_x*sum_x ; Q = n*sum_xy - sum_x*sum_y ; R = n*(sum_x_3+sum_x_y_2)-(sum_x_2+sum_y_2)*sum_x ; T = n*sum_y_2 - sum_y*sum_y ; W = n*(sum_y_x_2+sum_y_3)-(sum_x_2+sum_y_2)*sum_y ; if(abs(Q*Q-P*T)<eps) { printf("ERROE\n"); goto state; } else { double A = (W*Q-R*T)/(P*T-Q*Q); double B = (P*W-Q*R)/(Q*Q-P*T); double C = -(sum_x_2+sum_y_2+A*sum_x+B*sum_y)/n ; double r = 0.50*sqrt(A*A+B*B-4*C); printf("圆心坐标 :"); printf("(%f,%f)\n",A/2,B/2); printf("圆半径 : "); printf("%f\n",r) ; } system("pause"); return 0 ;}
2011年6月14日星期二
最小拟合圆
订阅:
博文评论 (Atom)
没有评论:
发表评论