eic-opticks/CSG/csg_intersect_leaf_sphere.h

0001 #pragma once
0002
0003 LEAF_FUNC
0004 float distance_leaf_sphere(const float3& pos, const quad& q0 )
0005 {
0006     float3 center = make_float3(q0.f);
0007     float radius = q0.f.w;
0008     float3 p = pos - center;
0009     float sd = length(p) - radius ;
0010     return sd ;
0011 }
0012
0013
0014 LEAF_FUNC
0015 void intersect_leaf_sphere( bool& valid_isect, float4& isect, const quad& q0, const float& t_min, const float3& ray_origin, const float3& ray_direction )
0016 {
0017     float3 center = make_float3(q0.f);
0018     float radius = q0.f.w;
0019
0020     float3 O = ray_origin - center;
0021     float3 D = ray_direction;
0022
0023     float b = dot(O, D);
0024     float c = dot(O, O)-radius*radius;
0025     float d = dot(D, D);
0026
0027 #ifdef CATASTROPHIC_SUBTRACTION_ROOTS
0028     float disc = b*b-d*c;   // when d*c small,  sdisc ~ b => catastrophic precision loss in root2 = (-b + sdisc)/d
0029     float sdisc = disc > 0.f ? sqrtf(disc) : 0.f ;   // repeated root for sdisc 0.f
0030     float root1 = (-b - sdisc)/d ;
0031     float root2 = (-b + sdisc)/d ;  // root2 > root1 always
0032 #else
0033     float root1, root2, disc, sdisc ;
0034     robust_quadratic_roots(root1, root2, disc, sdisc, d, b, c ) ; //  Solving:  d t^2 + 2 b t +  c = 0    root2 > root1
0035 #endif
0036
0037     float t_cand = sdisc > 0.f ? ( root1 > t_min ? root1 : root2 ) : t_min ;
0038
0039     valid_isect = t_cand > t_min ;
0040     if(valid_isect)
0041     {
0042         isect.x = (O.x + t_cand*D.x)/radius ;   // normalized by construction
0043         isect.y = (O.y + t_cand*D.y)/radius ;
0044         isect.z = (O.z + t_cand*D.z)/radius ;
0045         isect.w = t_cand ;
0046     }
0047
0048 #ifdef DEBUG
0049     //printf("//intersect_leaf_sphere valid_isect %d  isect ( %10.4f %10.4f %10.4f %10.4f)  \n", valid_isect, isect.x, isect.y, isect.z, isect.w );
0050     printf("//intersect_leaf_sphere valid %d radius %10.4f center (%10.4f, %10.4f, %10.4f) ray_ori (%10.4f, %10.4f, %10.4f)  \n",
0051        valid_isect,  radius, center.x, center.y, center.z, ray_origin.x, ray_origin.y, ray_origin.z  );
0052 #endif
0053
0054
0055 }
0056
0057