eic-opticks/CSG/csg_intersect_tree.h

0001 #pragma once
0002 /**
0003 csg_intersect_tree.h
0004 ======================
0005
0006 Flip off the dumpxyz code without changing DEBUG_PIDXYZ with vim replace::
0007
0008      :%s,if(dumpxyz),//if(dumpxyz),gc
0009
0010 **/
0011
0012 #if defined(__CUDACC__) || defined(__CUDABE__)
0013 #    define TREE_FUNC __forceinline__ __device__
0014 #else
0015 #    define TREE_FUNC inline
0016 #endif
0017
0018 #if defined(__CUDACC__) || defined(__CUDABE__)
0019 #include <stdio.h>
0020 #else
0021 #include <cmath>  // signbit
0022 using std::signbit ;
0023 #endif
0024
0025 #include "sflow.h"
0026 #include "csg_error.h"
0027 #include "csg_tranche.h"
0028 #include "csg_stack.h"
0029 #include "csg_postorder.h"
0030 #include "csg_pack.h"
0031 #include "csg_classify.h"
0032
0033 #include "f4_stack.h"
0034
0035 #ifdef DEBUG_RECORD
0036 #include "sc4u.h"
0037 #include "sbibit.h"
0038 #include "sbit_.h"
0039 #endif
0040
0041 /**
0042 distance_tree : gets given trees with numNode 3, 7, 15, ... where some nodes can be CSG_ZERO empties
0043 need to handle the tree without recursion, (the sdf approach of NNode.cpp in nunion::operator etc.. relies on recursion)
0044
0045
0046     0
0047   1   2
0048
0049
0050     0
0051   1    2
0052 3  4  5  6
0053
0054
0055            0
0056      1           2
0057   3    4     5      6
0058 7  8  9 10 11 12  13 14
0059
0060
0061 Should be a lot simpler than intersect_tree as CSG union/intersection/difference can be
0062 done simply by fminf fmaxf on the distances from lower levels.
0063
0064 Simply postorder traverse using a stack perhaps
0065
0066
0067
0068 **/
0069
0070
0071 TREE_FUNC
0072 float distance_tree( const float3& global_position, const CSGNode* node, const float4* plan0, const qat4* itra0 )
0073 {
0074     const int numNode = node->subNum();
0075     unsigned height = TREE_HEIGHT(numNode) ; // 1->0, 3->1, 7->2, 15->3, 31->4
0076
0077     // beginIdx, endIdx are 1-based level order tree indices, root:1 leftmost:1<<height    0:"parent" of root
0078     unsigned beginIdx = 1 << height ;  // leftmost
0079     unsigned endIdx   = 0 ;            // parent of root
0080     unsigned nodeIdx = beginIdx ;
0081
0082     F4_Stack stack = {} ;
0083     stack.curr = -1 ;
0084     float distance = 0.f ;
0085
0086     while( nodeIdx != endIdx )
0087     {
0088         unsigned depth = TREE_DEPTH(nodeIdx) ;
0089         unsigned elevation = height - depth ;
0090
0091         const CSGNode* nd = node + nodeIdx - 1 ;  // last *nd* in postorder traversal (the root) is *node*
0092         OpticksCSG_t typecode = (OpticksCSG_t)nd->typecode() ;
0093
0094         if( typecode == CSG_ZERO )
0095         {
0096             nodeIdx = POSTORDER_NEXT( nodeIdx, elevation ) ;
0097             continue ;
0098         }
0099
0100         bool node_or_leaf = typecode >= CSG_NODE ;
0101         if( node_or_leaf )
0102         {
0103             // could explicitly branch here on distance_node_contiguous or distance_leaf
0104             distance = distance_node(global_position, nd, node, plan0, itra0 ) ;
0105             stack.push(distance) ;
0106         }
0107         else
0108         {
0109             float lhs(0.f) ;
0110             float rhs(0.f) ;
0111             stack.pop2( lhs, rhs );
0112
0113             switch( typecode )
0114             {
0115                 case CSG_UNION:        distance = fminf( lhs,  rhs ) ; break ;
0116                 case CSG_INTERSECTION: distance = fmaxf( lhs,  rhs ) ; break ;
0117                 case CSG_DIFFERENCE:   distance = fmaxf( lhs, -rhs ) ; break ;
0118                 default:               distance = 0.f                ; break ;
0119             }
0120             stack.push(distance) ;
0121         }
0122         nodeIdx = POSTORDER_NEXT( nodeIdx, elevation ) ;
0123     }
0124     stack.pop(distance);
0125     return distance ;
0126 }
0127
0128 /**
0129 distance_list
0130 --------------
0131
0132 Using fminf so is implicitly a union list, could also use fmaxf and have intersection list
0133 controlled by ancillary type on the head node.
0134
0135 **/
0136
0137 TREE_FUNC
0138 float distance_list( const float3& global_position, const CSGNode* node, const float4* plan0, const qat4* itra0 )
0139 {
0140     const int numNode = node->subNum() ;
0141     float distance = RT_DEFAULT_MAX ;
0142     for(int nodeIdx=1 ; nodeIdx < numNode ; nodeIdx++)  // head node of list just carries subNum, so start from 1
0143     {
0144         const CSGNode* nd = node + nodeIdx ;
0145         float sd = distance_node(global_position, nd, node, plan0, itra0 ) ;
0146         distance = fminf( distance,  sd ) ;
0147     }
0148     return distance ;
0149 }
0150
0151
0152
0153 /**
0154 intersect_tree
0155 -----------------
0156
0157 The principal of the approach is described in a note by Andrew Kensler
0158 "Ray Tracing CSG Objects Using Single Hit Intersections",  available from
0159 the below link together with a discussion of the technique as used in the XRT raytracer.
0160
0161 * http://xrt.wikidot.com/doc:csg
0162 * http://xrt.wdfiles.com/local--files/doc%3Acsg/CSG.pdf
0163
0164
0165
0166                   +-----------------------+
0167                   |                     B |
0168        +------------------+               |
0169        | A        |       |               |
0170        |          |       |               |
0171        |          |       |               |
0172        |   0- - - 1 - - - 2 - - - - - - -[3]
0173        |          |       |               |
0174        |          |       |               |
0175        |          +-------|---------------+
0176        |                  |
0177        |                  |
0178        +------------------+
0179
0180 A ray is shot at each sub-object to find the nearest intersection, then the
0181 intersection with the sub-object is classified as one of entering, exiting or
0182 missing it. Based upon the combination of the two classifications, one of
0183 several actions is taken:
0184
0185 1. returning a hit
0186 2. returning a miss
0187 3. changing the starting point of the ray for one of the objects and then
0188    shooting this ray, classifying the intersection. In this case, the state
0189    machine enters a new loop.
0190
0191 For full details of the LUT (lookup table) and the single-intersect CSG implementation below
0192 see tests/CSGClassifyTest.cc and run the below::
0193
0194     CSGClassifyTest U
0195     CSGClassifyTest I
0196     CSGClassifyTest D
0197
0198
0199 Complete binary tree of height 4 (31 nodes) with 1-based nodeIdx in binary::
0200                                                                                                                                           depth    elevation
0201                                                                       1                                                                      0         4
0202
0203                                       10                                                            11                                       1         3
0204
0205                           100                        101                            110                           [111]                       2         2
0206
0207                  1000            1001          1010         1011             1100          1101            *1110*           1111              3         1
0208
0209              10000  10001    10010 10011    10100 10101   10110 10111     11000 11001   11010  11011   *11100* *11101*   11110   11111          4         0
0210
0211
0212
0213 CSG looping in the below implementation has been using the below complete binary tree slices(tranche)::
0214
0215     unsigned fullTree  = PACK4(  0,  0,  1 << height, 0 ) ;
0216
0217     unsigned leftIdx = 2*nodeIdx  ;    // left child of nodeIdx
0218     unsigned rightIdx = leftIdx + 1 ; // right child of nodeIdx
0219
0220     unsigned endTree   = PACK4(  0,  0,  nodeIdx,  endIdx  );
0221     unsigned leftTree  = PACK4(  0,  0,  leftIdx << (elevation-1), rightIdx << (elevation-1)) ;
0222     unsigned rightTree = PACK4(  0,  0,  rightIdx << (elevation-1), nodeIdx );
0223
0224
0225 1 << height
0226     leftmost, eg 10000
0227 0 = 1 >> 1
0228     one beyond root(1) in the sequence
0229
0230 nodeIdx
0231      node reached in the current slice of postorder sequence
0232 endIdx
0233      one beyond the last node in the current sequence (for fulltree that is 0)
0234
0235 leftTree
0236      consider example nodeIdx 111 which has elevation 2 in a height 4 tree
0237
0238      nodeIdx  :  111
0239      leftIdx  : 1110
0240      rightIdx : 1111
0241
0242      leftTree.start : leftIdx << (2-1)  : 11100
0243      leftTree.end   : rightIdx << (2-1) : 11110    one beyond the leftIdx subtree of three nodes in the postorder sequence
0244
0245 rightTree
0246     again consider nodeIdx 111
0247
0248     nodeIdx   :  111
0249     rightIdx  : 1111
0250
0251     rightTree.start : rightIdx << (2-1) : 11110     same one beyond end of leftTree is the start of the rightTree slice
0252     rightTree.end   :nodeIdx
0253
0254
0255 Now consider how different things would be with an unbalanced tree : the number of nodes traversed in a leftTree traverse
0256 of an unbalanced tree would be much more... the leftTree  would encompass the entirety of the postorder sequence
0257 up until the same end points as above.  The rightTree would not change.
0258
0259 Perhaps leftTreeOld should be replaced with leftTreeNew starting all the way from leftmost beginning of the postorder sequence::
0260
0261     unsigned leftTreeOld  = PACK4(  0,  0,  leftIdx << (elevation-1), rightIdx << (elevation-1)) ;
0262     unsigned leftTreeNew  = PACK4(  0,  0,  1 << height , rightIdx << (elevation-1)) ;
0263
0264 I suspect that when using balanced trees the leftTreeold can cause spurious intersects
0265 due to discontiguity from incomplete geometry as a result of not looping over the
0266 full prior postorder sequence.
0267
0268 Tried using leftTreeNew with a balanced tree and it still gives spurious intersects on internal boundariues,
0269 so it looks like tree balanching and the CSG algorithm as it stands are not compatible.
0270
0271 TODO: need to confirm exactly what is happening using CSGRecord, have not extinguished all hope of getting balanced
0272 to work yet, as it seems like it should be possible in principle.
0273
0274 **/
0275
0276 TREE_FUNC
0277 void intersect_tree( bool& valid_isect, float4& isect, const CSGNode* node, const float4* plan0, const qat4* itra0, const float t_min , const float3& ray_origin, const float3& ray_direction, bool dumpxyz )
0278 {
0279     const int numNode=node->subNum() ;   // SO THIS SHOULD NO LONGER EVER BE 1
0280     unsigned height = TREE_HEIGHT(numNode) ; // 1->0, 3->1, 7->2, 15->3, 31->4
0281     float propagate_epsilon = 0.0001f ;  // ?
0282     int ierr = 0 ;
0283
0284     LUT lut ;
0285     Tranche tr ;
0286     tr.curr = -1 ;
0287
0288     unsigned fullTree = PACK4(0,0, 1 << height, 0 ) ;  // leftmost: 1<<height,  root:1>>1 = 0 ("parent" of root)
0289
0290 #ifdef DEBUG
0291     printf("//intersect_tree  numNode(subNum) %d height %d fullTree(hex) %x \n", numNode, height, fullTree );
0292     assert( numNode > 0 );
0293 #endif
0294
0295 #ifdef DEBUG_PIDXYZ
0296     //if(dumpxyz) printf("//intersect_tree  numNode(subNum) %d height %d fullTree(hex) %x \n", numNode, height, fullTree );
0297 #endif
0298
0299
0300
0301
0302     tranche_push( tr, fullTree, t_min );
0303
0304     CSG_Stack csg ;
0305     csg.curr = -1 ;
0306 #ifdef DEBUG_RECORD
0307     int tloop = -1 ;
0308 #endif
0309
0310     while (tr.curr > -1)
0311     {
0312 #ifdef DEBUG_RECORD
0313         tloop++ ;
0314 #endif
0315         unsigned slice ;
0316         float tmin ;    // NB different from t_min when looping
0317
0318         ierr = tranche_pop(tr, slice, tmin );
0319         if(ierr) break ;
0320
0321         // nodeIdx, beginIdx, endIdx are 1-based level order tree indices, root:1 leftmost:1<<height
0322         unsigned beginIdx = UNPACK4_2(slice);
0323         unsigned endIdx   = UNPACK4_3(slice);
0324
0325 #ifdef DEBUG_RECORD
0326         if(CSGRecord::ENABLED)
0327             printf("// tranche_pop : tloop %d tmin %10.4f beginIdx %d endIdx %d  tr.curr %d csg.curr %d   \n", tloop, tmin, beginIdx, endIdx, tr.curr, csg.curr );
0328         assert( ierr == 0 );
0329 #endif
0330
0331         unsigned nodeIdx = beginIdx ;
0332         while( nodeIdx != endIdx )
0333         {
0334             unsigned depth = TREE_DEPTH(nodeIdx) ;
0335             unsigned elevation = height - depth ;
0336
0337             const CSGNode* nd = node + nodeIdx - 1 ;   // nodeIdx is 1-based level order index
0338             OpticksCSG_t typecode = (OpticksCSG_t)nd->typecode() ;
0339 #ifdef DEBUG
0340             printf("//intersect_tree  nodeIdx %d CSG::Name %10s depth %d elevation %d \n", nodeIdx, CSG::Name(typecode), depth, elevation );
0341 #endif
0342 #ifdef DEBUG_PIDXYZ
0343             //if(dumpxyz) printf("//intersect_tree nodeIdx %d typecode %d depth %d elevation %d \n", nodeIdx, typecode, depth, elevation );
0344 #endif
0345
0346             if( typecode == CSG_ZERO )
0347             {
0348                 nodeIdx = POSTORDER_NEXT( nodeIdx, elevation ) ;
0349                 continue ;
0350             }
0351             bool node_or_leaf = typecode >= CSG_NODE ;
0352 #ifdef DEBUG_PIDXYZ
0353             //if(dumpxyz) printf("//intersect_tree  nodeIdx %d node_or_leaf %d \n", nodeIdx, node_or_leaf );
0354 #endif
0355
0356             if(node_or_leaf)
0357             {
0358                 float4 nd_isect = make_float4(0.f, 0.f, 0.f, 0.f) ;
0359
0360                 bool nd_valid_isect(false);
0361                 intersect_node( nd_valid_isect, nd_isect, nd, node, plan0, itra0, tmin, ray_origin, ray_direction, dumpxyz );
0362                 // funny that hit or miss does not seem to be used here : presumably relying on the t value
0363
0364 #ifdef DEBUG_PIDXYZ
0365                 //if(dumpxyz) printf("//intersect_tree nd_valid_isect:%d \n", nd_valid_isect ) ;
0366 #endif
0367                 //if( !nd_valid_isect ) break ;   // BEHAVIOUR BREAKING : "SPHERE CLIPS" NOT VALID JUST MEANS A MISS
0368
0369                 nd_isect.w = copysignf( nd_isect.w, nodeIdx % 2 == 0 ? -1.f : 1.f );  // hijack t signbit, to record the side, LHS -ve
0370
0371 #ifdef DEBUG_PIDXYZ
0372                 //if(dumpxyz) printf("//intersect_tree  nodeIdx %d node_or_leaf %d nd_isect (%10.4f %10.4f %10.4f %10.4f) nd_valid_isect %d \n", nodeIdx, node_or_leaf, nd_isect.x, nd_isect.y, nd_isect.z, nd_isect.w, nd_valid_isect );
0373 #endif
0374                 ierr = csg_push(csg, nd_isect, nodeIdx );
0375
0376
0377 #ifdef DEBUG_RECORD
0378                 assert( ierr == 0 );
0379 #endif
0380                 if(ierr) break ;
0381             }
0382             else
0383             {
0384                 if(csg.curr < 1)  // curr 1 : 2 items
0385                 {
0386 #ifdef DEBUG
0387                     printf("//intersect_tree  ERROR_POP_EMPTY nodeIdx %4d typecode %d csg.curr %d \n", nodeIdx, typecode, csg.curr );
0388 #endif
0389                     ierr |= ERROR_POP_EMPTY ;
0390                     break ;
0391                 }
0392
0393                 // operator node : peek at the top of the stack
0394
0395                 bool firstLeft = signbit(csg.data[csg.curr].w) ;
0396                 bool secondLeft = signbit(csg.data[csg.curr-1].w) ;
0397
0398 #ifdef DEBUG_PIDXYZ
0399                //if(dumpxyz) printf("//intersect_tree  nodeIdx %d firstLeft %d secondLeft %d \n", nodeIdx, firstLeft, secondLeft );
0400 #endif
0401
0402
0403                 if(!(firstLeft ^ secondLeft))
0404                 {
0405 #ifdef DEBUG
0406                     printf("//intersect_tree ERROR_XOR_SIDE nodeIdx %4d typecode %d tl %10.3f tr %10.3f sl %d sr %d \n",
0407                              nodeIdx, typecode, csg.data[csg.curr].w, csg.data[csg.curr-1].w, firstLeft, secondLeft );
0408 #endif
0409                     ierr |= ERROR_XOR_SIDE ;
0410                     break ;
0411                 }
0412                 int left  = firstLeft ? csg.curr   : csg.curr-1 ;
0413                 int right = firstLeft ? csg.curr-1 : csg.curr   ;
0414
0415                 IntersectionState_t l_state = CSG_CLASSIFY( csg.data[left],  ray_direction, tmin );
0416                 IntersectionState_t r_state = CSG_CLASSIFY( csg.data[right], ray_direction, tmin );
0417
0418                 float t_left  = fabsf( csg.data[left].w );
0419                 float t_right = fabsf( csg.data[right].w );
0420
0421                 bool leftIsCloser = t_left <= t_right ;
0422
0423 #ifdef DEBUG_PIDXYZ
0424                 //  state :   0:Enter 1:Exit 2:Miss
0425                 //if(dumpxyz) printf("//intersect_tree  nodeIdx %d left %d right %d l_state %d r_state %d t_left %10.3f t_right %10.3f leftIsCloser %d \n", nodeIdx, left, right, l_state, r_state, t_left, t_right, leftIsCloser );
0426 #endif
0427
0428 #ifdef DEBUG_COS
0429                 {
0430                     // dot products of ray_direction and normals for the two isect
0431                     float l_cos = csg.data[left].x*ray_direction.x + csg.data[left].y*ray_direction.y + csg.data[left].z*ray_direction.z ;
0432                     float r_cos = csg.data[right].x*ray_direction.x + csg.data[right].y*ray_direction.y + csg.data[right].z*ray_direction.z ;
0433
0434                     printf("\n//intersect_tree nodeIdx %3d t_left %10.4f t_right %10.4f leftIsCloser %d  l_state %5s r_state %5s l_cos*1e6f %10.4f r_cos*1e6f %10.4f \n",
0435                              nodeIdx, t_left, t_right, leftIsCloser,
0436                              IntersectionState::Name(l_state),
0437                              IntersectionState::Name(r_state),
0438                              l_cos*1e6f,
0439                              r_cos*1e6f
0440                     );
0441
0442                 }
0443 #endif
0444                 // it is impossible to Miss a complemented (signaled by -0.f) solid as it is unbounded
0445                 // hence the below artificially changes leftIsCloser to reflect the unboundedness
0446                 // and sets the corresponding states to Exit
0447                 // see opticks/notes/issues/csg_complement.rst
0448                 // these settings are only valid (and only needed) for misses
0449
0450                 // Q: where are complements set into the intersects ?
0451                 // A: at the tail of intersect_leaf the normal is flipped for complemented solids
0452                 //    even (actually especially) for MISS
0453                 //
0454
0455                 bool l_complement = signbit(csg.data[left].x) ;
0456                 bool r_complement = signbit(csg.data[right].x) ;
0457
0458                 // unbounded values are only valid for CSG_THETACUT State_Miss
0459                 bool l_unbounded = signbit(csg.data[left].y) ;
0460                 bool r_unbounded = signbit(csg.data[right].y) ;
0461
0462                 bool l_promote_miss = l_state == State_Miss && ( l_complement || l_unbounded ) ;
0463                 bool r_promote_miss = r_state == State_Miss && ( r_complement || r_unbounded ) ;
0464
0465 #ifdef DEBUG_PIDXYZ
0466                 //  state :   0:Enter 1:Exit 2:Miss
0467                 //if(dumpxyz) printf("//intersect_tree  nodeIdx %d l/r_complement %d/%d l/r_unbounded %d/%d l/r_promote_miss %d/%d \n", nodeIdx, l_complement, r_complement, l_unbounded, r_unbounded, l_promote_miss, r_promote_miss );
0468 #endif
0469
0470                 if(r_promote_miss)
0471                 {
0472 #ifdef DEBUG_RECORD
0473                     if(CSGRecord::ENABLED)
0474                     {
0475                         printf("// %3d : r_promote_miss setting leftIsCloser %d to true and r_state %5s to Exit \n",
0476                                 nodeIdx, leftIsCloser, IntersectionState::Name(l_state)  );
0477                     }
0478 #endif
0479                     r_state = State_Exit ;
0480                     leftIsCloser = true ;
0481                }
0482
0483                 if(l_promote_miss)
0484                 {
0485 #ifdef DEBUG_RECORD
0486                     if(CSGRecord::ENABLED)
0487                     {
0488                         printf("// %3d : l_promote_miss setting leftIsCloser %d to false and l_state %5s to Exit \n",
0489                                 nodeIdx, leftIsCloser, IntersectionState::Name(r_state)  );
0490                     }
0491 #endif
0492                     l_state = State_Exit ;
0493                     leftIsCloser = false ;
0494                 }
0495
0496                 int ctrl = lut.lookup( typecode , l_state, r_state, leftIsCloser ) ;
0497
0498 #ifdef DEBUG
0499                 printf("// %3d : stack peeking : left %d right %d (stackIdx)  %15s  l:%5s %10.4f    r:%5s %10.4f     leftIsCloser %d -> %s \n",
0500                            nodeIdx,left,right,
0501                            CSG::Name(typecode),
0502                            IntersectionState::Name(l_state), t_left,
0503                            IntersectionState::Name(r_state), t_right,
0504                            leftIsCloser,
0505                            CTRL::Name(ctrl)  );
0506 #endif
0507
0508
0509 #ifdef DEBUG_RECORD
0510                 if(CSGRecord::ENABLED)
0511                 {
0512                     printf("// %3d CSG decision : left %d right %d (stackIdx)  %15s  l:%5s %10.4f    r:%5s %10.4f     leftIsCloser %d -> %s \n",
0513                            nodeIdx,left,right,
0514                            CSG::Name(typecode),
0515                            IntersectionState::Name(l_state), t_left,
0516                            IntersectionState::Name(r_state), t_right,
0517                            leftIsCloser,
0518                            CTRL::Name(ctrl)  );
0519
0520                     quad6 rec ;
0521                     rec.zero();
0522
0523                     rec.q0.f = csg.data[left] ;
0524                     rec.q1.f = csg.data[right] ;
0525
0526                     U4U uu ;   // union u4:4*uchar u:unsigned
0527                     uu.u4.x = sbibit_PACK4( typecode, l_state, r_state, leftIsCloser ) ;  // 4*2bit into 8bit uchar
0528                     uu.u4.y = sbit_rPACK8( l_promote_miss, l_complement, l_unbounded, false, r_promote_miss, r_complement, r_unbounded, false );  // 8 bits into uchar
0529                     uu.u4.z = tloop ;
0530                     uu.u4.w = nodeIdx ;
0531
0532                     rec.q2.u.x = uu.u ;
0533                     rec.q2.u.y = ctrl ;
0534                     rec.q2.u.z = 0u ;
0535                     rec.q2.u.w = 0u ;
0536
0537                     // HMM: tmin arriving here not the advanced one when just looping a leaf
0538                     rec.q3.f.x = tmin ;  // specific to this tranche, can be advanced compared to t_min
0539                     rec.q3.f.y = t_min ; // overall
0540                     rec.q3.f.z = 0.f ;   // set to tminAdvanced when looping below
0541                     rec.q3.f.w = 0.f ;
0542
0543                     // q4 is set below with result
0544                     // q5 currently unused
0545
0546                     CSGRecord::record.push_back(rec);
0547                 }
0548 #endif
0549
0550
0551                 unsigned act = UNDEFINED ;
0552
0553                 if(ctrl < CTRL_LOOP_A) // non-looping : CTRL_RETURN_MISS/CTRL_RETURN_A/CTRL_RETURN_B/CTRL_RETURN_FLIP_B "returning" with a push
0554                 {
0555                     float4 result = ctrl == CTRL_RETURN_MISS ?  make_float4(0.f, 0.f, 0.f, 0.f ) : csg.data[ctrl == CTRL_RETURN_A ? left : right] ;
0556                     if(ctrl == CTRL_RETURN_FLIP_B)
0557                     {
0558                         result.x = -result.x ;
0559                         result.y = -result.y ;
0560                         result.z = -result.z ;
0561                     }
0562                     result.w = copysignf( result.w , nodeIdx % 2 == 0 ? -1.f : 1.f );
0563                     // record left/right in sign of t
0564
0565                     ierr = csg_pop0(csg); if(ierr) break ;
0566                     ierr = csg_pop0(csg); if(ierr) break ;
0567                     ierr = csg_push(csg, result, nodeIdx );  if(ierr) break ;
0568
0569                     act = CONTINUE ;
0570
0571 #ifdef DEBUG_RECORD
0572                     if(CSGRecord::ENABLED)
0573                     {
0574                         quad6& rec = CSGRecord::record.back();
0575                         rec.q4.f = result ;
0576                     }
0577 #endif
0578
0579                 }
0580                 else   //   CTRL_LOOP_A/CTRL_LOOP_B
0581                 {
0582                     int loopside  = ctrl == CTRL_LOOP_A ? left : right ;
0583                     int otherside = ctrl == CTRL_LOOP_A ? right : left ;
0584
0585                     unsigned leftIdx = 2*nodeIdx ;
0586                     unsigned rightIdx = leftIdx + 1;
0587                     unsigned otherIdx = ctrl == CTRL_LOOP_A ? rightIdx : leftIdx ;
0588
0589                     float tminAdvanced = fabsf(csg.data[loopside].w) + propagate_epsilon ;
0590                     float4 other = csg.data[otherside] ;  // need tmp as pop about to invalidate indices
0591
0592                     ierr = csg_pop0(csg);                   if(ierr) break ;
0593                     ierr = csg_pop0(csg);                   if(ierr) break ;
0594                     ierr = csg_push(csg, other, otherIdx ); if(ierr) break ;
0595
0596                     // looping is effectively backtracking, pop both and put otherside back
0597
0598                     unsigned endTree       = PACK4(  0,  0,  nodeIdx,  endIdx  );
0599                     unsigned leftTree      = PACK4(  0,  0,  leftIdx << (elevation-1), rightIdx << (elevation-1)) ;
0600                     //unsigned leftTreeNew   = PACK4(  0,  0,  1 << height             , rightIdx << (elevation-1)) ;
0601                     unsigned rightTree     = PACK4(  0,  0,  rightIdx << (elevation-1), nodeIdx );
0602
0603                     unsigned loopTree  = ctrl == CTRL_LOOP_A ? leftTree : rightTree  ;
0604
0605 #ifdef DEBUG_RECORD
0606                     if(CSGRecord::ENABLED)
0607                     {
0608                         printf("// %3d : looping one side tminAdvanced %10.4f with eps %10.4f \n", nodeIdx, tminAdvanced, propagate_epsilon );
0609                         quad6& rec = CSGRecord::record.back();
0610                         rec.q3.f.z = tminAdvanced ;
0611                     }
0612 #endif
0613
0614 #ifdef DEBUG
0615                     printf("//intersect_tree nodeIdx %2d height %2d depth %2d elevation %2d endTree %8x leftTree %8x rightTree %8x \n",
0616                               nodeIdx,
0617                               height,
0618                               depth,
0619                               elevation,
0620                               endTree,
0621                               leftTree,
0622                               rightTree);
0623 #endif
0624
0625                    // push the tranche from here to endTree before pushing the backtracking tranche so known how to proceed after backtracking done
0626                    // (hmm: using tmin onwards to endTree looks a bit funny, maybe it should be advanced?)
0627
0628                     ierr = tranche_push( tr, endTree,  tmin );         if(ierr) break ;
0629                     ierr = tranche_push( tr, loopTree, tminAdvanced ); if(ierr) break ;
0630
0631                     act = BREAK  ;
0632
0633 #ifdef DEBUG_RECORD
0634                     if(CSGRecord::ENABLED)
0635                     printf("// %3d : looping :  act BREAK \n", nodeIdx );
0636 #endif
0637
0638                 }                      // "return" or "recursive call"
0639
0640
0641                 if(act == BREAK)
0642                 {
0643 #ifdef DEBUG_RECORD
0644                      if(CSGRecord::ENABLED)
0645                      printf("// %3d : break for backtracking \n", nodeIdx );
0646 #endif
0647                      break ;
0648                 }
0649             }                          // "primitive" or "operation"
0650             nodeIdx = POSTORDER_NEXT( nodeIdx, elevation ) ;
0651         }                     // node traversal
0652         if(ierr) break ;
0653     }                        // subtree tranches
0654
0655     ierr |= (( csg.curr !=  0)  ? ERROR_END_EMPTY : 0)  ;
0656
0657 #ifdef DEBUG_RECORD
0658     if(CSGRecord::ENABLED)
0659     printf("// intersect_tree ierr %d csg.curr %d \n", ierr, csg.curr );
0660 #endif
0661     if(csg.curr == 0)
0662     {
0663         const float4& ret = csg.data[0] ;
0664         isect.x = ret.x ;
0665         isect.y = ret.y ;
0666         isect.z = ret.z ;
0667         isect.w = ret.w ;
0668     }
0669
0670     valid_isect = isect.w > 0.f ;  // ?
0671     return ;
0672 }
0673
0674 /**
0675 intersect_prim
0676 ----------------
0677
0678 Canonically invoked from CSGOptiX/CSGOptiX7.cu:__intersection__is
0679 with object frame ray_origin and ray_direction
0680
0681 **/
0682
0683 TREE_FUNC
0684 bool intersect_prim( float4& isect, const CSGNode* node, const float4* plan, const qat4* itra, const float t_min , const float3& ray_origin, const float3& ray_direction, bool dumpxyz )
0685 {
0686     const unsigned typecode = node->typecode() ;
0687 #ifdef DEBUG
0688     printf("//intersect_prim typecode %u1 name %s \n", typecode, CSG::Name(typecode) );
0689 #endif
0690 #ifdef DEBUG_PIDXYZ
0691     //if(dumpxyz) printf("//intersect_prim typecode %u \n", typecode  );
0692 #endif
0693
0694
0695     bool valid_isect = false ;
0696     if( typecode >= CSG_LEAF )
0697     {
0698         intersect_leaf( valid_isect,  isect, node, plan, itra, t_min, ray_origin, ray_direction, dumpxyz ) ;
0699     }
0700     else if( typecode < CSG_NODE )
0701     {
0702         intersect_tree( valid_isect,  isect, node, plan, itra, t_min, ray_origin, ray_direction, dumpxyz ) ;
0703     }
0704 #ifdef WITH_CONTIGUOUS
0705     else if( typecode == CSG_CONTIGUOUS )
0706     {
0707         intersect_node_contiguous( valid_isect, isect, node, node, plan, itra, t_min, ray_origin, ray_direction, dumpxyz ) ;
0708     }
0709 #endif
0710     else if( typecode == CSG_DISCONTIGUOUS )
0711     {
0712         intersect_node_discontiguous( valid_isect, isect, node, node, plan, itra, t_min, ray_origin, ray_direction, dumpxyz ) ;
0713     }
0714     else if( typecode == CSG_OVERLAP )
0715     {
0716         intersect_node_overlap(       valid_isect, isect, node, node, plan, itra, t_min, ray_origin, ray_direction, dumpxyz ) ;
0717     }
0718     return valid_isect ;
0719 }
0720
0721
0722 TREE_FUNC
0723 float distance_prim( const float3& global_position, const CSGNode* node, const float4* plan, const qat4* itra )
0724 {
0725     //const int numNode = node->subNum();
0726     const unsigned typecode = node->typecode() ;
0727     float distance = -1.f ;
0728     if( typecode >= CSG_LEAF )       // leaf nodes
0729     {
0730         distance = distance_leaf(global_position, node, plan, itra );
0731     }
0732     else if( typecode < CSG_NODE )   // boolean tree nodes
0733     {
0734         distance = distance_tree(global_position, node, plan, itra );
0735     }
0736     else                             // list nodes
0737     {
0738         distance = distance_node_list( typecode, global_position, node, node, plan, itra ) ;
0739     }
0740     return distance ;
0741 }
0742
0743