2011年6月14日星期二

  最小拟合圆

#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 ;}

没有评论:

发表评论