# rays intersections octree

On 01/02/2015 at 05:06, xxxxxxxx wrote:

Hello,

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],
doc=doc)
if not obj: return                  #it was not possible to triangulate the object
return virtualop

#############################################
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)
doc.InsertObject(sphere)
sphere.SetAbsPos(erg[i])
#--------------------------------------------

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)

#--------------------------------------------

#grid definition
obj = poly
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_SUB] = 3
#'''sphere.InsertUnder(null_2)
#'''sphere.SetAbsPos(erg[i])
#--------------------------------------------

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,

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