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 Cylinder
13     {
14       const Vec3fa p0;  //!< start location
15       const Vec3fa p1;  //!< end position
16       const float rr;   //!< squared radius of cylinder
17 
CylinderCylinder18       __forceinline Cylinder(const Vec3fa& p0, const Vec3fa& p1, const float r)
19         : p0(p0), p1(p1), rr(sqr(r)) {}
20 
CylinderCylinder21       __forceinline Cylinder(const Vec3fa& p0, const Vec3fa& p1, const float rr, bool)
22         : p0(p0), p1(p1), rr(rr) {}
23 
intersectCylinder24       __forceinline bool intersect(const Vec3fa& org,
25                                    const Vec3fa& dir,
26                                    BBox1f& t_o,
27                                    float& u0_o, Vec3fa& Ng0_o,
28                                    float& u1_o, Vec3fa& Ng1_o) const
29       {
30         /* calculate quadratic equation to solve */
31         const float rl = rcp_length(p1-p0);
32         const Vec3fa P0 = p0, dP = (p1-p0)*rl;
33         const Vec3fa O = org-P0, dO = dir;
34 
35         const float dOdO = dot(dO,dO);
36         const float OdO = dot(dO,O);
37         const float OO = dot(O,O);
38         const float dOz = dot(dP,dO);
39         const float Oz = dot(dP,O);
40 
41         const float A = dOdO - sqr(dOz);
42         const float B = 2.0f * (OdO - dOz*Oz);
43         const float C = OO - sqr(Oz) - rr;
44 
45         /* we miss the cylinder if determinant is smaller than zero */
46         const float D = B*B - 4.0f*A*C;
47         if (D < 0.0f) {
48           t_o = BBox1f(pos_inf,neg_inf);
49           return false;
50         }
51 
52         /* special case for rays that are parallel to the cylinder */
53         const float eps = 16.0f*float(ulp)*max(abs(dOdO),abs(sqr(dOz)));
54         if (abs(A) < eps)
55         {
56           if (C <= 0.0f) {
57             t_o = BBox1f(neg_inf,pos_inf);
58             return true;
59           } else {
60             t_o = BBox1f(pos_inf,neg_inf);
61             return false;
62           }
63         }
64 
65         /* standard case for rays that are not parallel to the cylinder */
66         const float Q = sqrt(D);
67         const float rcp_2A = rcp(2.0f*A);
68         const float t0 = (-B-Q)*rcp_2A;
69         const float t1 = (-B+Q)*rcp_2A;
70 
71         /* calculates u and Ng for near hit */
72         {
73           u0_o = madd(t0,dOz,Oz)*rl;
74           const Vec3fa Pr = t0*dir;
75           const Vec3fa Pl = madd(u0_o,p1-p0,p0);
76           Ng0_o = Pr-Pl;
77         }
78 
79         /* calculates u and Ng for far hit */
80         {
81           u1_o = madd(t1,dOz,Oz)*rl;
82           const Vec3fa Pr = t1*dir;
83           const Vec3fa Pl = madd(u1_o,p1-p0,p0);
84           Ng1_o = Pr-Pl;
85         }
86 
87         t_o.lower = t0;
88         t_o.upper = t1;
89         return true;
90       }
91 
intersectCylinder92       __forceinline bool intersect(const Vec3fa& org_i, const Vec3fa& dir, BBox1f& t_o) const
93       {
94         float u0_o; Vec3fa Ng0_o;
95         float u1_o; Vec3fa Ng1_o;
96         return intersect(org_i,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o);
97       }
98 
verifyCylinder99       static bool verify(const size_t id, const Cylinder& cylinder, const RayHit& ray, bool shouldhit, const float t0, const float t1)
100       {
101         float eps = 0.001f;
102         BBox1f t; bool hit;
103         hit = cylinder.intersect(ray.org,ray.dir,t);
104 
105         bool failed = hit != shouldhit;
106         if (shouldhit) failed |= std::isinf(t0) ? t0 != t.lower : abs(t0-t.lower) > eps;
107         if (shouldhit) failed |= std::isinf(t1) ? t1 != t.upper : abs(t1-t.upper) > eps;
108         if (!failed) return true;
109         embree_cout << "Cylinder test " << id << " failed: cylinder = " << cylinder << ", ray = " << ray << ", hit = " << hit << ", t = " << t << embree_endl;
110         return false;
111       }
112 
113       /* verify cylinder class */
verifyCylinder114       static bool verify()
115       {
116         bool passed = true;
117         const Cylinder cylinder(Vec3fa(0.0f,0.0f,0.0f),Vec3fa(1.0f,0.0f,0.0f),1.0f);
118         passed &= verify(0,cylinder,RayHit(Vec3fa(-2.0f,1.0f,0.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),true,0.0f,2.0f);
119         passed &= verify(1,cylinder,RayHit(Vec3fa(+2.0f,1.0f,0.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),true,0.0f,2.0f);
120         passed &= verify(2,cylinder,RayHit(Vec3fa(+2.0f,1.0f,2.0f),Vec3fa( 0.0f,-1.0f,+0.0f),0.0f,float(inf)),false,0.0f,0.0f);
121         passed &= verify(3,cylinder,RayHit(Vec3fa(+0.0f,0.0f,0.0f),Vec3fa( 1.0f, 0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,pos_inf);
122         passed &= verify(4,cylinder,RayHit(Vec3fa(+0.0f,0.0f,0.0f),Vec3fa(-1.0f, 0.0f,+0.0f),0.0f,float(inf)),true,neg_inf,pos_inf);
123         passed &= verify(5,cylinder,RayHit(Vec3fa(+0.0f,2.0f,0.0f),Vec3fa( 1.0f, 0.0f,+0.0f),0.0f,float(inf)),false,pos_inf,neg_inf);
124         passed &= verify(6,cylinder,RayHit(Vec3fa(+0.0f,2.0f,0.0f),Vec3fa(-1.0f, 0.0f,+0.0f),0.0f,float(inf)),false,pos_inf,neg_inf);
125         return passed;
126       }
127 
128       /*! output operator */
129       friend __forceinline embree_ostream operator<<(embree_ostream cout, const Cylinder& c) {
130         return cout << "Cylinder { p0 = " << c.p0 << ", p1 = " << c.p1 << ", r = " << sqrtf(c.rr) << "}";
131       }
132     };
133 
134     template<int N>
135       struct CylinderN
136     {
137       const Vec3vf<N> p0;     //!< start location
138       const Vec3vf<N> p1;     //!< end position
139       const vfloat<N> rr;   //!< squared radius of cylinder
140 
CylinderNCylinderN141       __forceinline CylinderN(const Vec3vf<N>& p0, const Vec3vf<N>& p1, const vfloat<N>& r)
142         : p0(p0), p1(p1), rr(sqr(r)) {}
143 
CylinderNCylinderN144       __forceinline CylinderN(const Vec3vf<N>& p0, const Vec3vf<N>& p1, const vfloat<N>& rr, bool)
145         : p0(p0), p1(p1), rr(rr) {}
146 
147 
intersectCylinderN148       __forceinline vbool<N> intersect(const Vec3fa& org, const Vec3fa& dir,
149                                        BBox<vfloat<N>>& t_o,
150                                        vfloat<N>& u0_o, Vec3vf<N>& Ng0_o,
151                                        vfloat<N>& u1_o, Vec3vf<N>& Ng1_o) const
152       {
153         /* calculate quadratic equation to solve */
154         const vfloat<N> rl = rcp_length(p1-p0);
155         const Vec3vf<N> P0 = p0, dP = (p1-p0)*rl;
156         const Vec3vf<N> O = Vec3vf<N>(org)-P0, dO = dir;
157 
158         const vfloat<N> dOdO = dot(dO,dO);
159         const vfloat<N> OdO = dot(dO,O);
160         const vfloat<N> OO = dot(O,O);
161         const vfloat<N> dOz = dot(dP,dO);
162         const vfloat<N> Oz = dot(dP,O);
163 
164         const vfloat<N> A = dOdO - sqr(dOz);
165         const vfloat<N> B = 2.0f * (OdO - dOz*Oz);
166         const vfloat<N> C = OO - sqr(Oz) - rr;
167 
168         /* we miss the cylinder if determinant is smaller than zero */
169         const vfloat<N> D = B*B - 4.0f*A*C;
170         vbool<N> valid = D >= 0.0f;
171         if (none(valid)) {
172           t_o = BBox<vfloat<N>>(empty);
173           return valid;
174         }
175 
176         /* standard case for rays that are not parallel to the cylinder */
177         const vfloat<N> Q = sqrt(D);
178         const vfloat<N> rcp_2A = rcp(2.0f*A);
179         const vfloat<N> t0 = (-B-Q)*rcp_2A;
180         const vfloat<N> t1 = (-B+Q)*rcp_2A;
181 
182         /* calculates u and Ng for near hit */
183         {
184           u0_o = madd(t0,dOz,Oz)*rl;
185           const Vec3vf<N> Pr = t0*Vec3vf<N>(dir);
186           const Vec3vf<N> Pl = madd(u0_o,p1-p0,p0);
187           Ng0_o = Pr-Pl;
188         }
189 
190         /* calculates u and Ng for far hit */
191         {
192           u1_o = madd(t1,dOz,Oz)*rl;
193           const Vec3vf<N> Pr = t1*Vec3vf<N>(dir);
194           const Vec3vf<N> Pl = madd(u1_o,p1-p0,p0);
195           Ng1_o = Pr-Pl;
196         }
197 
198         t_o.lower = select(valid, t0, vfloat<N>(pos_inf));
199         t_o.upper = select(valid, t1, vfloat<N>(neg_inf));
200 
201         /* special case for rays that are parallel to the cylinder */
202         const vfloat<N> eps = 16.0f*float(ulp)*max(abs(dOdO),abs(sqr(dOz)));
203         vbool<N> validt = valid & (abs(A) < eps);
204         if (unlikely(any(validt)))
205         {
206           vbool<N> inside = C <= 0.0f;
207           t_o.lower = select(validt,select(inside,vfloat<N>(neg_inf),vfloat<N>(pos_inf)),t_o.lower);
208           t_o.upper = select(validt,select(inside,vfloat<N>(pos_inf),vfloat<N>(neg_inf)),t_o.upper);
209           valid &= !validt | inside;
210         }
211         return valid;
212       }
213 
intersectCylinderN214       __forceinline vbool<N> intersect(const Vec3fa& org_i, const Vec3fa& dir, BBox<vfloat<N>>& t_o) const
215       {
216         vfloat<N> u0_o; Vec3vf<N> Ng0_o;
217         vfloat<N> u1_o; Vec3vf<N> Ng1_o;
218         return intersect(org_i,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o);
219       }
220     };
221   }
222 }
223 
224