物体の体積と表面積

R cpp Rpackage 形・動き 数学

物体の体積と表面積を計算したい。
三角メッシュ化されていれば、表面の三角形を足せば表面積はどうにかなる。
問題は体積。
　
The GNU Triangulated Surface Library(GTS)というもので計算できるっぽい。

sudo apt-get install libgts-dev

　
Python でも使えるようになっている。しかし、pip でインストールできなかったのでソースからインストールする。
　
例によって説明など存在しない。
　
GTS が使うのはgts ファイルという、こんな感じのファイル。
一辺10の正四面体の頂点、辺、三角形の数と
座標
エッジグラフ
エッジで構成される三角形
をひたすら並べている。

8 18 12
 0  0  0
10  0  0
10 10  0
 0 10  0
 0  0 10
10  0 10
10 10 10
 0 10 10
1 2
2 3
3 4
4 1
5 6
6 7
7 8
8 5
1 5
2 6
3 7
4 8
1 6
2 7
3 8
4 5
1 3
5 7
13 1 10
13 5 9
14 2 11
14 6 10
15 3 12
15 7 11
16 4 9
16 8 12
17 2 1
17 4 3
18 5 6
18 7 8

　
これをpython gts で体積、表面積の計算をする。

# Python 2.7
import gts
f = "cube.gts"
g = open(f)
s = gts.read(g) # gts 読み込み

s.volume() # 体積
s.area()   # 表面積

　
問題はgts ファイルというやつで、obj ファイルからの変換がめんどくさい。
gtsの変換プログラム作ったよというメーリスも動かない。
　
Rvcg パッケージには、vcgGetEdge(mesh3d obj) という関数があって、三角形エッジをグラフにしてくれるが、これを使って上のcubeをいじってみると

library(Rvcg)
e <- vcgGetEdge(mesh.tri)
edge <- mapply(function(x) which(sapply(apply(t(mesh.tri$it), 1, setdiff, e[x, 1:2]), length) == 1), seq(nrow(e)))


r <- vector("list", max(edge))
for(i in seq(ncol(edge))){
  r[edge[,i]] <- lapply(r[edge[,i]], c, i)
}
r <- do.call(rbind.data.frame, r)
colnames(r) <- NULL

こんなファイルになる。番号が少し入れ替わる。

8 18 12
 0  0  0
10  0  0
10 10  0
 0 10  0
 0  0 10
10  0 10
10 10 10
 0 10 10
1 2
1 3
1 4
1 5
1 6
2 3
2 6
2 7
3 4
3 7
3 8
4 5
4 8
5 6
5 7
5 8
6 7
7 8
1 5 7
4 5 14
6 8 10
7 8 17
9 11 13
10 11 18
3 4 12
12 13 16
1 2 6
2 3 9
14 15 17
15 16 18

　
これを使ってpython gts をやってみると、area() はできるが volume() で

RuntimeError: Surface is not orientable

と言われてエラー。
　
それでは世の中に転がっているものから obj → gts 変換をしようとして、Hira 3D viewerというものでobj ファイル読み込みからのgts 出力をしてみるのだが、gts で頂点が重複して出力されるのでどうにもならん。
だがしかし、GUI で体積は計算できる。
　
がんばって変換スクリプト書くしかなさそう。
　
と思っていたら計算できるっぽい。これとかこれとか。
4点(0,0,0), (x_1, y_1, z_1), (x_2, y_2, z_2), (x_3, y_3, z_3)で作られる四面体をひたすら足せばいいっぽくて、この物体の体積は
\begin{vmatrix}x_1& y_1& z_1 \\x_2& y_2& z_2\\x_3& y_3& z_3 \end{vmatrix}/6
と、行列式で求められる。

sum(mapply(function(x) det(bun.v[bun.f[x,],])/6, seq(nrow(bun.f))))