任意の3点を通る球体の中心の求め方

質問

任意の3点、点Ａ(x1,y1,z1)、点Ｂ(x2,y2,z2)、点Ｃ(x3,y3,z3)を通る球体半径ｒの中心点Ｄ（a,b,c)を求める計算方式を教えてください。

マイナスとプラスで出てくると思いますが最終的にプラスの値を採用とします。

エクセルで常に計算出来るように書式を作りたいです。

ソルバーを使ってもなぜかうまくいきません。

私の計算式が間違っているのかもしれませんが。

http://detail.chiebukuro.yahoo.co.jp/qa/question_detail/q1230122222

回答No.7

noname#230359

回答(6)です。

下のほうで求めた中心座標の検証に使った、CAD(Rhinoceros)での求めかたの例も書いておきます。

1)3点を描く。
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100217062901.jpg

2)3点と中心とする球を描く。
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100217062916.jpg

3)球の交線の交点が答え。
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100217062932.jpg

簡単に求められて、確認も目視で分かりやすいですが、CADが高価だった一昔前でしたらとても贅沢な手段ですね。

Excelのソルバーは使ったことがなかったので、私も勉強になりました。
ついでに、任意の4点を通る球の中心と半径も、同じ方法で求められることを確認しました。
きっかけをいただき、ありがとうございました。

お礼

返答遅くなりました。回答何度もありがとうございます。
ソルバーまで検証して頂きありがとうございます。
私は2Ｄ-ＣＡＤしかなく3Ｄ検証が出来なく困っていました。

回答No.6

noname#230359

回答(5)です。

>ソルバーを使ってもなぜかうまくいきません。

ソルバーも試してみました。

質問で紹介された参考URLの内容そのままではダメで、
制約条件を設定したらうまくいきました。

仮に、中心に3点の最小値を入れて計算させると、
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100216221739.jpg

中心座標1が求まります。
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100216221804.jpg

同様に、中心に3点の最大値を入れて計算させると、中心座標2が求まります。
さすがExcel、すばらしいですね。

回答No.5

noname#230359

回答(3)です。回答(1)さんへ、

>青軸上ならどこでもかける
>
>任意の3点を通る最小の球体の中心　って縛りをつけないと

球の半径rを指定するので2つ以内に決まりますよ。
こんな感じ。
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100216073946.jpg

回答No.4

noname#230359

もう、殆ど答えは出ていますが、中心点Ｄ（a,b,c)で半径ｒの円を考えて、
点Ａを通る場合と点Ｂを通る場合と、点Ｃを通る場合とで式をたてると、
求まります。

お礼

返答遅くなりました。回答ありがとうございます。
計算するとなると大変ですね。

回答No.3

noname#230359

下記URLで、回答(2)さんの内容からもうちょっと先のところまでお話がされています。

(BIGLOBEなんでも相談室)
http://soudan1.biglobe.ne.jp/qa195295.html

Perlで確認してみました。
-----
# 3点
my($x1,$y1,$z1)=(-51,-21,-11);
my($x2,$y2,$z2)=( 62, 12, 22);
my($x3,$y3,$z3)=(-73, 83, 33);
# 半径
my $r1=100;

my $tm01=$x1**2-$x2**2+$y1**2-$y2**2+$z1**2-$z2**2;
my $tm02=$x1**2-$x3**2+$y1**2-$y3**2+$z1**2-$z3**2;
my $tm11=-2*($x1-$x2)*($z1-$z3)+2*($x1-$x3)*($z1-$z2);
my $tm12=-2*($y1-$y2)*($z1-$z3)+2*($y1-$y3)*($z1-$z2);
my $tm13=$tm01*($z1-$z3)-$tm02*($z1-$z2);
my $tm21=-2*($x1-$x2)*($y1-$y3)+2*($x1-$x3)*($y1-$y2);
my $tm22=-2*($z1-$z2)*($y1-$y3)+2*($z1-$z3)*($y1-$y2);
my $tm23=$tm01*($y1-$y3)-$tm02*($y1-$y2);
my $tma=1+$tm11**2/$tm12**2+$tm21**2/$tm22**2;
my $tmb=-2*$x1+2*($y1+$tm13/$tm12)*$tm11/$tm12+2*($z1+$tm23/$tm22)*$tm21/$tm22;
my $tmc=$x1**2+($y1+$tm13/$tm12)**2+($z1+$tm23/$tm22)**2-$r1**2;
my $xq1=(-$tmb+sqrt($tmb**2-4*$tma*$tmc))/2/$tma;
my $xq2=(-$tmb-sqrt($tmb**2-4*$tma*$tmc))/2/$tma;
my $yq1=-$tm13/$tm12-$tm11/$tm12*$xq1;
my $yq2=-$tm13/$tm12-$tm11/$tm12*$xq2;
my $zq1=-$tm23/$tm22-$tm21/$tm22*$xq1;
my $zq2=-$tm23/$tm22-$tm21/$tm22*$xq2;

printf("xq1=%9.3f yq1=%9.3f yq1=%9.3f\n",$xq1,$yq1,$zq1);
printf("xq2=%9.3f yq2=%9.3f yq2=%9.3f\n",$xq2,$yq2,$zq2);
-----

出た答え。
-----
xq1= -2.410 yq1= 63.052 yq1= -34.965
xq2= -20.432 yq2= 11.190 yq2= 78.607
-----

CADで作図したものと一致したので大丈夫と思います。

お疲れさまでした。

すみません、間違いがありました。
後ろのほうです。

+++++ 誤 +++++
printf("xq1=%9.3f yq1=%9.3f yq1=%9.3f\n",$xq1,$yq1,$zq1);
printf("xq2=%9.3f yq2=%9.3f yq2=%9.3f\n",$xq2,$yq2,$zq2);
-----

