三次元空間で任意の三点を通る円の中心Oの位置座標Ocに関する式を考えます。
辺cの中点の位置座標は(Tp1+Tp2)/2で、三角形AOBは、OAとOBの辺の長さが円の半径Rとなる二等辺三角形なので、ベクトル(Tp1-Tp2)とベクトル(Oc-(Tp1+Tp2)/2)は直交します。
よってその内積は、(Tp1-Tp2)・(Oc-(Tp1+Tp2)/2)=0
同じように、(Tp2-Tp3)・(Oc-(Tp2+Tp3)/2)=0、(Tp3-Tp1)・(Oc-(Tp3+Tp1)/2)=0 です。
また、法線ベクトルNvと(Tp1-Oc)も直交しているので、Nv・(Tp1-Oc)=0 同じように、Nv・(Tp2-Oc)=0、Nv・(Tp3-Oc)=0 です。
独立した式を三つ選んで三元連立方程式にします。
(Tp3-Tp1)・(Oc-(Tp3+Tp1)/2)=0
Nv・(Tp3-Oc)=0
(Tp1-Tp2)・(Oc-(Tp1+Tp2)/2)=0
この三元連立方程式の解をクラメルの公式により求めれば、円の中心Oの位置座標Ocが分かります。
一つ目の式を成分で表すと
(Tp3-Tp1)・Oc=(Tp3-Tp1)・(Tp3+Tp1)/2
(Tp3.x-Tp1.x)Oc.x+(Tp3.y-Tp1.y)Oc.y+(Tp3.z-Tp1.z)Oc.z=S1/2
S1=pow(Tp3.x,2)-pow(Tp1.x,2)+pow(Tp3.y,2)-pow(Tp1.y,2)+pow(Tp3.z,2)-pow(Tp1.z,2)
円の中心座標Ocに対する項をベクトルで表現すると
K1=<Tp3.x-Tp1.x,Tp3.y-Tp1.y,Tp3.z-Tp1.z,S1/2>
二つ目の式を成分で表すと、
Nv・Oc=Nv・Tp3
Nv.x*Oc.x+Nv.y*Oc.y+Nv.z*Oc.z=Nv.x* Tp3.x +Nv.y*Tp3.y+Nv.z*Tp3.z
円の中心座標Ocに対する項をベクトルで表現すると
K2=<Nv.x,Nv.y,Nv.z,Nv.x*Tp3.x+Nv.y*Tp3.y+Nv.z*Tp3.z>
三つ目の式を成分で表すと
(Tp1-Tp2)・Oc=(Tp1-Tp2)・(Tp1+Tp2)/2
(Tp1.x-Tp2.x)Oc.x+(Tp1.y-Tp2.y)Oc.y+(Tp1.z-Tp2.z)Oc.z=S3/2
S3=pow(Tp1.x,2)-pow(Tp2.x,2)+pow(Tp1.y,2)-pow(Tp2.y,2)+pow(Tp1.z,2)-pow(Tp2.z,2)
円の中心座標Ocに対する項をベクトルで表現すると
K3=<Tp1.x-Tp2.x,Tp1.y-Tp2.y,Tp1.z-Tp2.z,S3/2>
ベクトルで表現したK1,K2,K3を用いて、クラメルの公式により、円の中心座標Ocを求めます。
参考までにPOV-Rayで作成した、クラメルの公式を用いて円の中心座標を求めるマクロを載せておきます。
//**********************************************************************************
//input parameter Tp1, Tp2, Tp3
#macro Circle_center(Tp1,Tp2,Tp3)
#local Nv=Normal_vector(Tp2,Tp1,Tp3); //normal vector of the outside;
//determinant
#local S1=pow(Tp3.x,2)-pow(Tp1.x,2)+pow(Tp3.y,2)-pow(Tp1.y,2)+pow(Tp3.z,2)-pow(Tp1.z,2);
#local K1=<Tp3.x-Tp1.x,Tp3.y-Tp1.y,Tp3.z-Tp1.z,S1/2>;
#local K2=<Nv.x,Nv.y,Nv.z,Nv.x*Tp3.x+Nv.y*Tp3.y+Nv.z*Tp3.z>;
#local S3=pow(Tp1.x,2)-pow(Tp2.x,2)+pow(Tp1.y,2)-pow(Tp2.y,2)+pow(Tp1.z,2)-pow(Tp2.z,2);
#local K3=<Tp1.x-Tp2.x,Tp1.y-Tp2.y,Tp1.z-Tp2.z,S3/2>;
//cramer
#local D =(K1.x*K2.y*K3.z)+(K1.y*K2.z*K3.x)+(K1.z*K2.x*K3.y)-(K3.x*K2.y*K1.z)-(K1.x*K2.z*K3.y)-(K2.x*K1.y*K3.z);
#local D1=(K1.t*K2.y*K3.z)+(K1.y*K2.z*K3.t)+(K1.z*K2.t*K3.y)-(K3.t*K2.y*K1.z)-(K1.t*K2.z*K3.y)-(K2.t*K1.y*K3.z);
#local D2=(K1.x*K2.t*K3.z)+(K1.t*K2.z*K3.x)+(K1.z*K2.x*K3.t)-(K3.x*K2.t*K1.z)-(K1.x*K2.z*K3.t)-(K2.x*K1.t*K3.z);
#local D3=(K1.x*K2.y*K3.t)+(K1.y*K2.t*K3.x)+(K1.t*K2.x*K3.y)-(K3.x*K2.y*K1.t)-(K1.x*K2.t*K3.y)-(K2.x*K1.y*K3.t);
#local Oc=<D1/D,D2/D,D3/D>; //center
Oc // output parameter
#end
//**********************************************************************************
