epitaque_

Parametric Ray Octree Intersection Algorithm

May 26th, 2018
266
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. // Ray-Octree intersection based on
  2. //
  3. // An Efficient Parametric Algorithm for Octree Traversal (2000)
  4. // by J. Revelles , C. Ureña , M. Lastra
  5. //
  6. // Implementation by Jeroen Baert - www.forceflow.be - jeroen.baert@cs.kuleuven.be
  7. //
  8. // (Reformatted by Epitaque)
  9. //
  10. // Licensed under a Creative Commons Attribute-NonCommercial Sharealike 3.0 Unported license
  11. // which can be found at: http://creativecommons.org/licenses/by-nc-sa/3.0/
  12.  
  13. int first_node(double tx0, double ty0, double tz0, double txm, double tym, double tzm){
  14.     sbyte answer = 0;   // initialize to 00000000
  15.     // select the entry plane and set bits
  16.     if(tx0 > ty0){
  17.         if(tx0 > tz0){ // PLANE YZ
  18.             if(tym < tx0) answer|=2;    // set bit at position 1
  19.             if(tzm < tx0) answer|=1;    // set bit at position 0
  20.             if(txm < ty0) answer|=4;    // set bit at position 2
  21.             if(tzm < ty0) answer|=1;    // set bit at position 0
  22.             return (int) answer;
  23.         }
  24.     }
  25.     // PLANE XY
  26.     if(txm < tz0) answer|=4;    // set bit at position 2
  27.     if(tym < tz0) answer|=2;    // set bit at position 1
  28.     return (int) answer;
  29. }
  30.  
  31. int new_node(double txm, int x, double tym, int y, double tzm, int z){
  32.     if(txm < tym){
  33.         if(txm < tzm){return x;}  // YZ plane
  34.     }
  35.     else{
  36.         if(tym < tzm){return y;} // XZ plane
  37.     }
  38.     return z; // XY plane;
  39. }
  40.  
  41. void proc_subtree (Vector3 rayOrigin, Vector3 rayDirection, double tx0, double ty0, double tz0, double tx1, double ty1, double tz1, Node node, List<Node> intersectedNodes, sbyte a){
  42.     float txm, tym, tzm;
  43.     int currNode;
  44.  
  45.     if(tx1 < 0 || ty1 < 0 || tz1 < 0) return;  
  46.     if(node.IsLeaf){
  47.         intersectedNodes.Add(node);
  48.         //cout << "Reached leaf node " << node->debug_ID << endl;
  49.         return;
  50.     }
  51.     else{
  52.         //cout << "Reached node " << node->debug_ID << endl;
  53.     }  
  54.     txm = (float)(0.5*(tx0 + tx1));    
  55.     tym = (float)(0.5*(ty0 + ty1));    
  56.     tzm = (float)(0.5*(tz0 + tz1));    
  57.     currNode = first_node(tx0,ty0,tz0,txm,tym,tzm);    
  58.     do{        
  59.         switch (currNode) {        
  60.         case 0: {          
  61.             proc_subtree(rayOrigin, rayDirection, tx0,ty0,tz0,txm,tym,tzm,node.Children[a], intersectedNodes, a);
  62.             currNode = new_node(txm,4,tym,2,tzm,1);
  63.             break;}
  64.         case 1: {
  65.             proc_subtree(rayOrigin, rayDirection, tx0,ty0,tzm,txm,tym,tz1,node.Children[1^a], intersectedNodes, a);
  66.             currNode = new_node(txm,5,tym,3,tz1,8);
  67.             break;}
  68.         case 2: {
  69.             proc_subtree(rayOrigin, rayDirection, tx0,tym,tz0,txm,ty1,tzm,node.Children[2^a], intersectedNodes, a);
  70.             currNode = new_node(txm,6,ty1,8,tzm,3);
  71.             break;}
  72.         case 3: {
  73.             proc_subtree(rayOrigin, rayDirection, tx0,tym,tzm,txm,ty1,tz1,node.Children[3^a], intersectedNodes, a);
  74.             currNode = new_node(txm,7,ty1,8,tz1,8);
  75.             break;}
  76.         case 4: {
  77.             proc_subtree(rayOrigin, rayDirection, txm,ty0,tz0,tx1,tym,tzm,node.Children[4^a], intersectedNodes, a);
  78.             currNode = new_node(tx1,8,tym,6,tzm,5);
  79.             break;}
  80.         case 5: {
  81.             proc_subtree(rayOrigin, rayDirection, txm,ty0,tzm,tx1,tym,tz1,node.Children[5^a], intersectedNodes, a);
  82.             currNode = new_node(tx1,8,tym,7,tz1,8);
  83.             break;}
  84.         case 6: {
  85.             proc_subtree(rayOrigin, rayDirection, txm,tym,tz0,tx1,ty1,tzm,node.Children[6^a], intersectedNodes, a);
  86.             currNode = new_node(tx1,8,ty1,8,tzm,7);
  87.             break;}
  88.         case 7: {
  89.             proc_subtree(rayOrigin, rayDirection, txm,tym,tzm,tx1,ty1,tz1,node.Children[7^a], intersectedNodes, a);
  90.             currNode = 8;
  91.             break;}
  92.         }
  93.     } while (currNode < 8);
  94. }
  95.  
  96. public void ray_step(Node node, Vector3 rayOrigin, Vector3 rayDirection, List<Node> intersectedNodes)  {
  97.     Vector3 nodeMax = node.Min + Vector3.one * node.Size;
  98.     sbyte a = 0;
  99.  
  100.     if(rayDirection[0] < 0) {
  101.         rayOrigin[0] = - rayOrigin[0]; 
  102.         rayDirection[0] = - rayDirection[0];
  103.         a |= 4 ; //bitwise OR (latest bits are XYZ)
  104.     }
  105.     if(rayDirection[1] < 0){       
  106.         rayOrigin[1] = - rayOrigin[1];
  107.         rayDirection[1] = - rayDirection[1];
  108.         a |= 2 ;
  109.     }
  110.     if(rayDirection[2] < 0){       
  111.         rayOrigin[2] =  - rayOrigin[2];
  112.         rayDirection[2] = - rayDirection[2];
  113.         a |= 1 ;
  114.     }
  115.  
  116.     double divx = 1 / rayDirection[0]; // IEEE stability fix
  117.     double divy = 1 / rayDirection[1];
  118.     double divz = 1 / rayDirection[2];
  119.  
  120.     double tx0 = (node.Min[0] - rayOrigin[0]) * divx;
  121.     double tx1 = (nodeMax[0] - rayOrigin[0]) * divx;
  122.     double ty0 = (node.Min[1] - rayOrigin[1]) * divy;
  123.     double ty1 = (nodeMax[1] - rayOrigin[1]) * divy;
  124.     double tz0 = (node.Min[2] - rayOrigin[2]) * divz;
  125.     double tz1 = (nodeMax[2] - rayOrigin[2]) * divz;
  126.  
  127.     if(Mathd.Max(tx0,ty0,tz0) < Mathd.Min(tx1,ty1,tz1)){       
  128.         proc_subtree(rayOrigin, rayDirection, tx0,ty0,tz0,tx1,ty1,tz1,node,intersectedNodes, a);
  129.     }
  130. }
RAW Paste Data