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://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()  
      
    


  • 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


Log in to reply