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