三角形外接圆圆心

时间：2011-06-06 来源：LitIce

给出三个点求圆心，我之前有个很傻的做法。
分别求出两条边的中垂线，然后求交点。
为什么我会这么想到这里去呢？
而不是想到直接列三个方程解方程这么简单的做法呢？
我经过认真地思考原因就是高中老师告诉我们解圆心坐标应该求中垂线。
也的确
中点坐标直接求出，然后点斜式直接求出垂线方程。
然后解两个直线的交点。
这样可以简化计算复杂度
但是！！！
现在我是编程解答~计算交给计算机。

同一道题POJ1329

用前者的方法我打了100+行的代码。(当然也可能只是我太挫)
但是用后者我只用了60行代码。

不同时候不同条件用不同方法~~

point circumcenter (point p0 , point p1 , point p2)
{
point ret;
double a1=p1.x-p0.x,b1 = p1.y - p0.y,c1 = (sqr(a1) + sqr(b1)) / 2;
  double a2=p2.x-p0.x,b2 = p2.y - p0.y,c2 = (sqr(a2) + sqr(b2)) / 2;
  double d = a1 * b2 - a2 * b1;
  ret.x = p0.x + (c1 * b2 - c2 * b1) / d;
  ret.y = p0.y + (a1 * c2 - a2 * c1) / d;
return ret;
}
改函数读入p0,p1,p2三个点返回该三角形外接圆圆心。

列出方程
sqr(O.x-p0.x)+sqr(O.y-p0.y) = sqr(r); (1)
sqr(O.x-p0.x)+sqr(O.y-p0.y) = sqr(r); (2)
sqr(O.x-p0.x)+sqr(O.y-p0.y) = sqr(r); (3)
联立(1)(2)得
(p1.x-p0.x)*O.x+(p1.y-p0.y)*O.y = (sqr(p1.x)+sqr(p1.y)-sqr(p0.x)-sqr(p0.y))/2; (4)
联立(1)(3)的
(p2.x-p0.x)*O.x+(p2.y-p0.y)*O.y = (sqr(p2.x)+sqr(p2.y)-sqr(p0.x)-sqr(p0.y))/2; (5)
令O.x = p0.x+x;
O.y = p0.y+y;
带入(4)(5)后化简得到
(p1.x-p0.x)*x+(p1.y-p0.y)*y = (sqr(p1.x-p0.x)+sqr(p1.y-p0.y))/2;
(p2.x-p0.x)*x+(p2.y-p0.y)*y = (sqr(p2.x-p0.x)+sqr(p2.y-p0.y))/2;
分别解得(x,y)
p0+(x,y) = O;

#include<iostream>
#include<cstdio>
#include<cmath>

using namespace std;

struct  point{
        double x,y;
} p0,p1,p2,ret;

double sqr(double x){
        return x*x;
}
double r;

point circumcenter (point p0 , point p1 , point p2)  
{ 
        point ret;
        double a1=p1.x-p0.x,b1 = p1.y - p0.y,c1 = (sqr(a1) + sqr(b1)) / 2; 
        double a2=p2.x-p0.x,b2 = p2.y - p0.y,c2 = (sqr(a2) + sqr(b2)) / 2; 
        double d = a1 * b2 - a2 * b1; 
        ret.x = p0.x + (c1 * b2 - c2 * b1) / d; 
        ret.y = p0.y + (a1 * c2  - a2 * c1) / d; 
        return ret;
}  
void out(double x){
        if (x<0){
                printf(" + %.3lf",-x);
        }else
        if (x>0){
                printf(" - %.3lf",x);
        }
}
int main(){
        freopen("in.txt","r",stdin);
        freopen("out1.txt","w",stdout);
        int Test=0;
        while (cin >> p0.x >> p0.y >> p1.x >> p1.y >> p2.x >> p2.y){
        //      ++Test;
                //if (Test>1) 
        
                ret = circumcenter( p0,p1,p2 );
                r = sqrt( sqr(ret.x-p0.x) + sqr(ret.y-p0.y) );
             printf("(x");
             out(ret.x);
             printf(")^2 + (y");
             out(ret.y);
             printf(")^2 = %.3lf^2\n",r);
             printf("x^2 + y^2");
             out( ret.x*2 );
             printf("x");
             out( ret.y*2 );
             printf("y");
             out( sqr(r)-sqr(ret.x)-sqr(ret.y) );
             printf(" = 0\n");
                printf("\n");
        }
        return 0;
}