圆上三点求圆心和半径

下面的程序实现用到C++和OpenCV。

先定义一个用于存储圆的数据的结构体CircleData:

struct CircleData
{
	Point2f center;
	int radius;
};

假设不共线的三点分别为pt1,pt2,pt3,通过三个点求圆的圆心center和半径radius有两个方法：

方法一：通过其中两点的中垂线交点求圆心

CircleData findCircle1(Point2f pt1, Point2f pt2, Point2f pt3)
{
	//定义两个点，分别表示两个中点
	Point2f midpt1, midpt2;
	//求出点1和点2的中点
	midpt1.x = (pt2.x + pt1.x) / 2;
	midpt1.y = (pt2.y + pt1.y) / 2;
	//求出点3和点1的中点
	midpt2.x = (pt3.x + pt1.x) / 2;
	midpt2.y = (pt3.y + pt1.y) / 2;
	//求出分别与直线pt1pt2，pt1pt3垂直的直线的斜率
	float k1 = -(pt2.x - pt1.x) / (pt2.y - pt1.y);
	float k2 = -(pt3.x - pt1.x) / (pt3.y - pt1.y);
	//然后求出过中点midpt1，斜率为k1的直线方程（既pt1pt2的中垂线）：y - midPt1.y = k1( x - midPt1.x)
	//以及过中点midpt2，斜率为k2的直线方程（既pt1pt3的中垂线）：y - midPt2.y = k2( x - midPt2.x)
	//定义一个圆的数据的结构体对象CD
	CircleData CD;
	//连立两条中垂线方程求解交点得到：
	CD.center.x = (midpt2.y - midpt1.y - k2* midpt2.x + k1*midpt1.x) / (k1 - k2);
	CD.center.y = midpt1.y + k1*(midpt2.y - midpt1.y - k2*midpt2.x + k2*midpt1.x) / (k1 - k2);
	//用圆心和其中一个点求距离得到半径：
	CD.radius = sqrtf((CD.center.x - pt1.x)*(CD.center.x - pt1.x) + (CD.center.y - pt1.y)*(CD.center.y - pt1.y));
	return CD;
}

方法二：

通过三个点到圆心距离相等建立方程：

(pt1.x-center.x)²-(pt1.y-center.y)²=radius²     式子(1)
(pt2.x-center.x)²-(pt2.y-center.y)²=radius²     式子(2)
(pt3.x-center.x)²-(pt3.y-center.y)²=radius²     式子(3)

式子(1)-式子(2)得：

pt1.x²+pt2.y²-pt1.y²-pt2.x²+2*center.x*pt2.x-2*center.x*pt1.x+2*center.y*pt1.y-2*center.y*pt2.y-=0

式子(2)-式子(3)得：

pt2.x²+pt3.y²-pt2.y²-pt3.x²+2*center.x*pt3.x-2*center.x*pt2.x+2*center.y*pt2.y-2*center.y*pt3.y-=0

整理上面的两个式子得到：

(2*pt2.x-2*pt1.x)*center.x+(2*pt1.y-2*pt2.y)*center.y=pt1.y²+pt2.x²-pt1.x²-pt2.y²
(2*pt3.x-2*pt2.x)*center.x+(2*pt2.y-2*pt3.y)*center.y=pt2.y²+pt3.x²-pt2.x²-pt3.y²

令：
A1=2*pt2.x-2*pt1.x      B1=2*pt1.y-2*pt2.y       C1=pt1.y²+pt2.x²-pt1.x²-pt2.y²
A2=2*pt3.x-2*pt2.x      B2=2*pt2.y-2*pt3.y       C2=pt2.y²+pt3.x²-pt2.x²-pt3.y²

则上述方程组变成一下形式：
A1*center.x+B1*center.y=C1；
A2*center.x+B2*center.y=C2
联立以上方程组可以求出：
center.x = (C1*B2 - C2*B1) / A1*B2 - A2*B1;
center.y = (A1*C2 - A2*C1) / A1*B2 - A2*B1;
（为了方便编写程序，令temp = A1*B2 - A2*B1）

具体程序实现函数为：

CircleData findCircle2(Point2f pt1, Point2f pt2, Point2f pt3)
{

	//令：
	//A1 = 2 * pt2.x - 2 * pt1.x      B1 = 2 * pt1.y - 2 * pt2.y       C1 = pt1.y² + pt2.x² - pt1.x² - pt2.y²
	//A2 = 2 * pt3.x - 2 * pt2.x      B2 = 2 * pt2.y - 2 * pt3.y       C2 = pt2.y² + pt3.x² - pt2.x² - pt3.y²
	float A1, A2, B1, B2, C1, C2, temp;
	A1 = pt1.x - pt2.x;
	B1 = pt1.y - pt2.y;
	C1 = (pow(pt1.x, 2) - pow(pt2.x, 2) + pow(pt1.y, 2) - pow(pt2.y, 2)) / 2;
	A2 = pt3.x - pt2.x;
	B2 = pt3.y - pt2.y;
	C2 = (pow(pt3.x, 2) - pow(pt2.x, 2) + pow(pt3.y, 2) - pow(pt2.y, 2)) / 2;
	//为了方便编写程序，令temp = A1*B2 - A2*B1
	temp = A1*B2 - A2*B1;
	//定义一个圆的数据的结构体对象CD
	CircleData CD;
	//判断三点是否共线
	if (temp == 0){
		//共线则将第一个点pt1作为圆心
		CD.center.x = pt1.x;
		CD.center.y = pt1.y;
	}
	else{
		//不共线则求出圆心：
		//center.x = (C1*B2 - C2*B1) / A1*B2 - A2*B1;
		//center.y = (A1*C2 - A2*C1) / A1*B2 - A2*B1;
		CD.center.x = (C1*B2 - C2*B1) / temp;
		CD.center.y = (A1*C2 - A2*C1) / temp;
	}
	
	CD.radius = sqrtf((CD.center.x - pt1.x)*(CD.center.x - pt1.x) + (CD.center.y - pt1.y)*(CD.center.y - pt1.y));
	return CD;
}

测试主程序为：

int _tmain(int argc, _TCHAR* argv[])
{
	//初始化一张图片，全部为黑色
	Mat dst(500, 500, CV_8UC3, Scalar(0, 0, 0));
	//随便定义三个点
	Point2f p1, p2, p3;
	p1.x = 150;
	p1.y = 150;
	p2.x = 200;
	p2.y = 200;
	p3.x = 250;
	p3.y = 100;
	//画出三个点
	circle(dst, p1, 2, Scalar(0, 0, 255), 1, 8, 0);
	circle(dst, p2, 2, Scalar(0, 0, 255), 1, 8, 0);
	circle(dst, p3, 2, Scalar(0, 0, 255), 1, 8, 0);
	//定义一个圆的数据结构体对象
	CircleData CD;
	//用第一种方法找圆
	CD = findCircle1(p1, p2, p3);
	circle(dst, CD.center, 2, Scalar(255, 0, 0), 1, 8, 0);//画出圆心
	circle(dst, CD.center, CD.radius, Scalar(0, 255, 0), 1, 8, 0);//画出圆
	imshow("circle", dst);
	waitKey();
	//用第二种方法找圆
	CD = findCircle2(p1, p2, p3);
	circle(dst, CD.center, 2, Scalar(255, 255, 0), 1, 8, 0);//画出圆心
	circle(dst, CD.center, CD.radius, Scalar(0, 255, 255), 1, 8, 0);//画出圆
	imshow("circle", dst);
	waitKey();
	return 0;
}

运行结果：

第一种方法：

    

第二种方法：