1 // Copyright 2009-2021 Intel Corporation 2 // SPDX-License-Identifier: Apache-2.0 3 4 #pragma once 5 6 #include "../common/ray.h" 7 8 namespace embree 9 { 10 namespace isa 11 { 12 struct Cone 13 { 14 const Vec3fa p0; //!< start position of cone 15 const Vec3fa p1; //!< end position of cone 16 const float r0; //!< start radius of cone 17 const float r1; //!< end radius of cone 18 ConeCone19 __forceinline Cone(const Vec3fa& p0, const float r0, const Vec3fa& p1, const float r1) 20 : p0(p0), p1(p1), r0(r0), r1(r1) {} 21 intersectCone22 __forceinline bool intersect(const Vec3fa& org, const Vec3fa& dir, 23 BBox1f& t_o, 24 float& u0_o, Vec3fa& Ng0_o, 25 float& u1_o, Vec3fa& Ng1_o) const 26 { 27 /* calculate quadratic equation to solve */ 28 const Vec3fa v0 = p0-org; 29 const Vec3fa v1 = p1-org; 30 31 const float rl = rcp_length(v1-v0); 32 const Vec3fa P0 = v0, dP = (v1-v0)*rl; 33 const float dr = (r1-r0)*rl; 34 const Vec3fa O = -P0, dO = dir; 35 36 const float dOdO = dot(dO,dO); 37 const float OdO = dot(dO,O); 38 const float OO = dot(O,O); 39 const float dOz = dot(dP,dO); 40 const float Oz = dot(dP,O); 41 42 const float R = r0 + Oz*dr; 43 const float A = dOdO - sqr(dOz) * (1.0f+sqr(dr)); 44 const float B = 2.0f * (OdO - dOz*(Oz + R*dr)); 45 const float C = OO - (sqr(Oz) + sqr(R)); 46 47 /* we miss the cone if determinant is smaller than zero */ 48 const float D = B*B - 4.0f*A*C; 49 if (D < 0.0f) return false; 50 51 /* special case for rays that are "parallel" to the cone */ 52 const float eps = float(1<<8)*float(ulp)*max(abs(dOdO),abs(sqr(dOz))); 53 if (unlikely(abs(A) < eps)) 54 { 55 /* cylinder case */ 56 if (abs(dr) < 16.0f*float(ulp)) { 57 if (C <= 0.0f) { t_o = BBox1f(neg_inf,pos_inf); return true; } 58 else { t_o = BBox1f(pos_inf,neg_inf); return false; } 59 } 60 61 /* cone case */ 62 else 63 { 64 /* if we hit the negative cone there cannot be a hit */ 65 const float t = -C/B; 66 const float z0 = Oz+t*dOz; 67 const float z0r = r0+z0*dr; 68 if (z0r < 0.0f) return false; 69 70 /* test if we start inside or outside the cone */ 71 if (dOz*dr > 0.0f) t_o = BBox1f(t,pos_inf); 72 else t_o = BBox1f(neg_inf,t); 73 } 74 } 75 76 /* standard case for "non-parallel" rays */ 77 else 78 { 79 const float Q = sqrt(D); 80 const float rcp_2A = rcp(2.0f*A); 81 t_o.lower = (-B-Q)*rcp_2A; 82 t_o.upper = (-B+Q)*rcp_2A; 83 84 /* standard case where both hits are on same cone */ 85 if (likely(A > 0.0f)) { 86 const float z0 = Oz+t_o.lower*dOz; 87 const float z0r = r0+z0*dr; 88 if (z0r < 0.0f) return false; 89 } 90 91 /* special case where the hits are on the positive and negative cone */ 92 else 93 { 94 /* depending on the ray direction and the open direction 95 * of the cone we have a hit from inside or outside the 96 * cone */ 97 if (dOz*dr > 0) t_o.upper = pos_inf; 98 else t_o.lower = neg_inf; 99 } 100 } 101 102 /* calculates u and Ng for near hit */ 103 { 104 u0_o = (Oz+t_o.lower*dOz)*rl; 105 const Vec3fa Pr = t_o.lower*dir; 106 const Vec3fa Pl = v0 + u0_o*(v1-v0); 107 const Vec3fa R = normalize(Pr-Pl); 108 const Vec3fa U = (p1-p0)+(r1-r0)*R; 109 const Vec3fa V = cross(p1-p0,R); 110 Ng0_o = cross(V,U); 111 } 112 113 /* calculates u and Ng for far hit */ 114 { 115 u1_o = (Oz+t_o.upper*dOz)*rl; 116 const Vec3fa Pr = t_o.upper*dir; 117 const Vec3fa Pl = v0 + u1_o*(v1-v0); 118 const Vec3fa R = normalize(Pr-Pl); 119 const Vec3fa U = (p1-p0)+(r1-r0)*R; 120 const Vec3fa V = cross(p1-p0,R); 121 Ng1_o = cross(V,U); 122 } 123 return true; 124 } 125 intersectCone126 __forceinline bool intersect(const Vec3fa& org, const Vec3fa& dir, BBox1f& t_o) const 127 { 128 float u0_o; Vec3fa Ng0_o; float u1_o; Vec3fa Ng1_o; 129 return intersect(org,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o); 130 } 131 verifyCone132 static bool verify(const size_t id, const Cone& cone, const Ray& ray, bool shouldhit, const float t0, const float t1) 133 { 134 float eps = 0.001f; 135 BBox1f t; bool hit; 136 hit = cone.intersect(ray.org,ray.dir,t); 137 138 bool failed = hit != shouldhit; 139 if (shouldhit) failed |= std::isinf(t0) ? t0 != t.lower : (t0 == -1E6) ? t.lower > -1E6f : abs(t0-t.lower) > eps; 140 if (shouldhit) failed |= std::isinf(t1) ? t1 != t.upper : (t1 == +1E6) ? t.upper < +1E6f : abs(t1-t.upper) > eps; 141 if (!failed) return true; 142 embree_cout << "Cone test " << id << " failed: cone = " << cone << ", ray = " << ray << ", hit = " << hit << ", t = " << t << embree_endl; 143 return false; 144 } 145 146 /* verify cone class */ verifyCone147 static bool verify() 148 { 149 bool passed = true; 150 const Cone cone0(Vec3fa(0.0f,0.0f,0.0f),0.0f,Vec3fa(1.0f,0.0f,0.0f),1.0f); 151 passed &= verify(0,cone0,Ray(Vec3fa(-2.0f,1.0f,0.0f),Vec3fa(+1.0f,+0.0f,+0.0f),0.0f,float(inf)),true,3.0f,pos_inf); 152 passed &= verify(1,cone0,Ray(Vec3fa(+2.0f,1.0f,0.0f),Vec3fa(-1.0f,+0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,1.0f); 153 passed &= verify(2,cone0,Ray(Vec3fa(-1.0f,0.0f,2.0f),Vec3fa(+0.0f,+0.0f,-1.0f),0.0f,float(inf)),false,0.0f,0.0f); 154 passed &= verify(3,cone0,Ray(Vec3fa(+1.0f,0.0f,2.0f),Vec3fa(+0.0f,+0.0f,-1.0f),0.0f,float(inf)),true,1.0f,3.0f); 155 passed &= verify(4,cone0,Ray(Vec3fa(-1.0f,0.0f,0.0f),Vec3fa(+1.0f,+0.0f,+0.0f),0.0f,float(inf)),true,1.0f,pos_inf); 156 passed &= verify(5,cone0,Ray(Vec3fa(+1.0f,0.0f,0.0f),Vec3fa(-1.0f,+0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,1.0f); 157 passed &= verify(6,cone0,Ray(Vec3fa(+0.0f,0.0f,1.0f),Vec3fa(+0.0f,+0.0f,-1.0f),0.0f,float(inf)),true,1.0f,1.0f); 158 passed &= verify(7,cone0,Ray(Vec3fa(+0.0f,1.0f,0.0f),Vec3fa(-1.0f,-1.0f,+0.0f),0.0f,float(inf)),false,0.0f,0.0f); 159 passed &= verify(8,cone0,Ray(Vec3fa(+0.0f,1.0f,0.0f),Vec3fa(+1.0f,-1.0f,+0.0f),0.0f,float(inf)),true,0.5f,+1E6); 160 passed &= verify(9,cone0,Ray(Vec3fa(+0.0f,1.0f,0.0f),Vec3fa(-1.0f,+1.0f,+0.0f),0.0f,float(inf)),true,-1E6,-0.5f); 161 const Cone cone1(Vec3fa(0.0f,0.0f,0.0f),1.0f,Vec3fa(1.0f,0.0f,0.0f),0.0f); 162 passed &= verify(10,cone1,Ray(Vec3fa(-2.0f,1.0f,0.0f),Vec3fa(+1.0f,+0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,2.0f); 163 passed &= verify(11,cone1,Ray(Vec3fa(-1.0f,0.0f,2.0f),Vec3fa(+0.0f,+0.0f,-1.0f),0.0f,float(inf)),true,0.0f,4.0f); 164 const Cone cylinder(Vec3fa(0.0f,0.0f,0.0f),1.0f,Vec3fa(1.0f,0.0f,0.0f),1.0f); 165 passed &= verify(12,cylinder,Ray(Vec3fa(-2.0f,1.0f,0.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),true,0.0f,2.0f); 166 passed &= verify(13,cylinder,Ray(Vec3fa(+2.0f,1.0f,0.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),true,0.0f,2.0f); 167 passed &= verify(14,cylinder,Ray(Vec3fa(+2.0f,1.0f,2.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),false,0.0f,0.0f); 168 passed &= verify(15,cylinder,Ray(Vec3fa(+0.0f,0.0f,0.0f),Vec3fa( 1.0f, 0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,pos_inf); 169 passed &= verify(16,cylinder,Ray(Vec3fa(+0.0f,0.0f,0.0f),Vec3fa(-1.0f, 0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,pos_inf); 170 passed &= verify(17,cylinder,Ray(Vec3fa(+0.0f,2.0f,0.0f),Vec3fa( 1.0f, 0.0f,+0.0f),0.0f,float(inf)),false,pos_inf,neg_inf); 171 passed &= verify(18,cylinder,Ray(Vec3fa(+0.0f,2.0f,0.0f),Vec3fa(-1.0f, 0.0f,+0.0f),0.0f,float(inf)),false,pos_inf,neg_inf); 172 return passed; 173 } 174 175 /*! output operator */ 176 friend __forceinline embree_ostream operator<<(embree_ostream cout, const Cone& c) { 177 return cout << "Cone { p0 = " << c.p0 << ", r0 = " << c.r0 << ", p1 = " << c.p1 << ", r1 = " << c.r1 << "}"; 178 } 179 }; 180 181 template<int N> 182 struct ConeN 183 { 184 typedef Vec3<vfloat<N>> Vec3vfN; 185 186 const Vec3vfN p0; //!< start position of cone 187 const Vec3vfN p1; //!< end position of cone 188 const vfloat<N> r0; //!< start radius of cone 189 const vfloat<N> r1; //!< end radius of cone 190 ConeNConeN191 __forceinline ConeN(const Vec3vfN& p0, const vfloat<N>& r0, const Vec3vfN& p1, const vfloat<N>& r1) 192 : p0(p0), p1(p1), r0(r0), r1(r1) {} 193 194 __forceinline Cone operator[] (const size_t i) const 195 { 196 assert(i<N); 197 return Cone(Vec3fa(p0.x[i],p0.y[i],p0.z[i]),r0[i],Vec3fa(p1.x[i],p1.y[i],p1.z[i]),r1[i]); 198 } 199 intersectConeN200 __forceinline vbool<N> intersect(const Vec3fa& org, const Vec3fa& dir, 201 BBox<vfloat<N>>& t_o, 202 vfloat<N>& u0_o, Vec3vfN& Ng0_o, 203 vfloat<N>& u1_o, Vec3vfN& Ng1_o) const 204 { 205 /* calculate quadratic equation to solve */ 206 const Vec3vfN v0 = p0-Vec3vfN(org); 207 const Vec3vfN v1 = p1-Vec3vfN(org); 208 209 const vfloat<N> rl = rcp_length(v1-v0); 210 const Vec3vfN P0 = v0, dP = (v1-v0)*rl; 211 const vfloat<N> dr = (r1-r0)*rl; 212 const Vec3vfN O = -P0, dO = dir; 213 214 const vfloat<N> dOdO = dot(dO,dO); 215 const vfloat<N> OdO = dot(dO,O); 216 const vfloat<N> OO = dot(O,O); 217 const vfloat<N> dOz = dot(dP,dO); 218 const vfloat<N> Oz = dot(dP,O); 219 220 const vfloat<N> R = r0 + Oz*dr; 221 const vfloat<N> A = dOdO - sqr(dOz) * (vfloat<N>(1.0f)+sqr(dr)); 222 const vfloat<N> B = 2.0f * (OdO - dOz*(Oz + R*dr)); 223 const vfloat<N> C = OO - (sqr(Oz) + sqr(R)); 224 225 /* we miss the cone if determinant is smaller than zero */ 226 const vfloat<N> D = B*B - 4.0f*A*C; 227 vbool<N> valid = D >= 0.0f; 228 if (none(valid)) return valid; 229 230 /* special case for rays that are "parallel" to the cone */ 231 const vfloat<N> eps = float(1<<8)*float(ulp)*max(abs(dOdO),abs(sqr(dOz))); 232 const vbool<N> validt = valid & (abs(A) < eps); 233 const vbool<N> validf = valid & !(abs(A) < eps); 234 if (unlikely(any(validt))) 235 { 236 const vboolx validtt = validt & (abs(dr) < 16.0f*float(ulp)); 237 const vboolx validtf = validt & (abs(dr) >= 16.0f*float(ulp)); 238 239 /* cylinder case */ 240 if (unlikely(any(validtt))) 241 { 242 t_o.lower = select(validtt, select(C <= 0.0f, vfloat<N>(neg_inf), vfloat<N>(pos_inf)), t_o.lower); 243 t_o.upper = select(validtt, select(C <= 0.0f, vfloat<N>(pos_inf), vfloat<N>(neg_inf)), t_o.upper); 244 valid &= !validtt | C <= 0.0f; 245 } 246 247 /* cone case */ 248 if (any(validtf)) 249 { 250 /* if we hit the negative cone there cannot be a hit */ 251 const vfloat<N> t = -C/B; 252 const vfloat<N> z0 = Oz+t*dOz; 253 const vfloat<N> z0r = r0+z0*dr; 254 valid &= !validtf | z0r >= 0.0f; 255 256 /* test if we start inside or outside the cone */ 257 t_o.lower = select(validtf, select(dOz*dr > 0.0f, t, vfloat<N>(neg_inf)), t_o.lower); 258 t_o.upper = select(validtf, select(dOz*dr > 0.0f, vfloat<N>(pos_inf), t), t_o.upper); 259 } 260 } 261 262 /* standard case for "non-parallel" rays */ 263 if (likely(any(validf))) 264 { 265 const vfloat<N> Q = sqrt(D); 266 const vfloat<N> rcp_2A = 0.5f*rcp(A); 267 t_o.lower = select(validf, (-B-Q)*rcp_2A, t_o.lower); 268 t_o.upper = select(validf, (-B+Q)*rcp_2A, t_o.upper); 269 270 /* standard case where both hits are on same cone */ 271 const vbool<N> validft = validf & A>0.0f; 272 const vbool<N> validff = validf & !(A>0.0f); 273 if (any(validft)) { 274 const vfloat<N> z0 = Oz+t_o.lower*dOz; 275 const vfloat<N> z0r = r0+z0*dr; 276 valid &= !validft | z0r >= 0.0f; 277 } 278 279 /* special case where the hits are on the positive and negative cone */ 280 if (any(validff)) { 281 /* depending on the ray direction and the open direction 282 * of the cone we have a hit from inside or outside the 283 * cone */ 284 t_o.lower = select(validff, select(dOz*dr > 0.0f, t_o.lower, float(neg_inf)), t_o.lower); 285 t_o.upper = select(validff, select(dOz*dr > 0.0f, float(pos_inf), t_o.upper), t_o.upper); 286 } 287 } 288 289 /* calculates u and Ng for near hit */ 290 { 291 u0_o = (Oz+t_o.lower*dOz)*rl; 292 const Vec3vfN Pr = t_o.lower*Vec3vfN(dir); 293 const Vec3vfN Pl = v0 + u0_o*(v1-v0); 294 const Vec3vfN R = normalize(Pr-Pl); 295 const Vec3vfN U = (p1-p0)+(r1-r0)*R; 296 const Vec3vfN V = cross(p1-p0,R); 297 Ng0_o = cross(V,U); 298 } 299 300 /* calculates u and Ng for far hit */ 301 { 302 u1_o = (Oz+t_o.upper*dOz)*rl; 303 const Vec3vfN Pr = t_o.lower*Vec3vfN(dir); 304 const Vec3vfN Pl = v0 + u1_o*(v1-v0); 305 const Vec3vfN R = normalize(Pr-Pl); 306 const Vec3vfN U = (p1-p0)+(r1-r0)*R; 307 const Vec3vfN V = cross(p1-p0,R); 308 Ng1_o = cross(V,U); 309 } 310 return valid; 311 } 312 intersectConeN313 __forceinline vbool<N> intersect(const Vec3fa& org, const Vec3fa& dir, BBox<vfloat<N>>& t_o) const 314 { 315 vfloat<N> u0_o; Vec3vfN Ng0_o; vfloat<N> u1_o; Vec3vfN Ng1_o; 316 return intersect(org,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o); 317 } 318 }; 319 } 320 } 321 322