空間上の四面体の体積

質問

Ｃ言語で空間上の四面体の体積を求めるプログラムを作りたいんですが、どうすればいいのかわかりません。
構造体を使って空間上（三次元）の４点Ａ，Ｂ，Ｃ，Ｄの座標を定義するのですが。
四面体の求め方もわかりません。ヘロンの公式じゃ求められませんよね（＾＾；
空間上の三角形の面積を求めるプログラムを作ってみました。
少しでも参考になれば幸いです。
よろしくお願いします。

#include <stdio.h>
#include <math.h>
int main(void)
{
double AB, BC, CA, s, m;

struct zahyo
{
double x;
double y;
double z;
}A, B, C;

scanf("%lf %lf %lf",&A.x,&A.y,&A.z);
scanf("%lf %lf %lf",&B.x,&B.y,&B.z);
scanf("%lf %lf %lf",&C.x,&C.y,&C.z);

AB=sqrt((A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z));
BC=sqrt((B.x-C.x)*(B.x-C.x)+(B.y-C.y)*(B.y-C.y)+(B.z-C.z)*(B.z-C.z));
CA=sqrt((C.x-A.x)*(C.x-A.x)+(C.y-A.y)*(C.y-A.y)+(C.z-A.z)*(C.z-A.z));

printf("AB=%lf\nBC=%lf\nCA=%lf\n",AB,BC,CA);
s=(AB+BC+CA)/2;
m=sqrt(s*(s-AB)*(s-BC)*(s-CA));
printf("m=%lf\n",m);

return 0;
}

回答No.4

noname#5537 · ベストアンサー

> (1/6)*(AB*AC)*ADを率直に求めて実行したのですが、結果がでません。

#1 のはじめに少し書いたのですが（分かりにくかったかもしれません），
AB は A から B に向かうベクトルです。AC，AD についても同様です。
(AB×AC)・AD という風に '×' と '・' を使い分けているのには意味があります。

ベクトルの計算ですので C 言語でプログラムを書く場合，成分ごとに計算してやらないといけません。

今，AB の x 成分，y 成分，z 成分をそれぞれ，ABx, ABy, ABz としますと，
　ABx = (B の x 座標) - (A の x 座標)
　ABy = (B の y 座標) - (A の y 座標)
　ABz = (B の z 座標) - (A の z 座標)
です。AC，AD についても同様です。

このとき，
　(AB×AC)・AD = (ABy*ACz - ABz*ACy)*ADx
　　　　　　　+ (ABz*ACx - ABx*ACz)*ADy
　　　　　　　+ (ABx*ACy - ABy*ACx)*ADz
です。

これを使ってプログラムを書いてみてください。

お礼

再び回答ありがとうございます。
＞ABx = (B の x 座標) - (A の x 座標)
　ABy = (B の y 座標) - (A の y 座標)
　ABz = (B の z 座標) - (A の z 座標)
のそれぞれの成分を出すのを忘れていました。
たぶんこのまま出来ると思います。
アドバイスありがとうございました。

回答No.3

neKo_deux

別の方法の案です。

厳密解で無くても良いという条件で、せっかくコンピュータで計算するのなら、モンテカルロ法のような方法を使ってみては？

簡単化のため、ＡＢＣＤ点のｘｙｚの範囲を０～１とします。
乱数を使って０～１の中の適当な点を決めます。
その点が立体ＡＢＣＤの内側であるかどうか判定。
面ＢＣＤに対して、点Ａと同じ側か？
面ＡＣＤに対して、点Ｂと同じ側か？
面ＡＢＤに対して、点Ｃと同じ側か？
面ＡＢＣに対して、点Ｄと同じ側か？
ならは、立体の内側であるとしてhit回数をカウント。
以上をtry回数だけ繰り返して、hit/tryを立体の体積の近似値とする。

ある程度の精度を期待するのなら、それなりの回数繰り返す必要がありますが、プログラムは作りやすいと思います。
最近のＰＣだと、かなりのマシンパワーが稼げるので…。

--
モンテカルロ法で円周率を求めよう
http://www.f.waseda.jp/takezawa/math/Pi/monte.html

参考URL：

http://www.f.waseda.jp/takezawa/math/Pi/monte.html

お礼

回答ありがとうございます。
モンテカルロ法というのを初めて知りました。
参考ＵＲＬを見てもちょっと難しいですね（汗）
余力があればこちらの方法を使ってプログラミングしてみたいと思います。
アドバイスありがとうございました。

回答No.2

noname#5537

少し補足です。

スカラー三重積を計算するときは，AB と AC の外積を計算してから，AD との内積を取るよりも，
AB, AC, AD の x, y, z 成分から直接，スカラー三重積を計算するほうがプログラムしやすいと思います。
http://www5.plala.or.jp/h-fuchi/math/vector/4.htm

参考URL：

http://www5.plala.or.jp/h-fuchi/math/vector/4.htm

お礼

２回も回答して頂きありがとうございます。
(1/6)*(AB*AC)*ADを率直に求めて実行したのですが、結果がでません。
ちゃんとprintfはしてるんですが。
たぶん自分のしょーもいミスだと思うんですけど。
ＵＲＬはとても参考になりました。
スカラー三重積なんてあるんですね。
知りませんでした。
回答ありがとうございました。

回答No.1

noname#5537

三角錐の体積なので，底面積×高さ÷3 で求められます。

以下，AB，AC，AD はベクトルを表しているものとします。

まず ABC を底面とすると底面積は，
　(1/2)*|AB×AC|
で求められます。
ここで，AB×AC は AB と AC の外積です。

次に，面ABCの単位法線ベクトルを n とすると，
　n = (AB×AC)/|AB×AC|
であり，三角錐の高さは，
　n・AD　（n と AD の内積)
で求められます。

以上をまとめると三角錐の体積は，
　(1/6)*(AB×AC)・AD
で求められます。

要するに AB, AC, AD のスカラー三重積の 1/6 です。
http://spinman.phys.metro-u.ac.jp/lecture/Mech/Vector/Vector.html

A，B，C，D の位置関係によってはマイナスがつくかもしれないのでそのときは絶対値をとってください。

以上を素直にプログラムすればいいと思います。

参考URL：

http://spinman.phys.metro-u.ac.jp/lecture/Mech/Vector/Vector.html