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