eic-opticks/CSG/csg_intersect_leaf_convexpolyhedron.h

0001 #pragma once
0002
0003 LEAF_FUNC
0004 float distance_leaf_convexpolyhedron( const float3& pos, const CSGNode* node, const float4* plan )
0005 {
0006     unsigned planeIdx = node->planeIdx() ;
0007     unsigned planeNum = node->planeNum() ;
0008     float sd = 0.f ;
0009     for(unsigned i=0 ; i < planeNum ; i++)
0010     {
0011         const float4& plane = plan[planeIdx+i];
0012         float d = plane.w ;
0013         float3 n = make_float3(plane);
0014         float sd_plane = dot(pos, n) - d ;
0015         sd = i == 0 ? sd_plane : fmaxf( sd, sd_plane );
0016     }
0017     return sd ;
0018 }
0019
0020
0021 LEAF_FUNC
0022 void intersect_leaf_convexpolyhedron( bool& valid_isect, float4& isect, const CSGNode* node, const float4* plan, const float t_min , const float3& ray_origin, const float3& ray_direction )
0023 {
0024     float t0 = -CUDART_INF_F ;
0025     float t1 =  CUDART_INF_F ;
0026
0027     float3 t0_normal = make_float3(0.f);
0028     float3 t1_normal = make_float3(0.f);
0029
0030     unsigned planeIdx = node->planeIdx() ;
0031     unsigned planeNum = node->planeNum() ;
0032
0033     for(unsigned i=0 ; i < planeNum ; i++)
0034     {
0035         const float4& plane = plan[planeIdx+i];
0036         float3 n = make_float3(plane);
0037         float dplane = plane.w ;
0038
0039          // RTCD p199,
0040          //            n.X = dplane
0041          //
0042          //             n.(o+td) = dplane
0043          //            no + t nd = dplane
0044          //                    t = (dplane - no)/nd
0045          //
0046
0047         float nd = dot(n, ray_direction); // -ve: entering, +ve exiting halfspace
0048         float no = dot(n, ray_origin ) ;  //  distance from coordinate origin to ray origin in direction of plane normal
0049         float dist = no - dplane ;        //  subtract plane distance from origin to get signed distance from plane, -ve inside
0050         float t_cand = -dist/nd ;
0051
0052         bool parallel_inside = nd == 0.f && dist < 0.f ;   // ray parallel to plane and inside halfspace
0053         bool parallel_outside = nd == 0.f && dist > 0.f ;  // ray parallel to plane and outside halfspac
0054
0055         if(parallel_inside) continue ;       // continue to next plane
0056         if(parallel_outside)
0057         {
0058             valid_isect = false ;
0059             return ;  // <-- without early exit, this still works due to infinity handling
0060         }
0061
0062         //    NB ray parallel to plane and outside halfspace
0063         //         ->  t_cand = -inf
0064         //                 nd = 0.f
0065         //                t1 -> -inf
0066
0067         if( nd < 0.f)  // entering
0068         {
0069             if(t_cand > t0)
0070             {
0071                 t0 = t_cand ;
0072                 t0_normal = n ;
0073             }
0074         }
0075         else     // exiting
0076         {
0077             if(t_cand < t1)
0078             {
0079                 t1 = t_cand ;
0080                 t1_normal = n ;
0081             }
0082         }
0083     }
0084
0085     valid_isect = t0 < t1 ;
0086     if(valid_isect)
0087     {
0088         if( t0 > t_min )
0089         {
0090             isect.x = t0_normal.x ;
0091             isect.y = t0_normal.y ;
0092             isect.z = t0_normal.z ;
0093             isect.w = t0 ;
0094         }
0095         else if( t1 > t_min )
0096         {
0097             isect.x = t1_normal.x ;
0098             isect.y = t1_normal.y ;
0099             isect.z = t1_normal.z ;
0100             isect.w = t1 ;
0101         }
0102     }
0103 }
0104
0105
0106