rays intersections octree
-
On 01/02/2015 at 05:06, xxxxxxxx wrote:
Hello,
a few questions and answers referring also to this logged thread:
https://developers.maxon.net/forum/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
Martintest 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() -
On 02/02/2015 at 12:38, xxxxxxxx wrote:
Hello Martin,
I'm sorry, but we aren't allowed to disclose any implementation details of Cinema 4D.
Regarding the code examples, wouldn't it be a good idea to share those via Github or Google Code? It would probably be easier to find and access.
-
On 03/02/2015 at 01:32, xxxxxxxx wrote:
Hello Andreas,
thanks for the reply.
I see, would it be possible for you to answer question 3 and 4, as they are reffering to sdk handling?Regarding the code examples, I never thought about publishing it on another platform, cause I´m feeling just like an ambitious hobbyist and at this forum you might get a correction to learn from.
That shouldn´t mean that this forum is not professional, but in a way more collaborativ.
But yes I´ll think about it or maybe a blog like the cool kids have.
It´s a long way to scroll down
Best wishes
Martin -
On 11/03/2015 at 11:44, xxxxxxxx wrote:
Hi Martin,
I'm terribly sorry, I actually forgot there were open questions...
To be honest I have no idea, what would be the optimal way to do this in Python. I have never done such an implementation and so everything I could say would be mere guessing. Since nobody elese answered on these questions, I'm afraid you will have to test it out yourself
-
On 12/03/2015 at 13:18, xxxxxxxx wrote:
Hi Andreas,
no worries, all good.
I´ll try to wrap my head around this …Best wishes
Martin