diff --git a/common/math/linearspace3.h b/common/math/linearspace3.h index f6d2318fa0..6d1ad27b6a 100644 --- a/common/math/linearspace3.h +++ b/common/math/linearspace3.h @@ -113,13 +113,18 @@ namespace embree template __forceinline LinearSpace3 operator +( const LinearSpace3& a ) { return LinearSpace3(+a.vx,+a.vy,+a.vz); } template __forceinline LinearSpace3 rcp ( const LinearSpace3& a ) { return a.inverse(); } - /* constructs a coordinate frame form a normalized normal */ + /* constructs a coordinate frame form a normalized normal + following http://jcgt.org/published/0006/01/01/ */ template __forceinline LinearSpace3 frame(const T& N) { - const T dx0(0,N.z,-N.y); - const T dx1(-N.z,0,N.x); - const T dx = normalize(select(dot(dx0,dx0) > dot(dx1,dx1),dx0,dx1)); - const T dy = normalize(cross(N,dx)); + const typename T::Scalar sgn = sign(N.z); + const typename T::Scalar rs = N.z - sgn; + const typename T::Scalar a = rcp(rs); + const typename T::Scalar rxxs = N.x*N.x + sgn*rs; + const typename T::Scalar ryys = N.y*N.y + sgn*rs; + const typename T::Scalar rxy = N.x*N.y; + const T dx = T(rxxs*a,rxy*a,N.x); + const T dy = T(rxy*a,ryys*a,N.y); return LinearSpace3(dx,dy,N); } diff --git a/kernels/geometry/cone.h b/kernels/geometry/cone.h deleted file mode 100644 index 17429bab32..0000000000 --- a/kernels/geometry/cone.h +++ /dev/null @@ -1,321 +0,0 @@ -// Copyright 2009-2021 Intel Corporation -// SPDX-License-Identifier: Apache-2.0 - -#pragma once - -#include "../common/ray.h" - -namespace embree -{ - namespace isa - { - struct Cone - { - const Vec3fa p0; //!< start position of cone - const Vec3fa p1; //!< end position of cone - const float r0; //!< start radius of cone - const float r1; //!< end radius of cone - - __forceinline Cone(const Vec3fa& p0, const float r0, const Vec3fa& p1, const float r1) - : p0(p0), p1(p1), r0(r0), r1(r1) {} - - __forceinline bool intersect(const Vec3fa& org, const Vec3fa& dir, - BBox1f& t_o, - float& u0_o, Vec3fa& Ng0_o, - float& u1_o, Vec3fa& Ng1_o) const - { - /* calculate quadratic equation to solve */ - const Vec3fa v0 = p0-org; - const Vec3fa v1 = p1-org; - - const float rl = rcp_length(v1-v0); - const Vec3fa P0 = v0, dP = (v1-v0)*rl; - const float dr = (r1-r0)*rl; - const Vec3fa O = -P0, dO = dir; - - const float dOdO = dot(dO,dO); - const float OdO = dot(dO,O); - const float OO = dot(O,O); - const float dOz = dot(dP,dO); - const float Oz = dot(dP,O); - - const float R = r0 + Oz*dr; - const float A = dOdO - sqr(dOz) * (1.0f+sqr(dr)); - const float B = 2.0f * (OdO - dOz*(Oz + R*dr)); - const float C = OO - (sqr(Oz) + sqr(R)); - - /* we miss the cone if determinant is smaller than zero */ - const float D = B*B - 4.0f*A*C; - if (D < 0.0f) return false; - - /* special case for rays that are "parallel" to the cone */ - const float eps = float(1<<8)*float(ulp)*max(abs(dOdO),abs(sqr(dOz))); - if (unlikely(abs(A) < eps)) - { - /* cylinder case */ - if (abs(dr) < 16.0f*float(ulp)) { - if (C <= 0.0f) { t_o = BBox1f(neg_inf,pos_inf); return true; } - else { t_o = BBox1f(pos_inf,neg_inf); return false; } - } - - /* cone case */ - else - { - /* if we hit the negative cone there cannot be a hit */ - const float t = -C/B; - const float z0 = Oz+t*dOz; - const float z0r = r0+z0*dr; - if (z0r < 0.0f) return false; - - /* test if we start inside or outside the cone */ - if (dOz*dr > 0.0f) t_o = BBox1f(t,pos_inf); - else t_o = BBox1f(neg_inf,t); - } - } - - /* standard case for "non-parallel" rays */ - else - { - const float Q = sqrt(D); - const float rcp_2A = rcp(2.0f*A); - t_o.lower = (-B-Q)*rcp_2A; - t_o.upper = (-B+Q)*rcp_2A; - - /* standard case where both hits are on same cone */ - if (likely(A > 0.0f)) { - const float z0 = Oz+t_o.lower*dOz; - const float z0r = r0+z0*dr; - if (z0r < 0.0f) return false; - } - - /* special case where the hits are on the positive and negative cone */ - else - { - /* depending on the ray direction and the open direction - * of the cone we have a hit from inside or outside the - * cone */ - if (dOz*dr > 0) t_o.upper = pos_inf; - else t_o.lower = neg_inf; - } - } - - /* calculates u and Ng for near hit */ - { - u0_o = (Oz+t_o.lower*dOz)*rl; - const Vec3fa Pr = t_o.lower*dir; - const Vec3fa Pl = v0 + u0_o*(v1-v0); - const Vec3fa R = normalize(Pr-Pl); - const Vec3fa U = (p1-p0)+(r1-r0)*R; - const Vec3fa V = cross(p1-p0,R); - Ng0_o = cross(V,U); - } - - /* calculates u and Ng for far hit */ - { - u1_o = (Oz+t_o.upper*dOz)*rl; - const Vec3fa Pr = t_o.upper*dir; - const Vec3fa Pl = v0 + u1_o*(v1-v0); - const Vec3fa R = normalize(Pr-Pl); - const Vec3fa U = (p1-p0)+(r1-r0)*R; - const Vec3fa V = cross(p1-p0,R); - Ng1_o = cross(V,U); - } - return true; - } - - __forceinline bool intersect(const Vec3fa& org, const Vec3fa& dir, BBox1f& t_o) const - { - float u0_o; Vec3fa Ng0_o; float u1_o; Vec3fa Ng1_o; - return intersect(org,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o); - } - - static bool verify(const size_t id, const Cone& cone, const Ray& ray, bool shouldhit, const float t0, const float t1) - { - float eps = 0.001f; - BBox1f t; bool hit; - hit = cone.intersect(ray.org,ray.dir,t); - - bool failed = hit != shouldhit; - if (shouldhit) failed |= std::isinf(t0) ? t0 != t.lower : (t0 == -1E6) ? t.lower > -1E6f : abs(t0-t.lower) > eps; - if (shouldhit) failed |= std::isinf(t1) ? t1 != t.upper : (t1 == +1E6) ? t.upper < +1E6f : abs(t1-t.upper) > eps; - if (!failed) return true; - embree_cout << "Cone test " << id << " failed: cone = " << cone << ", ray = " << ray << ", hit = " << hit << ", t = " << t << embree_endl; - return false; - } - - /* verify cone class */ - static bool verify() - { - bool passed = true; - const Cone cone0(Vec3fa(0.0f,0.0f,0.0f),0.0f,Vec3fa(1.0f,0.0f,0.0f),1.0f); - 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); - 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); - 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); - 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); - 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); - 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); - 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); - 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); - 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); - 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); - const Cone cone1(Vec3fa(0.0f,0.0f,0.0f),1.0f,Vec3fa(1.0f,0.0f,0.0f),0.0f); - 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); - 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); - const Cone cylinder(Vec3fa(0.0f,0.0f,0.0f),1.0f,Vec3fa(1.0f,0.0f,0.0f),1.0f); - 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); - 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); - 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); - 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); - 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); - 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); - 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); - return passed; - } - - /*! output operator */ - friend __forceinline embree_ostream operator<<(embree_ostream cout, const Cone& c) { - return cout << "Cone { p0 = " << c.p0 << ", r0 = " << c.r0 << ", p1 = " << c.p1 << ", r1 = " << c.r1 << "}"; - } - }; - - template - struct ConeN - { - typedef Vec3> Vec3vfN; - - const Vec3vfN p0; //!< start position of cone - const Vec3vfN p1; //!< end position of cone - const vfloat r0; //!< start radius of cone - const vfloat r1; //!< end radius of cone - - __forceinline ConeN(const Vec3vfN& p0, const vfloat& r0, const Vec3vfN& p1, const vfloat& r1) - : p0(p0), p1(p1), r0(r0), r1(r1) {} - - __forceinline Cone operator[] (const size_t i) const - { - assert(i intersect(const Vec3fa& org, const Vec3fa& dir, - BBox>& t_o, - vfloat& u0_o, Vec3vfN& Ng0_o, - vfloat& u1_o, Vec3vfN& Ng1_o) const - { - /* calculate quadratic equation to solve */ - const Vec3vfN v0 = p0-Vec3vfN(org); - const Vec3vfN v1 = p1-Vec3vfN(org); - - const vfloat rl = rcp_length(v1-v0); - const Vec3vfN P0 = v0, dP = (v1-v0)*rl; - const vfloat dr = (r1-r0)*rl; - const Vec3vfN O = -P0, dO = dir; - - const vfloat dOdO = dot(dO,dO); - const vfloat OdO = dot(dO,O); - const vfloat OO = dot(O,O); - const vfloat dOz = dot(dP,dO); - const vfloat Oz = dot(dP,O); - - const vfloat R = r0 + Oz*dr; - const vfloat A = dOdO - sqr(dOz) * (vfloat(1.0f)+sqr(dr)); - const vfloat B = 2.0f * (OdO - dOz*(Oz + R*dr)); - const vfloat C = OO - (sqr(Oz) + sqr(R)); - - /* we miss the cone if determinant is smaller than zero */ - const vfloat D = B*B - 4.0f*A*C; - vbool valid = D >= 0.0f; - if (none(valid)) return valid; - - /* special case for rays that are "parallel" to the cone */ - const vfloat eps = float(1<<8)*float(ulp)*max(abs(dOdO),abs(sqr(dOz))); - const vbool validt = valid & (abs(A) < eps); - const vbool validf = valid & !(abs(A) < eps); - if (unlikely(any(validt))) - { - const vboolx validtt = validt & (abs(dr) < 16.0f*float(ulp)); - const vboolx validtf = validt & (abs(dr) >= 16.0f*float(ulp)); - - /* cylinder case */ - if (unlikely(any(validtt))) - { - t_o.lower = select(validtt, select(C <= 0.0f, vfloat(neg_inf), vfloat(pos_inf)), t_o.lower); - t_o.upper = select(validtt, select(C <= 0.0f, vfloat(pos_inf), vfloat(neg_inf)), t_o.upper); - valid &= !validtt | C <= 0.0f; - } - - /* cone case */ - if (any(validtf)) - { - /* if we hit the negative cone there cannot be a hit */ - const vfloat t = -C/B; - const vfloat z0 = Oz+t*dOz; - const vfloat z0r = r0+z0*dr; - valid &= !validtf | z0r >= 0.0f; - - /* test if we start inside or outside the cone */ - t_o.lower = select(validtf, select(dOz*dr > 0.0f, t, vfloat(neg_inf)), t_o.lower); - t_o.upper = select(validtf, select(dOz*dr > 0.0f, vfloat(pos_inf), t), t_o.upper); - } - } - - /* standard case for "non-parallel" rays */ - if (likely(any(validf))) - { - const vfloat Q = sqrt(D); - const vfloat rcp_2A = 0.5f*rcp(A); - t_o.lower = select(validf, (-B-Q)*rcp_2A, t_o.lower); - t_o.upper = select(validf, (-B+Q)*rcp_2A, t_o.upper); - - /* standard case where both hits are on same cone */ - const vbool validft = validf & A>0.0f; - const vbool validff = validf & !(A>0.0f); - if (any(validft)) { - const vfloat z0 = Oz+t_o.lower*dOz; - const vfloat z0r = r0+z0*dr; - valid &= !validft | z0r >= 0.0f; - } - - /* special case where the hits are on the positive and negative cone */ - if (any(validff)) { - /* depending on the ray direction and the open direction - * of the cone we have a hit from inside or outside the - * cone */ - t_o.lower = select(validff, select(dOz*dr > 0.0f, t_o.lower, float(neg_inf)), t_o.lower); - t_o.upper = select(validff, select(dOz*dr > 0.0f, float(pos_inf), t_o.upper), t_o.upper); - } - } - - /* calculates u and Ng for near hit */ - { - u0_o = (Oz+t_o.lower*dOz)*rl; - const Vec3vfN Pr = t_o.lower*Vec3vfN(dir); - const Vec3vfN Pl = v0 + u0_o*(v1-v0); - const Vec3vfN R = normalize(Pr-Pl); - const Vec3vfN U = (p1-p0)+(r1-r0)*R; - const Vec3vfN V = cross(p1-p0,R); - Ng0_o = cross(V,U); - } - - /* calculates u and Ng for far hit */ - { - u1_o = (Oz+t_o.upper*dOz)*rl; - const Vec3vfN Pr = t_o.lower*Vec3vfN(dir); - const Vec3vfN Pl = v0 + u1_o*(v1-v0); - const Vec3vfN R = normalize(Pr-Pl); - const Vec3vfN U = (p1-p0)+(r1-r0)*R; - const Vec3vfN V = cross(p1-p0,R); - Ng1_o = cross(V,U); - } - return valid; - } - - __forceinline vbool intersect(const Vec3fa& org, const Vec3fa& dir, BBox>& t_o) const - { - vfloat u0_o; Vec3vfN Ng0_o; vfloat u1_o; Vec3vfN Ng1_o; - return intersect(org,dir,t_o,u0_o,Ng0_o,u1_o,Ng1_o); - } - }; - } -} - diff --git a/kernels/geometry/coneline_intersector.h b/kernels/geometry/coneline_intersector.h index 696ea41ebc..9b8621e577 100644 --- a/kernels/geometry/coneline_intersector.h +++ b/kernels/geometry/coneline_intersector.h @@ -10,16 +10,16 @@ namespace embree { namespace isa { - namespace __coneline_internal + namespace __coneline_internal { template static __forceinline bool intersectCone(const vbool& valid_i, - const Vec3vf& ray_org_in, const Vec3vf& ray_dir, + const Vec3vf& ray_org_in, const Vec3vf& ray_dir, const vfloat& ray_tnear, const ray_tfar_func& ray_tfar, const Vec4vf& v0, const Vec4vf& v1, const vbool& cL, const vbool& cR, const Epilog& epilog) - { + { vbool valid = valid_i; /* move ray origin closer to make calculations numerically stable */ @@ -32,10 +32,10 @@ namespace embree const Vec3vf dP = v1.xyz() - v0.xyz(); const Vec3vf p0 = ray_org - v0.xyz(); const Vec3vf p1 = ray_org - v1.xyz(); - + const vfloat dPdP = sqr(dP); const vfloat dP0 = dot(p0,dP); - const vfloat dP1 = dot(p1,dP); + const vfloat dP1 = dot(p1,dP); const vfloat dOdP = dot(ray_dir,dP); // intersect cone body @@ -45,11 +45,11 @@ namespace embree const vfloat OO = sqr(p0); const vfloat dPdP2 = sqr(dPdP); const vfloat dPdPr0 = dPdP*v0.w; - + const vfloat A = dPdP2 - sqr(dOdP)*hy; const vfloat B = dPdP2*dO0 - dP0*dOdP*hy + dPdPr0*(dr*dOdP); const vfloat C = dPdP2*OO - sqr(dP0)*hy + dPdPr0*(2.0f*dr*dP0 - dPdPr0); - + const vfloat D = B*B - A*C; valid &= D >= 0.0f; if (unlikely(none(valid))) { @@ -78,8 +78,8 @@ namespace embree const vfloat t_disk_upper = max(t_disk0, t_disk1); const vfloat t_lower = min(t_cone_lower, t_disk_lower); - const vfloat t_upper = max(t_cone_upper, select(t_lower==t_disk_lower, - select(t_disk_upper==vfloat(pos_inf),neg_inf,t_disk_upper), + const vfloat t_upper = max(t_cone_upper, select(t_lower==t_disk_lower, + select(t_disk_upper==vfloat(pos_inf),neg_inf,t_disk_upper), select(t_disk_lower==vfloat(pos_inf),neg_inf,t_disk_lower))); const vbool valid_lower = valid & ray_tnear <= dt+t_lower & dt+t_lower <= ray_tfar() & t_lower != vfloat(pos_inf); @@ -110,7 +110,7 @@ namespace embree const vbool valid_second = valid_lower & valid_upper & (dt+t_upper <= ray_tfar()); if (unlikely(none(valid_second))) return is_hit_first; - + /* invoke intersection filter for second hit */ const vbool cone_hit_second = t_second == t_cone_lower | t_second == t_cone_upper; const vbool disk0_hit_second = t_second == t_disk0; @@ -119,7 +119,7 @@ namespace embree hit = RoundLineIntersectorHitM(u_second,zero,dt+t_second,Ng_second); const bool is_hit_second = epilog(valid_second, hit); - + return is_hit_first | is_hit_second; } } @@ -128,7 +128,7 @@ namespace embree struct ConeLineIntersectorHitM { __forceinline ConeLineIntersectorHitM() {} - + __forceinline ConeLineIntersectorHitM(const vfloat& u, const vfloat& v, const vfloat& t, const Vec3vf& Ng) : vu(u), vv(v), vt(t), vNg(Ng) {} @@ -144,12 +144,12 @@ namespace embree vfloat vt; Vec3vf vNg; }; - + template struct ConeCurveIntersector1 { typedef CurvePrecalculations1 Precalculations; - + struct ray_tfar { Ray& ray; __forceinline ray_tfar(Ray& ray) : ray(ray) {} @@ -174,19 +174,19 @@ namespace embree return __coneline_internal::intersectCone(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray),v0,v1,cL,cR,epilog); } }; - + template struct ConeCurveIntersectorK { typedef CurvePrecalculationsK Precalculations; - + struct ray_tfar { RayK& ray; size_t k; __forceinline ray_tfar(RayK& ray, size_t k) : ray(ray), k(k) {} __forceinline vfloat operator() () const { return ray.tfar[k]; }; }; - + template static __forceinline bool intersect(const vbool& valid_i, RayK& ray, size_t k, diff --git a/kernels/geometry/roundline_intersector.h b/kernels/geometry/roundline_intersector.h index a83dd72a7f..f3ffc2e292 100644 --- a/kernels/geometry/roundline_intersector.h +++ b/kernels/geometry/roundline_intersector.h @@ -8,7 +8,6 @@ /* - This file implements the intersection of a ray with a round linear curve segment. We define the geometry of such a round linear curve segment from point p0 with radius r0 to point p1 with radius r1 @@ -17,6 +16,14 @@ p0/r0 to p1/r1 with cone(p0,r0,p1,r1) and the cone plus the ending sphere with cone_sphere(p0,r0,p1,r1). + The method follows Reshetov and Hart "Modeling Hair Strands with + Roving Capsules", SIGGRAPH 2024, + https://dl.acm.org/doi/10.1145/3641519.3657450 + See also https://www.shadertoy.com/view/4ffXWs + The surface is directly described by a sphere (with changing radius) + swept along a line. First we solve for u, where the ray will intersect + the sphere, followed by a sphere intersection. + For multiple connected round linear curve segments this construction yield a proper shape when viewed from the outside. Using the following CSG we can also handle the interior in most common cases: @@ -29,7 +36,7 @@ combined curve. This approach works as long as geometry of the current cone_sphere penetrates into direct neighbor segments only, and not into segments further away. - + To construct a cone that touches two spheres at p0 and p1 with r0 and r1, one has to increase the cone radius at r0 and r1 to obtain larger radii w0 and w1, such that the infinite cone properly touches @@ -72,7 +79,7 @@ namespace embree struct RoundLineIntersectorHitM { __forceinline RoundLineIntersectorHitM() {} - + __forceinline RoundLineIntersectorHitM(const vfloat& u, const vfloat& v, const vfloat& t, const Vec3vf& Ng) : vu(u), vv(v), vt(t), vNg(Ng) {} @@ -85,14 +92,14 @@ namespace embree __forceinline Vec2vf uv() const { return Vec2vf(vu,vv); } __forceinline vfloat t () const { return vt; } __forceinline Vec3vf Ng() const { return vNg; } - + public: vfloat vu; vfloat vv; vfloat vt; Vec3vf vNg; }; - + namespace __roundline_internal { template @@ -100,116 +107,75 @@ namespace embree { ConeGeometry (const Vec4vf& a, const Vec4vf& b) : p0(a.xyz()), p1(b.xyz()), dP(p1-p0), dPdP(dot(dP,dP)), r0(a.w), sqr_r0(sqr(r0)), r1(b.w), dr(r1-r0), drdr(dr*dr), r0dr (r0*dr), g(dPdP - drdr) {} - - /* - - This function tests if a point is accepted by first cone - clipping plane. - - First, we need to project the point onto the line p0->p1: - - Y = (p-p0)*(p1-p0)/length(p1-p0) - - This value y is the distance to the projection point from - p0. The clip distances are calculated as: - - Y0 = - r0 * (r1-r0) / length(p1-p0) - Y1 = length(p1-p0) - r1 * (r1-r0) / length(p1-p0) - - Thus to test if the point p is accepted by the first - clipping plane we need to test Y > Y0 and to test if it - is accepted by the second clipping plane we need to test - Y < Y1. - - By multiplying the calculations with length(p1-p0) these - calculation can get simplied to: - - y = (p-p0)*(p1-p0) - y0 = - r0 * (r1-r0) - y1 = (p1-p0)^2 - r1 * (r1-r0) - - and the test y > y0 and y < y1. - - */ - - __forceinline vbool isClippedByPlane (const vbool& valid_i, const Vec3vf& p) const - { - const Vec3vf p0p = p - p0; - const vfloat y = dot(p0p,dP); - const vfloat cap0 = -r0dr; - const vbool inside_cone = y > cap0; - return valid_i & (p0.x != vfloat(inf)) & (p1.x != vfloat(inf)) & inside_cone; - } - - /* - + + /* This function tests whether a point lies inside the capped cone tangential to its ending spheres. Therefore one has to check if the point is inside the region defined by the cone clipping planes, which is performed similar as in the previous function. - + To perform the inside cone test we need to project the point onto the line p0->p1: - + dP = p1-p0 Y = (p-p0)*dP/length(dP) - + This value Y is the distance to the projection point from p0. To obtain a parameter value u going from 0 to 1 along the line p0->p1 we calculate: - + U = Y/length(dP) - + The radii to use at points p0 and p1 are: - + w0 = sr * r0 w1 = sr * r1 dw = w1-w0 - + Using these radii and u one can directly test if the point lies inside the cone using the formula dP*dP < wy*wy with: - + wy = w0 + u*dw py = p0 + u*dP - p - + By multiplying the calculations with length(p1-p0) and inserting the definition of w can obtain simpler equations: - + y = (p-p0)*dP ry = r0 + y/dP^2 * dr - wy = sr*ry + wy = sr*ry py = p0 + y/dP^2*dP - p y0 = - r0 * dr y1 = dP^2 - r1 * dr - + Thus for the in-cone test we get: - + py^2 < wy^2 <=> py^2 < sr^2 * ry^2 <=> py^2 * ( dP^2 - dr^2 ) < dP^2 * ry^2 - + This can further get simplified to: - - (p0-p)^2 * (dP^2 - dr^2) - y^2 < dP^2 * r0^2 + 2.0f*r0*dr*y; - + + (p0-p)^2 * (dP^2 - dr^2) - y^2 < dP^2 * r0^2 + 2.0f*r0*dr*y; */ - - __forceinline vbool isInsideCappedCone (const vbool& valid_i, const Vec3vf& p) const + + __forceinline vbool isInsideCappedCone(const vbool& valid_i, const vfloat& t) const { - const Vec3vf p0p = p - p0; + Vec3vf p0p = -p0; + p0p.z += t; const vfloat y = dot(p0p,dP); - const vfloat cap0 = -r0dr+vfloat(ulp); - const vfloat cap1 = -r1*dr + dPdP; - + const vfloat cap0 = vfloat(ulp) - r0dr; + const vfloat cap1 = dPdP - r1*dr; + vbool inside_cone = valid_i & (p0.x != vfloat(inf)) & (p1.x != vfloat(inf)); inside_cone &= y > cap0; // start clipping plane - inside_cone &= y < cap1; // end clipping plane + inside_cone &= y < cap1; // end clipping plane inside_cone &= sqr(p0p)*g - sqr(y) < dPdP * sqr_r0 + 2.0f*r0dr*y; // in cone test return inside_cone; } - + protected: Vec3vf p0; Vec3vf p1; @@ -223,429 +189,201 @@ namespace embree vfloat r0dr; vfloat g; }; - + template struct ConeGeometryIntersector : public ConeGeometry { using ConeGeometry::p0; - using ConeGeometry::p1; using ConeGeometry::dP; - using ConeGeometry::dPdP; using ConeGeometry::r0; using ConeGeometry::sqr_r0; - using ConeGeometry::r1; using ConeGeometry::dr; using ConeGeometry::r0dr; + using ConeGeometry::drdr; using ConeGeometry::g; - - ConeGeometryIntersector (const Vec3vf& ray_org, const Vec3vf& ray_dir, const vfloat& dOdO, const vfloat& rcp_dOdO, const Vec4vf& a, const Vec4vf& b) - : ConeGeometry(a,b), org(ray_org), O(ray_org-p0), dO(ray_dir), dOdO(dOdO), rcp_dOdO(rcp_dOdO), OdP(dot(dP,O)), dOdP(dot(dP,dO)), yp(OdP + r0dr) {} - + + ConeGeometryIntersector(const Vec4vf& a, const Vec4vf& b) + : ConeGeometry(a,b) { + G = sqr_r0 - sqr(p0.x) - sqr(p0.y); + C = sqr(dP.x) + sqr(dP.y) - drdr; + E = sqr(dP.z) + C; + F = r0dr - p0.x*dP.x - p0.y*dP.y; + } + /* - - This function intersects a ray with a cone that touches a - start sphere p0/r0 and end sphere p1/r1. - - To find this ray/cone intersections one could just - calculate radii w0 and w1 as described above and use a - standard ray/cone intersection routine with these - radii. However, it turns out that calculations can get - simplified when deriving a specialized ray/cone - intersection for this special case. We perform - calculations relative to the cone origin p0 and define: - - O = ray_org - p0 - dO = ray_dir - dP = p1-p0 - dr = r1-r0 - dw = w1-w0 - - For some t we can compute the potential hit point h = O + t*dO and - project it onto the cone vector dP to obtain u = (h*dP)/(dP*dP). In - case of an intersection, the squared distance from the hit point - projected onto the cone center line to the hit point should be equal - to the squared cone radius at u: - - (u*dP - h)^2 = (w0 + u*dw)^2 - - Inserting the definition of h, u, w0, and dw into this formula, then - factoring out all terms, and sorting by t^2, t^1, and t^0 terms - yields a quadratic equation to solve. - - Inserting u: - ( (h*dP)*dP/dP^2 - h )^2 = ( w0 + (h*dP)*dw/dP^2 )^2 - - Multiplying by dP^4: - ( (h*dP)*dP - h*dP^2 )^2 = ( w0*dP^2 + (h*dP)*dw )^2 - - Inserting w0 and dw: - ( (h*dP)*dP - h*dP^2 )^2 = ( r0*dP^2 + (h*dP)*dr )^2 / (1-dr^2/dP^2) - ( (h*dP)*dP - h*dP^2 )^2 *(dP^2 - dr^2) = dP^2 * ( r0*dP^2 + (h*dP)*dr )^2 - - Now one can insert the definition of h, factor out, and presort by t: - ( ((O + t*dO)*dP)*dP - (O + t*dO)*dP^2 )^2 *(dP^2 - dr^2) = dP^2 * ( r0*dP^2 + ((O + t*dO)*dP)*dr )^2 - ( (O*dP)*dP-O*dP^2 + t*( (dO*dP)*dP - dO*dP^2 ) )^2 *(dP^2 - dr^2) = dP^2 * ( r0*dP^2 + (O*dP)*dr + t*(dO*dP)*dr )^2 - - Factoring out further and sorting by t^2, t^1 and t^0 yields: - - 0 = t^2 * [ ((dO*dP)*dP - dO-dP^2)^2 * (dP^2 - dr^2) - dP^2*(dO*dP)^2*dr^2 ] - + 2*t^1 * [ ((O*dP)*dP - O*dP^2) * ((dO*dP)*dP - dO*dP^2) * (dP^2 - dr^2) - dP^2*(r0*dP^2 + (O*dP)*dr)*(dO*dP)*dr ] - + t^0 * [ ( (O*dP)*dP - O*dP^2)^2 * (dP^2-dr^2) - dP^2*(r0*dP^2 + (O*dP)*dr)^2 ] - - This can be simplified to: - - 0 = t^2 * [ (dP^2 - dr^2)*dO^2 - (dO*dP)^2 ] - + 2*t^1 * [ (dP^2 - dr^2)*(O*dO) - (dO*dP)*(O*dP + r0*dr) ] - + t^0 * [ (dP^2 - dr^2)*O^2 - (O*dP)^2 - r0^2*dP^2 - 2.0f*r0*dr*(O*dP) ] - - Solving this quadratic equation yields the values for t at which the - ray intersects the cone. - + + This function finds u where the ray will intersect the swept sphere + */ - - __forceinline bool intersectCone(vbool& valid, vfloat& lower, vfloat& upper) + + __forceinline bool intersectConeU(vbool& valid, vfloat& u_front, vfloat& u_back) { - /* return no hit by default */ - lower = pos_inf; - upper = neg_inf; - - /* compute quadratic equation A*t^2 + B*t + C = 0 */ - const vfloat OO = dot(O,O); - const vfloat OdO = dot(dO,O); - const vfloat A = g * dOdO - sqr(dOdP); - const vfloat B = 2.0f * (g*OdO - dOdP*yp); - const vfloat C = g*OO - sqr(OdP) - sqr_r0*dPdP - 2.0f*r0dr*OdP; - /* we miss the cone if determinant is smaller than zero */ - const vfloat D = B*B - 4.0f*A*C; - valid &= (D >= 0.0f & g > 0.0f); // if g <= 0 then the cone is inside a sphere end - - /* When rays are parallel to the cone surface, then the - * ray may be inside or outside the cone. We just assume a - * miss in that case, which is fine as rays inside the - * cone would anyway hit the ending spheres in that - * case. */ - valid &= abs(A) > min_rcp_input; - if (unlikely(none(valid))) { + const vfloat D = (F*F + G*C) * rcp(E); + valid &= (D >= 0.0f) | (g <= 0.0f); // or inside a sphere end + + if (unlikely(none(valid))) return false; - } - - /* compute distance to front and back hit */ - const vfloat Q = sqrt(D); - const vfloat rcp_2A = rcp(2.0f*A); - t_cone_front = (-B-Q)*rcp_2A; - y_cone_front = yp + t_cone_front*dOdP; - lower = select( (y_cone_front > -(float)ulp) & (y_cone_front <= g) & (g > 0.0f), t_cone_front, vfloat(pos_inf)); + + const vfloat Q = dP.z * sqrt(D); + const vfloat rcp_C = rcp(C); + vfloat u0 = (F + Q) * rcp_C; + vfloat u1 = (F - Q) * rcp_C; + // flip negative caps to other side + u0 = select(r0 + u0*dr < 0.0f, 1.0f-u0, u0); + u1 = select(r0 + u1*dr < 0.0f, 1.0f-u1, u1); + const vfloat uend = select(dr >= 0.0f, vfloat(one), vfloat(zero)); + const vbool swap = (dP.z >= 0.0) != (u1 > u0); + // already clamp for end sphere (always present) + u_front = min(1.0f, select(g > 0.0f, select(swap, u1, u0), uend)); #if !defined (EMBREE_BACKFACE_CULLING_CURVES) - t_cone_back = (-B+Q)*rcp_2A; - y_cone_back = yp + t_cone_back *dOdP; - upper = select( (y_cone_back > -(float)ulp) & (y_cone_back <= g) & (g > 0.0f), t_cone_back , vfloat(neg_inf)); -#endif - return true; - } - - /* - This function intersects the ray with the end sphere at - p1. We already clip away hits that are inside the - neighboring cone segment. - - */ - - __forceinline void intersectEndSphere(vbool& valid, - const ConeGeometry& coneR, - vfloat& lower, vfloat& upper) - { - /* calculate front and back hit with end sphere */ - const Vec3vf O1 = org - p1; - const vfloat O1dO = dot(O1,dO); - const vfloat h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r1)); - const vfloat rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat(neg_inf) ); - - /* clip away front hit if it is inside next cone segment */ - t_sph1_front = (-O1dO - rhs1)*rcp_dOdO; - const Vec3vf hit_front = org + t_sph1_front*dO; - vbool valid_sph1_front = h2 >= 0.0f & yp + t_sph1_front*dOdP > g & !coneR.isClippedByPlane (valid, hit_front); - lower = select(valid_sph1_front, t_sph1_front, vfloat(pos_inf)); - -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - /* clip away back hit if it is inside next cone segment */ - t_sph1_back = (-O1dO + rhs1)*rcp_dOdO; - const Vec3vf hit_back = org + t_sph1_back*dO; - vbool valid_sph1_back = h2 >= 0.0f & yp + t_sph1_back*dOdP > g & !coneR.isClippedByPlane (valid, hit_back); - upper = select(valid_sph1_back, t_sph1_back, vfloat(neg_inf)); -#else - upper = vfloat(neg_inf); + u_back = min(1.0f, select(g > 0.0f, select(swap, u0, u1), uend)); #endif + return true; } - __forceinline void intersectBeginSphere(const vbool& valid, - vfloat& lower, vfloat& upper) - { - /* calculate front and back hit with end sphere */ - const Vec3vf O1 = org - p0; - const vfloat O1dO = dot(O1,dO); - const vfloat h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r0)); - const vfloat rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat(neg_inf) ); - - /* clip away front hit if it is inside next cone segment */ - t_sph0_front = (-O1dO - rhs1)*rcp_dOdO; - vbool valid_sph1_front = valid & h2 >= 0.0f & yp + t_sph0_front*dOdP < 0; - lower = select(valid_sph1_front, t_sph0_front, vfloat(pos_inf)); - -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - /* clip away back hit if it is inside next cone segment */ - t_sph0_back = (-O1dO + rhs1)*rcp_dOdO; - vbool valid_sph1_back = valid & h2 >= 0.0f & yp + t_sph0_back*dOdP < 0; - upper = select(valid_sph1_back, t_sph0_back, vfloat(neg_inf)); -#else - upper = vfloat(neg_inf); -#endif - } - - /* - - This function calculates the geometry normal of some cone hit. - - For a given hit point h (relative to p0) with a cone - starting at p0 with radius w0 and ending at p1 with - radius w1 one normally calculates the geometry normal by - first calculating the parmetric u hit location along the - cone: - - u = dot(h,dP)/dP^2 - - Using this value one can now directly calculate the - geometry normal by bending the connection vector (h-u*dP) - from hit to projected hit with some cone dependent value - dw/sqrt(dP^2) * normalize(dP): - - Ng = normalize(h-u*dP) - dw/length(dP) * normalize(dP) - - The length of the vector (h-u*dP) can also get calculated - by interpolating the radii as w0+u*dw which yields: - - Ng = (h-u*dP)/(w0+u*dw) - dw/dP^2 * dP - - Multiplying with (w0+u*dw) yield a scaled Ng': - - Ng' = (h-u*dP) - (w0+u*dw)*dw/dP^2*dP - - Inserting the definition of w0 and dw and refactoring - yield a further scaled Ng'': - - Ng'' = (dP^2 - dr^2) (h-q) - (r0+u*dr)*dr*dP - - Now inserting the definition of u gives and multiplying - with the denominator yields: - - Ng''' = (dP^2-dr^2)*(dP^2*h-dot(h,dP)*dP) - (dP^2*r0+dot(h,dP)*dr)*dr*dP - - Factoring out, cancelling terms, dividing by dP^2, and - factoring again yields finally: - - Ng'''' = (dP^2-dr^2)*h - dP*(dot(h,dP) + r0*dr) - + /* + This function intersects the ray with the swept sphere centered at u. */ - - __forceinline Vec3vf Ng_cone(const vbool& front_hit) const - { -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - const vfloat y = select(front_hit, y_cone_front, y_cone_back); - const vfloat t = select(front_hit, t_cone_front, t_cone_back); - const Vec3vf h = O + t*dO; - return g*h-dP*y; -#else - const Vec3vf h = O + t_cone_front*dO; - return g*h-dP*y_cone_front; -#endif - } - - /* compute geometry normal of sphere hit as the difference - * vector from hit point to sphere center */ - - __forceinline Vec3vf Ng_sphere1(const vbool& front_hit) const + template + __forceinline void intersectSphere(vbool& valid, const vfloat& u, vfloat& t) { -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - const vfloat t_sph1 = select(front_hit, t_sph1_front, t_sph1_back); - return org+t_sph1*dO-p1; -#else - return org+t_sph1_front*dO-p1; -#endif + const vfloat D = G + u * (2.0 * F - u * C); + valid &= D >= 0.0f; + const vfloat d = sqrt(D); + t = p0.z + u * dP.z; + if (front) + t -= d; + else + t += d; } - __forceinline Vec3vf Ng_sphere0(const vbool& front_hit) const - { -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - const vfloat t_sph0 = select(front_hit, t_sph0_front, t_sph0_back); - return org+t_sph0*dO-p0; -#else - return org+t_sph0_front*dO-p0; -#endif - } - - /* - This function calculates the u coordinate of a - hit. Therefore we use the hit distance y (which is zero - at the first cone clipping plane) and divide by distance - g between the clipping planes. - + /* + Compute geometry normal of hit as the difference vector from hit + point to projected hit, which is the interpolated sphere center. + This normal is valid for the cone as well for the end cap spheres. + + Ng = h - (p0 + u*dP) + h = (0,0,t) */ - - __forceinline vfloat u_cone(const vbool& front_hit) const + + __forceinline Vec3vf Ng(const vfloat& u, const vfloat& t) const { -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - const vfloat y = select(front_hit, y_cone_front, y_cone_back); - return clamp(y*rcp(g)); -#else - return clamp(y_cone_front*rcp(g)); -#endif + Vec3vf n = -(p0 + u*dP); + n.z += t; + return n; } - - private: - Vec3vf org; - Vec3vf O; - Vec3vf dO; - vfloat dOdO; - vfloat rcp_dOdO; - vfloat OdP; - vfloat dOdP; - - /* for ray/cone intersection */ - private: - vfloat yp; - vfloat y_cone_front; - vfloat t_cone_front; -#if !defined (EMBREE_BACKFACE_CULLING_CURVES) - vfloat y_cone_back; - vfloat t_cone_back; -#endif - - /* for ray/sphere intersection */ + private: - vfloat t_sph1_front; - vfloat t_sph0_front; -#if !defined (EMBREE_BACKFACE_CULLING_CURVES) - vfloat t_sph1_back; - vfloat t_sph0_back; -#endif + vfloat C; + vfloat E; + vfloat F; + vfloat G; }; - - + template static __forceinline bool intersectConeSphere(const vbool& valid_i, - const Vec3vf& ray_org_in, const Vec3vf& ray_dir, + const Vec3vf& ray_org, const Vec3vf& ray_dir_in, const vfloat& ray_tnear, const ray_tfar_func& ray_tfar, - const Vec4vf& v0, const Vec4vf& v1, - const Vec4vf& vL, const Vec4vf& vR, + const Vec4vf& v0_in, const Vec4vf& v1_in, + const Vec4vf& vL_in, const Vec4vf& vR_in, const Epilog& epilog) - { + { vbool valid = valid_i; - - /* move ray origin closer to make calculations numerically stable */ - const vfloat dOdO = sqr(ray_dir); - const vfloat rcp_dOdO = rcp(dOdO); - const Vec3vf center = vfloat(0.5f)*(v0.xyz()+v1.xyz()); - const vfloat dt = dot(center-ray_org_in,ray_dir)*rcp_dOdO; - const Vec3vf ray_org = ray_org_in + dt*ray_dir; - - /* intersect with cone from v0 to v1 */ - vfloat t_cone_lower, t_cone_upper; - ConeGeometryIntersector cone (ray_org, ray_dir, dOdO, rcp_dOdO, v0, v1); - vbool validCone = valid; - cone.intersectCone(validCone, t_cone_lower, t_cone_upper); - valid &= (validCone | (cone.g <= 0.0f)); // if cone is entirely in sphere end - check sphere - if (unlikely(none(valid))) + // normalize ray dir, keep length + const vfloat rcp_rd = rcp_length(ray_dir_in); + const vfloat rd = rcp(rcp_rd); + const Vec3vf& ray_dir = ray_dir_in*rcp_rd; + + // transform into ray space + const LinearSpace3vf m = frame(ray_dir); + const Vec4vf v0 = Vec4vf(m*(Vec3vf(v0_in)-ray_org), v0_in.w); + const Vec4vf v1 = Vec4vf(m*(Vec3vf(v1_in)-ray_org), v1_in.w); + + /* intersect with cone from v0 to v1 */ + ConeGeometryIntersector cone(v0, v1); + vfloat u_front, u_back; + if (unlikely(!cone.intersectConeU(valid, u_front, u_back))) return false; - - /* cone hits inside the neighboring capped cones are inside the geometry and thus ignored */ + + /* calculate front hit */ + vbool validFront = valid & ((u_front >= 0.0f) | (vL_in[0] == vfloat(pos_inf))); + u_front = max(u_front, 0.0f); + vfloat t_lower; + cone.template intersectSphere(validFront, u_front, t_lower); + + // for normal + Vec3vf dP = v0_in.xyz() - v1_in.xyz(); + Vec3vf P = ray_org - v0_in.xyz(); + +#if !defined (EMBREE_BACKFACE_CULLING_CURVES) + /* hits inside the neighboring capped cones are inside the geometry and thus ignored */ + const Vec4vf vL = Vec4vf(m*(Vec3vf(vL_in)-ray_org), vL_in.w); const ConeGeometry coneL (v0, vL); + const Vec4vf vR = Vec4vf(m*(Vec3vf(vR_in)-ray_org), vR_in.w); const ConeGeometry coneR (v1, vR); -#if !defined(EMBREE_BACKFACE_CULLING_CURVES) - const Vec3vf hit_lower = ray_org + t_cone_lower*ray_dir; - const Vec3vf hit_upper = ray_org + t_cone_upper*ray_dir; - t_cone_lower = select (!coneL.isInsideCappedCone (validCone, hit_lower) & !coneR.isInsideCappedCone (validCone, hit_lower), t_cone_lower, vfloat(pos_inf)); - t_cone_upper = select (!coneL.isInsideCappedCone (validCone, hit_upper) & !coneR.isInsideCappedCone (validCone, hit_upper), t_cone_upper, vfloat(neg_inf)); -#endif + validFront &= !coneL.isInsideCappedCone(validFront, t_lower) & !coneR.isInsideCappedCone(validFront, t_lower); - /* intersect ending sphere */ - vfloat t_sph1_lower, t_sph1_upper; - vfloat t_sph0_lower = vfloat(pos_inf); - vfloat t_sph0_upper = vfloat(neg_inf); - cone.intersectEndSphere(valid, coneR, t_sph1_lower, t_sph1_upper); + /* calculate back hit */ + vbool validBack = valid & ((u_back >= 0.0f) | (vL_in[0] == vfloat(pos_inf))); + u_back = max(u_back, 0.0f); + vfloat t_upper; + cone.template intersectSphere(validBack, u_back, t_upper); + validBack &= !coneL.isInsideCappedCone(validBack, t_upper) & !coneR.isInsideCappedCone(validBack, t_upper); + + // scale back + t_lower *= rd; + t_upper *= rd; - const vbool isBeginPoint = valid & (vL[0] == vfloat(pos_inf)); - if (unlikely(any(isBeginPoint))) { - cone.intersectBeginSphere (isBeginPoint, t_sph0_lower, t_sph0_upper); - } - - /* CSG union of cone and end sphere */ - vfloat t_sph_lower = min(t_sph0_lower, t_sph1_lower); - vfloat t_cone_sphere_lower = min(t_cone_lower, t_sph_lower); -#if !defined (EMBREE_BACKFACE_CULLING_CURVES) - vfloat t_sph_upper = max(t_sph0_upper, t_sph1_upper); - vfloat t_cone_sphere_upper = max(t_cone_upper, t_sph_upper); - /* filter out hits that are not in tnear/tfar range */ - const vbool valid_lower = valid & ray_tnear <= dt+t_cone_sphere_lower & dt+t_cone_sphere_lower <= ray_tfar() & t_cone_sphere_lower != vfloat(pos_inf); - const vbool valid_upper = valid & ray_tnear <= dt+t_cone_sphere_upper & dt+t_cone_sphere_upper <= ray_tfar() & t_cone_sphere_upper != vfloat(neg_inf); - + const vbool valid_lower = validFront & ray_tnear <= t_lower & t_lower <= ray_tfar(); + const vbool valid_upper = validBack & ray_tnear <= t_upper & t_upper <= ray_tfar(); + /* check if there is a first hit */ const vbool valid_first = valid_lower | valid_upper; if (unlikely(none(valid_first))) return false; - + /* construct first hit */ - const vfloat t_first = select(valid_lower, t_cone_sphere_lower, t_cone_sphere_upper); - const vbool cone_hit_first = t_first == t_cone_lower | t_first == t_cone_upper; - const vbool sph0_hit_first = t_first == t_sph0_lower | t_first == t_sph0_upper; - const Vec3vf Ng_first = select(cone_hit_first, cone.Ng_cone(valid_lower), select (sph0_hit_first, cone.Ng_sphere0(valid_lower), cone.Ng_sphere1(valid_lower))); - const vfloat u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat(zero), vfloat(one))); + const vfloat t_first = select(valid_lower, t_lower, t_upper); + const vfloat u_first = select(valid_lower, u_front, u_back); /* invoke intersection filter for first hit */ - RoundLineIntersectorHitM hit(u_first,zero,dt+t_first,Ng_first); + RoundLineIntersectorHitM hit(u_first,zero,t_first,t_first*ray_dir_in + P + u_first*dP); const bool is_hit_first = epilog(valid_first, hit); - + /* check for possible second hits before potentially accepted hit */ - const vfloat t_second = t_cone_sphere_upper; - const vbool valid_second = valid_lower & valid_upper & (dt+t_cone_sphere_upper <= ray_tfar()); + const vbool valid_second = valid_lower & valid_upper & t_upper <= ray_tfar(); if (unlikely(none(valid_second))) return is_hit_first; - - /* invoke intersection filter for second hit */ - const vbool cone_hit_second = t_second == t_cone_lower | t_second == t_cone_upper; - const vbool sph0_hit_second = t_second == t_sph0_lower | t_second == t_sph0_upper; - const Vec3vf Ng_second = select(cone_hit_second, cone.Ng_cone(false), select (sph0_hit_second, cone.Ng_sphere0(false), cone.Ng_sphere1(false))); - const vfloat u_second = select(cone_hit_second, cone.u_cone(false), select (sph0_hit_second, vfloat(zero), vfloat(one))); - hit = RoundLineIntersectorHitM(u_second,zero,dt+t_second,Ng_second); + /* invoke intersection filter for second hit */ + hit = RoundLineIntersectorHitM(u_back,zero,t_upper,t_upper*ray_dir_in + P + u_back*dP); const bool is_hit_second = epilog(valid_second, hit); - + return is_hit_first | is_hit_second; #else + // scale back + t_lower *= rd; + /* filter out hits that are not in tnear/tfar range */ - const vbool valid_lower = valid & ray_tnear <= dt+t_cone_sphere_lower & dt+t_cone_sphere_lower <= ray_tfar() & t_cone_sphere_lower != vfloat(pos_inf); - + const vbool valid_lower = validFront & ray_tnear <= t_lower & t_lower <= ray_tfar(); + /* check if there is a valid hit */ if (unlikely(none(valid_lower))) return false; - - /* construct first hit */ - const vbool cone_hit_first = t_cone_sphere_lower == t_cone_lower | t_cone_sphere_lower == t_cone_upper; - const vbool sph0_hit_first = t_cone_sphere_lower == t_sph0_lower | t_cone_sphere_lower == t_sph0_upper; - const Vec3vf Ng_first = select(cone_hit_first, cone.Ng_cone(valid_lower), select (sph0_hit_first, cone.Ng_sphere0(valid_lower), cone.Ng_sphere1(valid_lower))); - const vfloat u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat(zero), vfloat(one))); - /* invoke intersection filter for first hit */ - RoundLineIntersectorHitM hit(u_first,zero,dt+t_cone_sphere_lower,Ng_first); + /* construct first hit and invoke intersection filter for first hit */ + RoundLineIntersectorHitM hit(u_front,zero,t_lower,t_lower*ray_dir_in + P + u_front*dP); const bool is_hit_first = epilog(valid_lower, hit); - + return is_hit_first; #endif } - + } // end namespace __roundline_internal - + template struct RoundLinearCurveIntersector1 { @@ -678,19 +416,19 @@ namespace embree return __roundline_internal::intersectConeSphere(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray),v0,v1,vL,vR,epilog); } }; - + template struct RoundLinearCurveIntersectorK { typedef CurvePrecalculationsK Precalculations; - + struct ray_tfar { RayK& ray; size_t k; __forceinline ray_tfar(RayK& ray, size_t k) : ray(ray), k(k) {} __forceinline vfloat operator() () const { return ray.tfar[k]; }; }; - + template static __forceinline bool intersect(const vbool& valid_i, RayK& ray, size_t k,