不規則だが凸面体の体積を計算する方法をまとめようとしています。三角形分割を使用して多面体を複数のサブ四面体 (シンプレックス) に分割し、体積を個別に計算してから、すべてのサブボリューム値を合計します。
ただし、私のテストでは、以下のユニット - キューブに対して奇妙な結果が得られます。バグがどこにあるのか誰にも分かりませんか?
class Simplex(object):
def __init__(self,coordinates):
if not len(coordinates) == 4:
raise RuntimeError('You must provide only 4 coordinates!')
self.coordinates = coordinates
def volume(self):
'''
volume: Return volume of simplex. Formula from http://de.wikipedia.org/wiki/Tetraeder
'''
import numpy
vA = numpy.array(self.coordinates[1]) - numpy.array(self.coordinates[0])
vB = numpy.array(self.coordinates[2]) - numpy.array(self.coordinates[0])
vC = numpy.array(self.coordinates[3]) - numpy.array(self.coordinates[0])
return numpy.abs(numpy.dot(numpy.cross(vA,vB),vC)) / 6.0
'''
Old code that did not work
class Polyeder(object):
def __init__(self,coordinates):
if len(coordinates) < 4:
raise RuntimeError('You must provide at least 4 coordinates!')
self.coordinates = coordinates
def volume(self):
pivotCoordinate = self.coordinates[0]
volumeSum = 0
for i in xrange(1,len(self.coordinates)-3):
newCoordinates = [pivotCoordinate]
for j in xrange(i,i+3):
newCoordinates.append(self.coordinates[j])
simplex = Simplex(newCoordinates)
volumeSum += simplex.volume()
return volumeSum
'''
class Polyeder(object):
def __init__(self,coordinates):
'''
Constructor
'''
if len(coordinates) < 4:
raise RuntimeError('You must provide at least 4 coordinates!')
self.coordinates = coordinates
def volume(self):
from pyhull.delaunay import DelaunayTri
delaunay = DelaunayTri(self.coordinates,joggle=True)
volume = 0
for vertices in delaunay.vertices:
coords = [self.coordinates[i] for i in vertices]
simplex = Simplex(coords)
volume += simplex.volume()
return volume
coords = []
coords.append([0,0,0])
coords.append([1,0,0])
coords.append([0,1,0])
coords.append([0,0,1])
s = Simplex(coords)
print s.volume()
coords.append([0,1,1])
coords.append([1,0,1])
coords.append([1,1,0])
coords.append([1,1,1])
p = Polyeder(coords)
print p.volume()
古い結果の出力は次のとおりです。
0.166666666667
0.666666666667
正四面体の値は 1/6 (正しい) ですが、単位立方体の値は 1 です。
新しい結果: 0.166666666667 1.0