3点から円の中心と半径を求める
3点を通る円
円は一直線上ではない3点の座標があれば一意に決定します。
下図を参照してください。ここで、3点の座標を、
(x1, y1), (x2, y2), (x3, y3)
求める中心座標を、
(Cx, Cy)
求める半径を、
r
とします。
ごく普通に3つの連立方程式を解いていきます。
逆行列で方程式を解く
基本的には3つの連立方程式を一般的に解いてプログラム化すればよいのですが、できるだけ簡単なプログラムになるように工夫してみます。
[math]{ left( { x }_{ 1 }-c_{ x } right) }^{ 2 }+{ left( y_{ 1 }-c_{ y } right) }^{ 2 }={ r }^{ 2 }….(1)\ { left( { x }_{ 2 }-c_{ x } right) }^{ 2 }+{ left( y_{ 2 }-c_{ y } right) }^{ 2 }={ r }^{ 2 }….(2)\ { left( { x }_{ 3 }-c_{ x } right) }^{ 2 }+{ left( y_{ 3 }-c_{ y } right) }^{ 2 }={ r }^{ 2 }….(3)[/math](1)と(2)および(1)と(3)から以下が導けます
[math]c_{ x }left( 2cdot left( { x }_{ 2 }-{ x }_{ 1 } right) right) +{ c }_{ y }left( 2cdot left( y_{ 2 }-{ y }_{ 1 } right) right) ={ x }_{ 2 }^{ 2 }+{ y }_{ 2 }^{ 2 }-{ x }_{ 1 }^{ 2 }-{ y }_{ 1 }^{ 2 }…..(4)\ { c }_{ x }left( 2cdot left( { x }_{ 3 }-{ x }_{ 1 } right) right) +{ c }_{ y }left( 2cdot left( y_{ 3 }-{ y }_{ 1 } right) right) ={ x }_{ 3 }^{ 2 }+{ y }_{ 3 }^{ 2 }-{ x }_{ 1 }^{ 2 }-{ y }_{ 1 }^{ 2 }…..(5)[/math]ここで以下のように定義します。
（これがポイントです）
(4)(5)は上記定義と行列形式で以下のように記述できます。
[math]begin{bmatrix} 2a & 2b \ 2c & 2d end{bmatrix}begin{bmatrix} { c }_{ x } \ { c }_{ y } end{bmatrix}=begin{bmatrix} { x }_{ 2 }^{ 2 }+{ y }_{ 2 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 } \ { x }_{ 3 }^{ 2 }+{ y }_{ 3 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 } end{bmatrix}quad …..(6)[/math]ここで、Ax=bの形式なのでAの逆行列を求めればよいことになります。
Aのデターミナントは、
ですから、逆行列は、
[math]{ A }^{ -1 }=begin{bmatrix} frac { 2d }{ det(A) } & frac { -2b }{ det(A) } \ frac { -2c }{ det(A) } & frac { 2a }{ det(A) } end{bmatrix}quad =frac { 1 }{ 2(ad-bc) } begin{bmatrix} d & -b \ -c & a end{bmatrix}[/math]となります。従って解くべき、Cx,Cyは、
[math]begin{bmatrix} { c }_{ x } \ { c }_{ y } end{bmatrix}=frac { 1 }{ 2(ad-bc) } begin{bmatrix} d & -b \ -c & a end{bmatrix}begin{bmatrix} { x }_{ 2 }^{ 2 }+{ y }_{ 2 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 } \ { x }_{ 3 }^{ 2 }+{ y }_{ 3 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 } end{bmatrix}quad …..(7)[/math]となります。右辺のごちゃごちゃをa,b,c,dで整理してみると、
[math]{ x }_{ 2 }^{ 2 }+{ y }_{ 2 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 }\ =({ x }_{ 2 }-{ x }_{ 1 })({ x }_{ 2 }-{ x }_{ 1 }+2{ x }_{ 1 })+({ y }_{ 2 }-y_{ 1 })(y_{ 2 }-{ y }_{ 1 }+2{ y }_{ 1 })\ =a(a+2{ x }_{ 1 })+b(b+2y_{ 1 })\ ={ a }^{ 2 }+b^{ 2 }+2(a{ x }_{ 1 }+b{ y_{ 1 } })\ \ { x }_{ 3 }^{ 2 }+{ y }_{ 3 }^{ 2 }-{ x }_{ 1 }^{ 2 }-y_{ 1 }^{ 2 }\ ={ c }^{ 2 }+d^{ 2 }+2(c{ x }_{ 1 }+dy_{ 1 })[/math]となるので、これで、Cxを解いてみると、
[math]{ c }_{ x }=frac { d({ a }^{ 2 }+b^{ 2 })-b({ c }^{ 2 }+d^{ 2 }) }{ 2(ad-bc) } +frac { 1 }{ 2(ad-bc) } cdot (dcdot 2cdot (a{ x }_{ 1 }+b{ y }_{ 1 })-bcdot 2cdot (c{ x }_{ 1 }+d{ y }_{ 1 }))\ =frac { d({ a }^{ 2 }+b^{ 2 })-b({ c }^{ 2 }+d^{ 2 }) }{ 2(ad-bc) } +frac { 1 }{ ad-bc } cdot (ad{ x }_{ 1 }+bd{ y }_{ 1 }-bc{ x }_{ 1 }-bd{ y }_{ 1 })\ =frac { d({ a }^{ 2 }+b^{ 2 })-b({ c }^{ 2 }+d^{ 2 }) }{ 2(ad-bc) } +frac { 1 }{ ad-bc } { x }_{ 1 }(ad-bc)\ ={ frac { d({ a }^{ 2 }+b^{ 2 })-b({ c }^{ 2 }+d^{ 2 }) }{ 2(ad-bc) } +{ x }_{ 1 } }[/math]となるわけです。ポイントはごちゃごちゃしていた右側の部分がX1で置き換えられることですね。
中心座標のXがこれで計算できるようになりましたので、Yも簡単に求められます。
中心座標が求まれば、円周上点が既知なので半径は簡単です。
プログラム
では早速プログラムにしてみます。
中心Y座標を求めるところでちょっと工夫しているのに注目してください。
また、半径も誤差を考慮して、中心と既知の円周上3点全て求めて平均化しているところにも注目です。
//-----------------------------------------------------------------------------------
// ３点から円を作成する
//-----------------------------------------------------------------------------------
int calcCircleOf2Point(int x1, int y1, int x2, int y2, int x3, int y3, int *cx, int *cy, int *r)
{
long ox, oy, a, b, c, d ;
long r1, r2, r3 ;
int stat ;
a = x2 - x1 ;
b = y2 - y1 ;
c = x3 - x1 ;
d = y3 - y1 ;
stat = -1 ;
if ((a && d) || (b && c)) {
ox = x1 + (d * (a * a + b * b) - b * (c * c + d * d)) / (a * d - b * c) / 2 ;
if (b) {
oy = (a * (x1 + x2 - ox - ox) + b * (y1 + y2)) / b / 2 ;
} else {
oy = (c * (x1 + x3 - ox - ox) + d * (y1 + y3)) / d / 2 ;
}
r1 = sqrt((ox - x1) * (ox - x1) + (oy - y1) * (oy - y1)) ;
r2 = sqrt((ox - x2) * (ox - x2) + (oy - y2) * (oy - y2)) ;
r3 = sqrt((ox - x3) * (ox - x3) + (oy - y3) * (oy - y3)) ;
*cx = ox ;
*cy = oy ;
*r = (r1 + r2 + r3) / 3 ;
stat = 0 ;
}
return stat ;
}