| 1 | // Copyright 2009-2021 Intel Corporation | 
| 2 | // SPDX-License-Identifier: Apache-2.0 | 
| 3 |  | 
| 4 | #pragma once | 
| 5 |  | 
| 6 | #include "../common/ray.h" | 
| 7 | #include "curve_intersector_precalculations.h" | 
| 8 |  | 
| 9 | namespace embree | 
| 10 | { | 
| 11 |   namespace isa | 
| 12 |   { | 
| 13 |     namespace __coneline_internal  | 
| 14 |     { | 
| 15 |       template<int M, typename Epilog, typename ray_tfar_func> | 
| 16 |         static __forceinline bool intersectCone(const vbool<M>& valid_i, | 
| 17 |                                                 const Vec3vf<M>& ray_org_in, const Vec3vf<M>& ray_dir,  | 
| 18 |                                                 const vfloat<M>& ray_tnear, const ray_tfar_func& ray_tfar, | 
| 19 |                                                 const Vec4vf<M>& v0, const Vec4vf<M>& v1, | 
| 20 |                                                 const vbool<M>& cL, const vbool<M>& cR, | 
| 21 |                                                 const Epilog& epilog) | 
| 22 |       {    | 
| 23 |         vbool<M> valid = valid_i; | 
| 24 |  | 
| 25 |         /* move ray origin closer to make calculations numerically stable */ | 
| 26 |         const vfloat<M> dOdO = sqr(ray_dir); | 
| 27 |         const vfloat<M> rcp_dOdO = rcp(dOdO); | 
| 28 |         const Vec3vf<M> center = vfloat<M>(0.5f)*(v0.xyz()+v1.xyz()); | 
| 29 |         const vfloat<M> dt = dot(center-ray_org_in,ray_dir)*rcp_dOdO; | 
| 30 |         const Vec3vf<M> ray_org = ray_org_in + dt*ray_dir; | 
| 31 |  | 
| 32 |         const Vec3vf<M> dP = v1.xyz() - v0.xyz(); | 
| 33 |         const Vec3vf<M> p0 = ray_org - v0.xyz(); | 
| 34 |         const Vec3vf<M> p1 = ray_org - v1.xyz(); | 
| 35 |          | 
| 36 |         const vfloat<M> dPdP  = sqr(dP); | 
| 37 |         const vfloat<M> dP0   = dot(p0,dP); | 
| 38 |         const vfloat<M> dP1   = dot(p1,dP);  | 
| 39 |         const vfloat<M> dOdP  = dot(ray_dir,dP); | 
| 40 |  | 
| 41 |         // intersect cone body | 
| 42 |         const vfloat<M> dr  = v0.w - v1.w; | 
| 43 |         const vfloat<M> hy  = dPdP + sqr(dr); | 
| 44 |         const vfloat<M> dO0 = dot(ray_dir,p0); | 
| 45 |         const vfloat<M> OO  = sqr(p0); | 
| 46 |         const vfloat<M> dPdP2 = sqr(dPdP); | 
| 47 |         const vfloat<M> dPdPr0 = dPdP*v0.w; | 
| 48 |          | 
| 49 |         const vfloat<M> A = dPdP2     - sqr(dOdP)*hy; | 
| 50 |         const vfloat<M> B = dPdP2*dO0 - dP0*dOdP*hy   + dPdPr0*(dr*dOdP); | 
| 51 |         const vfloat<M> C = dPdP2*OO  - sqr(dP0)*hy   + dPdPr0*(2.0f*dr*dP0 - dPdPr0); | 
| 52 |          | 
| 53 |         const vfloat<M> D = B*B - A*C; | 
| 54 |         valid &= D >= 0.0f; | 
| 55 |         if (unlikely(none(valid))) { | 
| 56 |           return false; | 
| 57 |         } | 
| 58 |  | 
| 59 |         /* standard case for "non-parallel" rays */ | 
| 60 |         const vfloat<M> Q = sqrt(D); | 
| 61 |         const vfloat<M> rcp_A = rcp(A); | 
| 62 |         /* special case for rays that are "parallel" to the cone - assume miss */ | 
| 63 |         const vbool<M> isParallel = abs(A) <= min_rcp_input; | 
| 64 |  | 
| 65 |         vfloat<M> t_cone_lower = select (isParallel, neg_inf, (-B-Q)*rcp_A); | 
| 66 |         vfloat<M> t_cone_upper = select (isParallel, pos_inf, (-B+Q)*rcp_A); | 
| 67 |         const vfloat<M> y_lower = dP0 + t_cone_lower*dOdP; | 
| 68 |         const vfloat<M> y_upper = dP0 + t_cone_upper*dOdP; | 
| 69 |         t_cone_lower = select(valid & y_lower > 0.0f & y_lower < dPdP, t_cone_lower, pos_inf); | 
| 70 |         t_cone_upper = select(valid & y_upper > 0.0f & y_upper < dPdP, t_cone_upper, neg_inf); | 
| 71 |  | 
| 72 |         const vbool<M> hitDisk0 = valid & cL; | 
| 73 |         const vbool<M> hitDisk1 = valid & cR; | 
| 74 |         const vfloat<M> rcp_dOdP = rcp(dOdP); | 
| 75 |         const vfloat<M> t_disk0 = select (hitDisk0, select (sqr(p0*dOdP-ray_dir*dP0)<(sqr(v0.w)*sqr(dOdP)), -dP0*rcp_dOdP, pos_inf), pos_inf); | 
| 76 |         const vfloat<M> t_disk1 = select (hitDisk1, select (sqr(p1*dOdP-ray_dir*dP1)<(sqr(v1.w)*sqr(dOdP)), -dP1*rcp_dOdP, pos_inf), pos_inf); | 
| 77 |         const vfloat<M> t_disk_lower = min(t_disk0, t_disk1); | 
| 78 |         const vfloat<M> t_disk_upper = max(t_disk0, t_disk1); | 
| 79 |  | 
| 80 |         const vfloat<M> t_lower = min(t_cone_lower, t_disk_lower); | 
| 81 |         const vfloat<M> t_upper = max(t_cone_upper, select(t_lower==t_disk_lower,  | 
| 82 |                                                       select(t_disk_upper==vfloat<M>(pos_inf),neg_inf,t_disk_upper),  | 
| 83 |                                                       select(t_disk_lower==vfloat<M>(pos_inf),neg_inf,t_disk_lower))); | 
| 84 |  | 
| 85 |         const vbool<M> valid_lower = valid & ray_tnear <= dt+t_lower & dt+t_lower <= ray_tfar() & t_lower != vfloat<M>(pos_inf); | 
| 86 |         const vbool<M> valid_upper = valid & ray_tnear <= dt+t_upper & dt+t_upper <= ray_tfar() & t_upper != vfloat<M>(neg_inf); | 
| 87 |  | 
| 88 |         const vbool<M> valid_first = valid_lower | valid_upper; | 
| 89 |         if (unlikely(none(valid_first))) | 
| 90 |           return false; | 
| 91 |  | 
| 92 |         const vfloat<M> t_first = select(valid_lower, t_lower, t_upper); | 
| 93 |         const vfloat<M> y_first = select(valid_lower, y_lower, y_upper); | 
| 94 |  | 
| 95 |         const vfloat<M> rcp_dPdP = rcp(dPdP); | 
| 96 |         const Vec3vf<M> dP2drr0dP = dPdP*dr*v0.w*dP; | 
| 97 |         const Vec3vf<M> dPhy = dP*hy; | 
| 98 |         const vbool<M> cone_hit_first = valid & (t_first == t_cone_lower | t_first == t_cone_upper); | 
| 99 |         const vbool<M> disk0_hit_first = valid & (t_first == t_disk0); | 
| 100 |         const Vec3vf<M> Ng_first = select(cone_hit_first, dPdP2*(p0+t_first*ray_dir)+dP2drr0dP-dPhy*y_first, select(disk0_hit_first, -dP, dP)); | 
| 101 |         const vfloat<M> u_first = select(cone_hit_first, y_first*rcp_dPdP, select(disk0_hit_first, vfloat<M>(zero), vfloat<M>(one))); | 
| 102 |  | 
| 103 |         /* invoke intersection filter for first hit */ | 
| 104 |         RoundLineIntersectorHitM<M> hit(u_first,zero,dt+t_first,Ng_first); | 
| 105 |         const bool is_hit_first = epilog(valid_first, hit); | 
| 106 |  | 
| 107 |         /* check for possible second hits before potentially accepted hit */ | 
| 108 |         const vfloat<M> t_second = t_upper; | 
| 109 |         const vfloat<M> y_second = y_upper; | 
| 110 |         const vbool<M> valid_second = valid_lower & valid_upper & (dt+t_upper <= ray_tfar()); | 
| 111 |         if (unlikely(none(valid_second))) | 
| 112 |           return is_hit_first; | 
| 113 |          | 
| 114 |         /* invoke intersection filter for second hit */ | 
| 115 |         const vbool<M> cone_hit_second = t_second == t_cone_lower | t_second == t_cone_upper; | 
| 116 |         const vbool<M> disk0_hit_second = t_second == t_disk0; | 
| 117 |         const Vec3vf<M> Ng_second = select(cone_hit_second, dPdP2*(p0+t_second*ray_dir)+dP2drr0dP-dPhy*y_second, select(disk0_hit_second, -dP, dP)); | 
| 118 |         const vfloat<M> u_second = select(cone_hit_second, y_second*rcp_dPdP, select(disk0_hit_first, vfloat<M>(zero), vfloat<M>(one))); | 
| 119 |  | 
| 120 |         hit = RoundLineIntersectorHitM<M>(u_second,zero,dt+t_second,Ng_second); | 
| 121 |         const bool is_hit_second = epilog(valid_second, hit); | 
| 122 |          | 
| 123 |         return is_hit_first | is_hit_second; | 
| 124 |       } | 
| 125 |     } | 
| 126 |  | 
| 127 |     template<int M> | 
| 128 |       struct ConeLineIntersectorHitM | 
| 129 |       { | 
| 130 |         __forceinline ConeLineIntersectorHitM() {} | 
| 131 |          | 
| 132 |         __forceinline ConeLineIntersectorHitM(const vfloat<M>& u, const vfloat<M>& v, const vfloat<M>& t, const Vec3vf<M>& Ng) | 
| 133 |           : vu(u), vv(v), vt(t), vNg(Ng) {} | 
| 134 | 	 | 
| 135 |         __forceinline void finalize() {} | 
| 136 | 	 | 
| 137 |         __forceinline Vec2f uv (const size_t i) const { return Vec2f(vu[i],vv[i]); } | 
| 138 |         __forceinline float t  (const size_t i) const { return vt[i]; } | 
| 139 |         __forceinline Vec3fa Ng(const size_t i) const { return Vec3fa(vNg.x[i],vNg.y[i],vNg.z[i]); } | 
| 140 | 	 | 
| 141 |       public: | 
| 142 |         vfloat<M> vu; | 
| 143 |         vfloat<M> vv; | 
| 144 |         vfloat<M> vt; | 
| 145 |         Vec3vf<M> vNg; | 
| 146 |       }; | 
| 147 |      | 
| 148 |     template<int M> | 
| 149 |       struct ConeCurveIntersector1 | 
| 150 |       { | 
| 151 |         typedef CurvePrecalculations1 Precalculations; | 
| 152 |          | 
| 153 |         struct ray_tfar { | 
| 154 |           Ray& ray; | 
| 155 |           __forceinline ray_tfar(Ray& ray) : ray(ray) {} | 
| 156 |           __forceinline vfloat<M> operator() () const { return ray.tfar; }; | 
| 157 |         }; | 
| 158 |  | 
| 159 |         template<typename Epilog> | 
| 160 |         static __forceinline bool intersect(const vbool<M>& valid_i, | 
| 161 |                                             Ray& ray, | 
| 162 |                                             IntersectContext* context, | 
| 163 |                                             const LineSegments* geom, | 
| 164 |                                             const Precalculations& pre, | 
| 165 |                                             const Vec4vf<M>& v0i, const Vec4vf<M>& v1i, | 
| 166 |                                             const vbool<M>& cL, const vbool<M>& cR, | 
| 167 |                                             const Epilog& epilog) | 
| 168 |         { | 
| 169 |           const Vec3vf<M> ray_org(ray.org.x, ray.org.y, ray.org.z); | 
| 170 |           const Vec3vf<M> ray_dir(ray.dir.x, ray.dir.y, ray.dir.z); | 
| 171 |           const vfloat<M> ray_tnear(ray.tnear()); | 
| 172 |           const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i); | 
| 173 |           const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i); | 
| 174 |           return  __coneline_internal::intersectCone<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray),v0,v1,cL,cR,epilog); | 
| 175 |         } | 
| 176 |       }; | 
| 177 |      | 
| 178 |     template<int M, int K> | 
| 179 |       struct ConeCurveIntersectorK | 
| 180 |       { | 
| 181 |         typedef CurvePrecalculationsK<K> Precalculations; | 
| 182 |          | 
| 183 |         struct ray_tfar { | 
| 184 |           RayK<K>& ray; | 
| 185 |           size_t k; | 
| 186 |           __forceinline ray_tfar(RayK<K>& ray, size_t k) : ray(ray), k(k) {} | 
| 187 |           __forceinline vfloat<M> operator() () const { return ray.tfar[k]; }; | 
| 188 |         }; | 
| 189 |          | 
| 190 |         template<typename Epilog> | 
| 191 |         static __forceinline bool intersect(const vbool<M>& valid_i, | 
| 192 |                                             RayK<K>& ray, size_t k, | 
| 193 |                                             IntersectContext* context, | 
| 194 |                                             const LineSegments* geom, | 
| 195 |                                             const Precalculations& pre, | 
| 196 |                                             const Vec4vf<M>& v0i, const Vec4vf<M>& v1i, | 
| 197 |                                             const vbool<M>& cL, const vbool<M>& cR, | 
| 198 |                                             const Epilog& epilog) | 
| 199 |         { | 
| 200 |           const Vec3vf<M> ray_org(ray.org.x[k], ray.org.y[k], ray.org.z[k]); | 
| 201 |           const Vec3vf<M> ray_dir(ray.dir.x[k], ray.dir.y[k], ray.dir.z[k]); | 
| 202 |           const vfloat<M> ray_tnear = ray.tnear()[k]; | 
| 203 |           const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i); | 
| 204 |           const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i); | 
| 205 |           return __coneline_internal::intersectCone<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray,k),v0,v1,cL,cR,epilog); | 
| 206 |         } | 
| 207 |       }; | 
| 208 |   } | 
| 209 | } | 
| 210 |  |