On 01/02/2015 at 05:06, xxxxxxxx wrote:
Hello,
a few questions and answers referring also to this logged thread:
https://plugincafe.maxon.net/topic/8445/11017_3x3-matrix
(the test scripts below shows how to take advantage of c4d.Matrix while calculating 3x3Matrices. / The Intersect_Mp function)
Testing what will be the fastest approach to get triangle ray intersections with Python in c4d, I compared several algorithms: With the first attached script 16 million polygons and the second 16 million rays and some questions pop up:
1.The results for 16 million rays lead me to the assumption, that c4d also uses the Möller-Trumbore algorithm, because c4d.RayCollider is as fast as the Möller-Trumbore algorithm. Additionally, one can read barycentric coordinates from c4d.RayCollider which are calculated by Möller-Trumbore, too.
Is this true?
2.For 16 million polygons c4d.RayCollider is by far the fastest one, which lets me think that it uses an octree for reducing the areas to calculate.
Is this true?
3.For rebuilding an octree in Python what will be the best way of storage?
c4d.storage.ByteSeq()?
4.What is the best way to set pointers to polygons in parent nodes in an octree with python in c4d?
Any help is very much appreciated!
Thanks
Martin
test scripts:
The first one is for testing a lot of polygons:
you´ll need a Spline object which is called LINE and object for collision.
For the second one an empty scene is enough but the resolution of the ray grid should be adjusted with the grid variables in main.
Have fun sending rays!!
1___
import c4d,time
from c4d import utils
#############################################
def Triangulate(op) :
#convert and triangulate
doc = op.GetDocument()
if op.GetType() !=c4d.Opolygon:
virtualop = utils.SendModelingCommand(command=c4d.MCOMMAND_CURRENTSTATETOOBJECT,
list=[op.GetClone()],
doc=doc)
if not virtualop: return #it was not possible to convert the object
else:
virtualop = [op.GetClone()]
obj = utils.SendModelingCommand(command=c4d.MCOMMAND_TRIANGULATE,
list=[virtualop[0]],
doc=doc)
if not obj: return #it was not possible to triangulate the object
return virtualop[0]
#############################################
def GlobalPoints(op) :
#translate points to global coordinates
matr = op.GetMg()
pcnt = op.GetPointCount()
points = op.GetAllPoints()
for i in xrange(pcnt) :
points[i] = op.GetPoint(i)*matr
return points
#############################################
def Intersect_Mt(a,b,c,p0,p1,epsi) :
case = 0
result = None
ldir = p1-p0
##check if parallel
e1 = b - a
e2 = c - a
s1 = ldir % e2
divisor = s1.Dot(e1)
if 0.0-epsi <= divisor <= 0.0+epsi:
#print "case1/ parallel to triangle"
case = 1
return
#--------------------------------------------
## Compute barycentric coordinates
inv_divisor = 1.0 / divisor
d = p0 - a
b1 = d.Dot(s1) * inv_divisor
if (b1 < 0.0 or b1 > 1) :
#print "case2/ to high or to low"
case = 2
#--------------------------------------------
s2 = d % e1
b2 = ldir.Dot(s2) * inv_divisor
if (b2 < 0.0 or (b1 + b2) > 1.0) :
#print "case3/b2 < to left or to right"
case = 3
#--------------------------------------------
## Compute distance(p0-intersection) and position from intersection point
if case == 0:
t = e2.Dot(s2) * inv_divisor
result = p0+t*ldir
#print "intersection"
if result is None:
return
else:
return result
#############################################
def Intersect_Mp(a,b,c,p0,p1,epsi) :
#intersection with Plane
dif = p0-a
rdir = p0-p1
ldir = p1-p0
e1 = b - a
e2 = c - a
##calc Matrix
v1 = c4d.Vector(rdir.x, e1.x, e2.x)
v2 = c4d.Vector(rdir.y, e1.y, e2.y)
v3 = c4d.Vector(rdir.z, e1.z, e2.z)
off = c4d.Vector(0,0,0)
D = c4d.Matrix(off,v1,v2,v3)
##invert Matrix
InD = ~D
t = dif.Dot(InD.v1)
u = dif.Dot(InD.v2)
v = dif.Dot(InD.v3)
result = p0+t*ldir
if (u+v>1 or t>1 or t<0 or v<0 or u<0) :
return
else:
return result
#############################################
def Intersect_Rc(op,p0,p1,epsi) :
ray = utils.GeRayCollider()
ray.Init(op, True)
matr = op.GetMg()
pos = p0 * matr
ldir = p1-p0
direction = ldir.GetNormalized()
raylength = ldir.GetLength()
CollisionState = ray.Intersect(pos, direction, raylength)
erg = []
count= ray.GetIntersectionCount()
print count
if count >0:
for i in xrange(count) :
result = ray.GetIntersection(i)["hitpos"]
erg.append(result)
return erg
else: return
#############################################
def main() :
#
t = time.time()
epsi = 0.000000001
erg = [] #intersections
#--------------------------------------------
#using Möller-Trumbore / Intersect_Mt()
#algo = 0
#using Matrix projection / Intersect_Mp()
algo = 1
#using c4d.RayCollider / Intersect_Rc()
#algo = 2
#--------------------------------------------
#validate object to calculate
if not op:return
##triangulate
obj = Triangulate(op)
if not obj:return
if not obj.IsInstanceOf(c4d.Opolygon) :return #R16 else GetType()
polies = obj.GetAllPolygons()
if polies == []:return
##globalize points
points = GlobalPoints(obj)
#--------------------------------------------
#validate line/ray to check
LINE=doc.SearchObject("LINE")
if not LINE: return
p0l = LINE.GetPoint(0)
p1l = LINE.GetPoint(1)
Lmatr = LINE.GetMg()
##global points line
p0 = p0l*Lmatr
p1 = p1l*Lmatr
#--------------------------------------------
#print timing for preparing
t1 = time.time() - t
print "time for preparing: "+ str(t1) + " sec"
#--------------------------------------------
#testing the algorithms
t = time.time()
if algo == 0 or algo == 1:
for i in xrange(obj.GetPolygonCount()) :
poly = obj.GetPolygon(i)
a = points[poly.a]
b = points[poly.b]
c = points[poly.c]
if algo == 0:
result = Intersect_Mt(a,b,c,p0,p1,epsi)
if result:
erg.append(result)
else :
result = Intersect_Mp(a,b,c,p0,p1,epsi)
if result:
erg.append(result)
else:
erg = Intersect_Rc(obj,p0,p1,epsi)
#--------------------------------------------
#print algorithm test timing
t1 = time.time() - t
if algo == 0:
print str(len(polies))+" polygons calculated in "+ str(t1) + " sec with Möller-Trumbore"
elif algo == 1:
print str(len(polies))+" polygons calculated in "+ str(t1) + " sec with Matrix"
else:
print str(len(polies))+" polygons calculated in "+ str(t1) + " sec with c4d.RayCollider"
#--------------------------------------------
#insert test objects
for i in xrange(len(erg)) :
sphere = c4d.BaseObject(c4d.Osphere)
sphere[c4d.PRIM_SPHERE_RAD] = 10
doc.InsertObject(sphere)
sphere.SetAbsPos(erg[i])
c4d.EventAdd()
#--------------------------------------------
if __name__=='__main__':
main()
2___
import c4d,time, math
from c4d import utils
#############################################
def Intersect_Mt(a,b,c,p0,p1,epsi) :
case = 0
result = None
ldir = p1-p0
##check if parallel
e1 = b - a
e2 = c - a
s1 = ldir % e2
divisor = s1.Dot(e1)
if 0.0-epsi <= divisor <= 0.0+epsi:
#print "case1/ parallel to triangle"
case = 1
return
#--------------------------------------------
## Compute barycentric coordinates
inv_divisor = 1.0 / divisor
d = p0 - a
b1 = d.Dot(s1) * inv_divisor
if (b1 < 0.0 or b1 > 1) :
#print "case2/ to high or to low"
case = 2
#--------------------------------------------
s2 = d % e1
b2 = ldir.Dot(s2) * inv_divisor
if (b2 < 0.0 or (b1 + b2) > 1.0) :
#print "case3/b2 < to left or to right"
case = 3
#--------------------------------------------
## Compute distance(p0-intersection) and position from intersection point
if case == 0:
t = e2.Dot(s2) * inv_divisor
result = p0+t*ldir
#print "intersection"
if result is None:
return
else:
return result
#############################################
def Intersect_Mp(a,b,c,p0,p1,epsi) :
#intersection with Plane
dif = p0-a
rdir = p0-p1
ldir = p1-p0
e1 = b - a
e2 = c - a
##calc Matrix
v1 = c4d.Vector(rdir.x, e1.x, e2.x)
v2 = c4d.Vector(rdir.y, e1.y, e2.y)
v3 = c4d.Vector(rdir.z, e1.z, e2.z)
off = c4d.Vector(0,0,0)
D = c4d.Matrix(off,v1,v2,v3)
##invert Matrix
InD = ~D
t = dif.Dot(InD.v1)
u = dif.Dot(InD.v2)
v = dif.Dot(InD.v3)
result = p0+t*ldir
if (u+v>1 or t>1 or t<0 or v<0 or u<0) :
return
else:
return result
#############################################
def Intersect_Rc(ray,p0,p1,epsi) :
ldir = p1-p0
direction = ldir.GetNormalized()
raylength = ldir.GetLength()
CollisionState = ray.Intersect(p0, direction, raylength)
erg = []
if CollisionState:
result = ray.GetIntersection(0)["hitpos"]
return result
else: return
#############################################
def main() :
#grid variables
xdim = 50
ydim = 50
#--------------------------------------------
#
t = time.time()
epsi = 0.000000001
erg = []
#--------------------------------------------
#using Möller-Trumbore / Intersect_Mt()
#algo = 0
#using Matrix projection / Intersect_Mp()
algo = 1
#using c4d.RayCollider / Intersect_Rc()
#algo = 2
#--------------------------------------------
#create a triangle
poly = c4d.BaseObject(c4d.Opolygon)
poly.ResizeObject(3,1)
a = c4d.Vector(-1500,0,100)
b = c4d.Vector(500,0,100)
c = c4d.Vector(0,1500,-1000)
poly.SetPoint(0,a)
poly.SetPoint(1,b)
poly.SetPoint(2,c)
poly. SetPolygon(0,c4d.CPolygon(0,1,2))
poly.SetAbsPos(c4d.Vector(0,1050,1500))
poly.SetName("raytestobject")
poly.Message(c4d.MSG_UPDATE)
doc.InsertObject(poly)
c4d.EventAdd()
#--------------------------------------------
#grid definition
obj = poly
Bbox = obj.GetRad()
BboxC = obj.GetMp()
BboxEx = Bbox*2
matr = obj.GetMg()
CheckXcount= BboxEx.x/xdim
CeilXcount= math.ceil(CheckXcount)
if c4d.utils.CompareFloatTolerant(CheckXcount+1.0, CeilXcount)== True:
Xcount= CheckXcount
else:
Xcount= CeilXcount
CheckYcount= BboxEx.y/ydim
CeilYcount= math.ceil(CheckYcount)
if c4d.utils.CompareFloatTolerant(CheckYcount+1.0, CeilYcount)== True:
Ycount= CheckYcount
else:
Ycount= CeilYcount
start = ( Bbox + epsi)
Zstart = start.z
Zend = Zstart - (Bbox.z +epsi)*2
Xorg= (Xcount*xdim)/2-xdim/2
Yorg= (Ycount*ydim)/2-ydim/2
direction = c4d.Vector(0,0,-1)
#--------------------------------------------
#check Grid
#'''null = c4d.BaseObject(c4d.Onull)
#'''null.SetName("rays")
#'''doc.InsertObject(null)
#'''for py in xrange(int(Ycount)) :
#'''y= float(ydim*py)
#'''for px in xrange(int(Xcount)) :
#'''x= float(xdim*px)
#'''line = c4d.BaseObject(c4d.Ospline)
#'''line.ResizeObject(2,1)
#'''line.SetPoint(0,c4d.Vector((Xorg-x),(Yorg-y),Zstart)+BboxC)
#'''line.SetPoint(1,c4d.Vector((Xorg-x),(Yorg-y),Zend)+BboxC)
#'''line.SetSegment(0,2,0)
#'''line.SetAbsPos(matr.off)
#'''line.InsertUnder(null)
#--------------------------------------------
#print timing for preparing
t1 = time.time() - t
print "time for preparing: "+ str(t1) + " sec"
#--------------------------------------------
#testing the algorithms
t = time.time()
if algo == 0 :
for py in xrange(int(Ycount)) :
y = float(ydim*py)
for px in xrange(int(Xcount)) :
x = float(xdim*px)
p0 = c4d.Vector((Xorg-x),(Yorg-y),Zstart)+BboxC
p1 = c4d.Vector((Xorg-x),(Yorg-y),Zend)+BboxC
result = Intersect_Mt(a,b,c,p0,p1,epsi)
if result:
erg.append(result * matr)
elif algo == 1 :
for py in xrange(int(Ycount)) :
y = float(ydim*py)
for px in xrange(int(Xcount)) :
x = float(xdim*px)
p0 = c4d.Vector((Xorg-x),(Yorg-y),Zstart)+BboxC
p1 = c4d.Vector((Xorg-x),(Yorg-y),Zend)+BboxC
result = Intersect_Mp(a,b,c,p0,p1,epsi)
if result:
erg.append(result * matr)
else:
ray = utils.GeRayCollider()
ray.Init(obj, True)
for py in xrange(int(Ycount)) :
y = float(ydim*py)
for px in xrange(int(Xcount)) :
x = float(xdim*px)
p0 = c4d.Vector((Xorg-x),(Yorg-y),Zstart)+BboxC
p1 = c4d.Vector((Xorg-x),(Yorg-y),Zend)+BboxC
result = Intersect_Rc(ray,p0,p1,epsi)
if result:
erg.append(result * matr)
#--------------------------------------------
#print algorithm test timing
count = Xcount*Ycount
t1 = time.time() - t
if algo == 0:
print str(count)+" rays calculated in "+ str(t1) + " sec with Möller-Trumbore"
elif algo == 1:
print str(count)+" rays calculated in "+ str(t1) + " sec with Matrix"
else:
print str(count)+" rays calculated in "+ str(t1) + " sec with c4d.RayCollider"
#--------------------------------------------
#insert test objects
#'''null_2 = c4d.BaseObject(c4d.Onull)
#'''null_2.SetName("hitpoints")
#'''doc.InsertObject(null_2)
#'''for i in xrange(len(erg)) :
#'''sphere = c4d.BaseObject(c4d.Osphere)
#'''sphere[c4d.PRIM_SPHERE_RAD] = 20
#'''sphere[c4d.PRIM_SPHERE_SUB] = 3
#'''sphere.InsertUnder(null_2)
#'''sphere.SetAbsPos(erg[i])
#'''c4d.EventAdd()
#--------------------------------------------
if __name__=='__main__':
main()