出た答え。
-----
xq1= -2.410 yq1= 63.052 yq1= -34.965
xq2= -20.432 yq2= 11.190 yq2= 78.607
-----
++++++++++++++

+++++ 正 +++++
printf("xq1=%9.3f yq1=%9.3f zq1=%9.3f\n",$xq1,$yq1,$zq1);
printf("xq2=%9.3f yq2=%9.3f zq2=%9.3f\n",$xq2,$yq2,$zq2);
-----

出た答え。
-----
xq1= -2.410 yq1= 63.052 zq1= -34.965
xq2= -20.432 yq2= 11.190 zq2= 78.607
-----
++++++++++++++

参考URL：

http://soudan1.biglobe.ne.jp/qa195295.html

回答No.2

noname#230359

問題の条件から次の三式が得られます。
(x1-a)^2+(y1-b)^2+(z1-c)^2=r^2
(x2-a)^2+(y2-b)^2+(z2-c)^2=r^2
(x3-a)^2+(y3-b)^2+(z3-c)^2=r^2
変数の二乗の項が含まれていますが、それぞれの式を展開して、引き算をすれば、二乗の項は消えて、a，b，cを変数とする連立三元一次方程式になりますので、これを解けば球体の中心座標を求めることができます。ここまでは確認しましたが、その後は面倒くさくて挑戦する気になりませんでした。

お礼

返答遅くなりました。回答ありがとうございます。
この連立三次方程式を解くのが自力では大変ですよね。

補足

3元2次連立方程式ですね。間違えていました。

回答No.1

noname#230359

点A　点B　をX軸　とし　点C　で　Yを通る方向を求め　新たな座標系を作る
そぞ座標系で　点A　点B　点Cの値を変換する

http://www.geocities.co.jp/Technopolis-Mars/8897/FIG/313r/3point_r.htm
で　中心O　を求め

元の座標系に戻す

うーん、考えただけでめんどくさそうだ

http://www.app-pc-soft.jp/file10_10.html
絵で書くと簡単だけどね

http://www.app-pc-soft.jp/file10_10.html
こっちから途中まで解いてみた

X-Y平面上だけど　あとで変換する
点A（XA,YA,0)
点B（XB,YB,0)
点C（XC,YC,0)

点Aを原点とする
点A'（A'X= 0 , A'Y= 0 , 0)
点B'（B'X=XB-XA , B'Y=YB-YA , 0)
点C'（C'X=XC-XA , C'Y=YC-YA , 0)

直行時の点Bの角度
BΘ=ATAN(B'Y/B'X)*180/PI()

直行時の点Cの角度
CΘ=ATAN(C'Y/C'X)*180/PI()

Z軸をBΘ回転したときの
点C　の　角度　ちょっと（符号が怪しい）
CBΘ　=　CΘ　- BΘ　

点A-C　の距離
R　＝　sqrt(C'X^2＋C'Y~2)

Z軸をBΘ回転したときの
座標系
点A''(A''X= 0 , A''Y=0 ,A''Z = 0）
点B''(B''X=1/COS(BΘ *PI()/180)*B'X , B''Y=0 ,B''Z = 0)
点C''(C''X= COS(CBΘ*PI()/180)*R , C''Y=SIN(CBΘ*PI()/180)*R,C''Z=0)

点A''-点B''の中点を通り　点A''-点B''に直行
X　=　B''X　/　2

ここまで解いた
あとは　点A''-点C''の中点を通り　点A''-点C''の直行式を求め
交わるところが中心

２D平面上だがあとで３D変換すれば求められるはず
検算してないけど

あとは　点A''-点C''の中点を通るグラフの式
Y=(C''Y-A''Y)/(C''X-A''X) * X

その角度
Θ＝ATAN(Y)*180/PI()

それに直角で交わる角度
Θ１＝Θ-90

Θ１の傾きを持ち点A''-点C''の中点を通るグラフの　a

a　=　TAN(B20*PI()/180)

そのb
b=(C''Y/2) - a * (C''X/2)

Θ１の傾きを持ち点A''-点C''の中点を通るグラフの式

y=ax+b

点A''-点B''の中点を通り　点A''-点B''に直行
X　=　B''X　/　2

と
Θ１の傾きを持ち点A''-点C''の中点を通るグラフの式
y=ax+b

を交わる　点を出せば　平面上の　中心点が求められる

で
はじめに戻って
ワールドから　三点の平面を持つ座標の変換をして

今の計算をして

結果を　ワールド座標に変換すれば出来上がり
検算まだだけどね

>>その角度
>>Θ＝ATAN(Y)*180/PI()

この部分で90度以上になると成り立たないや

点ABCを
はじめに　∠AB－AC　が　鈍角になる
AX<BX
とか
決めないといけないみたい

決めれば成り立つ

ふと思ったんだが
任意の3点を通る球体の中心
球体だと解けないぞ

上記で点O（3点を通る円の中心）
求め　その点O　を原点とする　座標系を　作る

点A（XA,YA,0)
点B（XB,YB,0)
点C（XC,YC,0)

点O(0,0,0)

点Oを中心に球をかけるが
点(0,0,100）でも点ABCを通る球は書けてしまう
Z 軸上にどの点でも書けてしまうよ

こんな感じ
http://mcnc.hp.infoseek.co.jp/cgi-bin/img-box/img20100215205559.jpg

青軸上ならどこでもかける

任意の3点を通る最小の球体の中心　って縛りをつけないと

お礼

返答遅くなりました。回答ありがとうございます。
3Ｄソフトがあると便利ですね。私はもっていなく2Ｄ-ＣＡＤでの検討となっています。