// Copyright 2009-2021 Intel Corporation // SPDX-License-Identifier: Apache-2.0 #pragma once #include "../common/ray.h" #include "curve_intersector_precalculations.h" /* 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 using the cone that touches spheres p0/r0 and p1/r1 tangentially plus the sphere p1/r1. We denote the tangentially touching cone from 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). 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 interiour in most common cases: round_linear_curve(pl,rl,p0,r0,p1,r1,pr,rr) = cone_sphere(p0,r0,p1,r1) - cone(pl,rl,p0,r0) - cone(p1,r1,pr,rr) Thus by subtracting the neighboring cone geometries, we cut away parts of the center cone_sphere surface which lie inside the 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 the spheres. From the paper "Ray Tracing Generalized Tube Primitives: Method and Applications" (https://www.researchgate.net/publication/334378683_Ray_Tracing_Generalized_Tube_Primitives_Method_and_Applications) one can derive the following equations for these increased radii: sr = 1.0f / sqrt(1-sqr(dr)/sqr(p1-p0)) w0 = sr*r0 w1 = sr*r1 Further, we want the cone to start where it touches the sphere at p0 and to end where it touches sphere at p1. Therefore, we need to construct clipping locations y0 and y1 for the start and end of the cone. These start and end clipping location of the cone can get calculated as: Y0 = - r0 * (r1-r0) / length(p1-p0) Y1 = length(p1-p0) - r1 * (r1-r0) / length(p1-p0) Where the cone starts a distance Y0 and ends a distance Y1 away of point p0 along the cone center. The distance between Y1-Y0 can get calculated as: dY = length(p1-p0) - (r1-r0)^2 / length(p1-p0) In the code below, Y will always be scaled by length(p1-p0) to obtain y and you will find the terms r0*(r1-r0) and (p1-p0)^2-(r1-r0)^2. */ namespace embree { namespace isa { template<int M> struct RoundLineIntersectorHitM { __forceinline RoundLineIntersectorHitM() {} __forceinline RoundLineIntersectorHitM(const vfloat<M>& u, const vfloat<M>& v, const vfloat<M>& t, const Vec3vf<M>& Ng) : vu(u), vv(v), vt(t), vNg(Ng) {} __forceinline void finalize() {} __forceinline Vec2f uv (const size_t i) const { return Vec2f(vu[i],vv[i]); } __forceinline float t (const size_t i) const { return vt[i]; } __forceinline Vec3fa Ng(const size_t i) const { return Vec3fa(vNg.x[i],vNg.y[i],vNg.z[i]); } __forceinline Vec2vf<M> uv() const { return Vec2vf<M>(vu,vv); } __forceinline vfloat<M> t () const { return vt; } __forceinline Vec3vf<M> Ng() const { return vNg; } public: vfloat<M> vu; vfloat<M> vv; vfloat<M> vt; Vec3vf<M> vNg; }; namespace __roundline_internal { template<int M> struct ConeGeometry { ConeGeometry (const Vec4vf<M>& a, const Vec4vf<M>& 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<M> isClippedByPlane (const vbool<M>& valid_i, const Vec3vf<M>& p) const { const Vec3vf<M> p0p = p - p0; const vfloat<M> y = dot(p0p,dP); const vfloat<M> cap0 = -r0dr; const vbool<M> inside_cone = y > cap0; return valid_i & (p0.x != vfloat<M>(inf)) & (p1.x != vfloat<M>(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 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; */ __forceinline vbool<M> isInsideCappedCone (const vbool<M>& valid_i, const Vec3vf<M>& p) const { const Vec3vf<M> p0p = p - p0; const vfloat<M> y = dot(p0p,dP); const vfloat<M> cap0 = -r0dr+vfloat<M>(ulp); const vfloat<M> cap1 = -r1*dr + dPdP; vbool<M> inside_cone = valid_i & (p0.x != vfloat<M>(inf)) & (p1.x != vfloat<M>(inf)); inside_cone &= y > cap0; // start 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<M> p0; Vec3vf<M> p1; Vec3vf<M> dP; vfloat<M> dPdP; vfloat<M> r0; vfloat<M> sqr_r0; vfloat<M> r1; vfloat<M> dr; vfloat<M> drdr; vfloat<M> r0dr; vfloat<M> g; }; template<int M> struct ConeGeometryIntersector : public ConeGeometry<M> { using ConeGeometry<M>::p0; using ConeGeometry<M>::p1; using ConeGeometry<M>::dP; using ConeGeometry<M>::dPdP; using ConeGeometry<M>::r0; using ConeGeometry<M>::sqr_r0; using ConeGeometry<M>::r1; using ConeGeometry<M>::dr; using ConeGeometry<M>::r0dr; using ConeGeometry<M>::g; ConeGeometryIntersector (const Vec3vf<M>& ray_org, const Vec3vf<M>& ray_dir, const vfloat<M>& dOdO, const vfloat<M>& rcp_dOdO, const Vec4vf<M>& a, const Vec4vf<M>& b) : ConeGeometry<M>(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) {} /* 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. */ __forceinline bool intersectCone(vbool<M>& valid, vfloat<M>& lower, vfloat<M>& upper) { /* return no hit by default */ lower = pos_inf; upper = neg_inf; /* compute quadratic equation A*t^2 + B*t + C = 0 */ const vfloat<M> OO = dot(O,O); const vfloat<M> OdO = dot(dO,O); const vfloat<M> A = g * dOdO - sqr(dOdP); const vfloat<M> B = 2.0f * (g*OdO - dOdP*yp); const vfloat<M> C = g*OO - sqr(OdP) - sqr_r0*dPdP - 2.0f*r0dr*OdP; /* we miss the cone if determinant is smaller than zero */ const vfloat<M> 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))) { return false; } /* compute distance to front and back hit */ const vfloat<M> Q = sqrt(D); const vfloat<M> 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<M>(pos_inf)); #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<M>(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<M>& valid, const ConeGeometry<M>& coneR, vfloat<M>& lower, vfloat<M>& upper) { /* calculate front and back hit with end sphere */ const Vec3vf<M> O1 = org - p1; const vfloat<M> O1dO = dot(O1,dO); const vfloat<M> h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r1)); const vfloat<M> rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat<M>(neg_inf) ); /* clip away front hit if it is inside next cone segment */ t_sph1_front = (-O1dO - rhs1)*rcp_dOdO; const Vec3vf<M> hit_front = org + t_sph1_front*dO; vbool<M> 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<M>(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<M> hit_back = org + t_sph1_back*dO; vbool<M> 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<M>(neg_inf)); #else upper = vfloat<M>(neg_inf); #endif } __forceinline void intersectBeginSphere(const vbool<M>& valid, vfloat<M>& lower, vfloat<M>& upper) { /* calculate front and back hit with end sphere */ const Vec3vf<M> O1 = org - p0; const vfloat<M> O1dO = dot(O1,dO); const vfloat<M> h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r0)); const vfloat<M> rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat<M>(neg_inf) ); /* clip away front hit if it is inside next cone segment */ t_sph0_front = (-O1dO - rhs1)*rcp_dOdO; vbool<M> valid_sph1_front = valid & h2 >= 0.0f & yp + t_sph0_front*dOdP < 0; lower = select(valid_sph1_front, t_sph0_front, vfloat<M>(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<M> valid_sph1_back = valid & h2 >= 0.0f & yp + t_sph0_back*dOdP < 0; upper = select(valid_sph1_back, t_sph0_back, vfloat<M>(neg_inf)); #else upper = vfloat<M>(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 furhter 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) */ __forceinline Vec3vf<M> Ng_cone(const vbool<M>& front_hit) const { #if !defined(EMBREE_BACKFACE_CULLING_CURVES) const vfloat<M> y = select(front_hit, y_cone_front, y_cone_back); const vfloat<M> t = select(front_hit, t_cone_front, t_cone_back); const Vec3vf<M> h = O + t*dO; return g*h-dP*y; #else const Vec3vf<M> 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<M> Ng_sphere1(const vbool<M>& front_hit) const { #if !defined(EMBREE_BACKFACE_CULLING_CURVES) const vfloat<M> 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 } __forceinline Vec3vf<M> Ng_sphere0(const vbool<M>& front_hit) const { #if !defined(EMBREE_BACKFACE_CULLING_CURVES) const vfloat<M> 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. */ __forceinline vfloat<M> u_cone(const vbool<M>& front_hit) const { #if !defined(EMBREE_BACKFACE_CULLING_CURVES) const vfloat<M> y = select(front_hit, y_cone_front, y_cone_back); return clamp(y*rcp(g)); #else return clamp(y_cone_front*rcp(g)); #endif } private: Vec3vf<M> org; Vec3vf<M> O; Vec3vf<M> dO; vfloat<M> dOdO; vfloat<M> rcp_dOdO; vfloat<M> OdP; vfloat<M> dOdP; /* for ray/cone intersection */ private: vfloat<M> yp; vfloat<M> y_cone_front; vfloat<M> t_cone_front; #if !defined (EMBREE_BACKFACE_CULLING_CURVES) vfloat<M> y_cone_back; vfloat<M> t_cone_back; #endif /* for ray/sphere intersection */ private: vfloat<M> t_sph1_front; vfloat<M> t_sph0_front; #if !defined (EMBREE_BACKFACE_CULLING_CURVES) vfloat<M> t_sph1_back; vfloat<M> t_sph0_back; #endif }; template<int M, typename Epilog, typename ray_tfar_func> static __forceinline bool intersectConeSphere(const vbool<M>& valid_i, const Vec3vf<M>& ray_org_in, const Vec3vf<M>& ray_dir, const vfloat<M>& ray_tnear, const ray_tfar_func& ray_tfar, const Vec4vf<M>& v0, const Vec4vf<M>& v1, const Vec4vf<M>& vL, const Vec4vf<M>& vR, const Epilog& epilog) { vbool<M> valid = valid_i; /* move ray origin closer to make calculations numerically stable */ const vfloat<M> dOdO = sqr(ray_dir); const vfloat<M> rcp_dOdO = rcp(dOdO); const Vec3vf<M> center = vfloat<M>(0.5f)*(v0.xyz()+v1.xyz()); const vfloat<M> dt = dot(center-ray_org_in,ray_dir)*rcp_dOdO; const Vec3vf<M> ray_org = ray_org_in + dt*ray_dir; /* intersect with cone from v0 to v1 */ vfloat<M> t_cone_lower, t_cone_upper; ConeGeometryIntersector<M> cone (ray_org, ray_dir, dOdO, rcp_dOdO, v0, v1); vbool<M> 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))) return false; /* cone hits inside the neighboring capped cones are inside the geometry and thus ignored */ const ConeGeometry<M> coneL (v0, vL); const ConeGeometry<M> coneR (v1, vR); #if !defined(EMBREE_BACKFACE_CULLING_CURVES) const Vec3vf<M> hit_lower = ray_org + t_cone_lower*ray_dir; const Vec3vf<M> 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<M>(pos_inf)); t_cone_upper = select (!coneL.isInsideCappedCone (validCone, hit_upper) & !coneR.isInsideCappedCone (validCone, hit_upper), t_cone_upper, vfloat<M>(neg_inf)); #endif /* intersect ending sphere */ vfloat<M> t_sph1_lower, t_sph1_upper; vfloat<M> t_sph0_lower = vfloat<M>(pos_inf); vfloat<M> t_sph0_upper = vfloat<M>(neg_inf); cone.intersectEndSphere(valid, coneR, t_sph1_lower, t_sph1_upper); const vbool<M> isBeginPoint = valid & (vL[0] == vfloat<M>(pos_inf)); if (unlikely(any(isBeginPoint))) { cone.intersectBeginSphere (isBeginPoint, t_sph0_lower, t_sph0_upper); } /* CSG union of cone and end sphere */ vfloat<M> t_sph_lower = min(t_sph0_lower, t_sph1_lower); vfloat<M> t_cone_sphere_lower = min(t_cone_lower, t_sph_lower); #if !defined (EMBREE_BACKFACE_CULLING_CURVES) vfloat<M> t_sph_upper = max(t_sph0_upper, t_sph1_upper); vfloat<M> t_cone_sphere_upper = max(t_cone_upper, t_sph_upper); /* filter out hits that are not in tnear/tfar range */ const vbool<M> valid_lower = valid & ray_tnear <= dt+t_cone_sphere_lower & dt+t_cone_sphere_lower <= ray_tfar() & t_cone_sphere_lower != vfloat<M>(pos_inf); const vbool<M> valid_upper = valid & ray_tnear <= dt+t_cone_sphere_upper & dt+t_cone_sphere_upper <= ray_tfar() & t_cone_sphere_upper != vfloat<M>(neg_inf); /* check if there is a first hit */ const vbool<M> valid_first = valid_lower | valid_upper; if (unlikely(none(valid_first))) return false; /* construct first hit */ const vfloat<M> t_first = select(valid_lower, t_cone_sphere_lower, t_cone_sphere_upper); const vbool<M> cone_hit_first = t_first == t_cone_lower | t_first == t_cone_upper; const vbool<M> sph0_hit_first = t_first == t_sph0_lower | t_first == t_sph0_upper; const Vec3vf<M> 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<M> u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat<M>(zero), vfloat<M>(one))); /* invoke intersection filter for first hit */ RoundLineIntersectorHitM<M> hit(u_first,zero,dt+t_first,Ng_first); const bool is_hit_first = epilog(valid_first, hit); /* check for possible second hits before potentially accepted hit */ const vfloat<M> t_second = t_cone_sphere_upper; const vbool<M> valid_second = valid_lower & valid_upper & (dt+t_cone_sphere_upper <= ray_tfar()); if (unlikely(none(valid_second))) return is_hit_first; /* invoke intersection filter for second hit */ const vbool<M> cone_hit_second = t_second == t_cone_lower | t_second == t_cone_upper; const vbool<M> sph0_hit_second = t_second == t_sph0_lower | t_second == t_sph0_upper; const Vec3vf<M> Ng_second = select(cone_hit_second, cone.Ng_cone(false), select (sph0_hit_second, cone.Ng_sphere0(false), cone.Ng_sphere1(false))); const vfloat<M> u_second = select(cone_hit_second, cone.u_cone(false), select (sph0_hit_second, vfloat<M>(zero), vfloat<M>(one))); hit = RoundLineIntersectorHitM<M>(u_second,zero,dt+t_second,Ng_second); const bool is_hit_second = epilog(valid_second, hit); return is_hit_first | is_hit_second; #else /* filter out hits that are not in tnear/tfar range */ const vbool<M> valid_lower = valid & ray_tnear <= dt+t_cone_sphere_lower & dt+t_cone_sphere_lower <= ray_tfar() & t_cone_sphere_lower != vfloat<M>(pos_inf); /* check if there is a valid hit */ if (unlikely(none(valid_lower))) return false; /* construct first hit */ const vbool<M> cone_hit_first = t_cone_sphere_lower == t_cone_lower | t_cone_sphere_lower == t_cone_upper; const vbool<M> sph0_hit_first = t_cone_sphere_lower == t_sph0_lower | t_cone_sphere_lower == t_sph0_upper; const Vec3vf<M> 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<M> u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat<M>(zero), vfloat<M>(one))); /* invoke intersection filter for first hit */ RoundLineIntersectorHitM<M> hit(u_first,zero,dt+t_cone_sphere_lower,Ng_first); const bool is_hit_first = epilog(valid_lower, hit); return is_hit_first; #endif } } // end namespace __roundline_internal template<int M> struct RoundLinearCurveIntersector1 { typedef CurvePrecalculations1 Precalculations; template<typename Ray> struct ray_tfar { Ray& ray; __forceinline ray_tfar(Ray& ray) : ray(ray) {} __forceinline vfloat<M> operator() () const { return ray.tfar; }; }; template<typename Ray, typename Epilog> static __forceinline bool intersect(const vbool<M>& valid_i, Ray& ray, IntersectContext* context, const LineSegments* geom, const Precalculations& pre, const Vec4vf<M>& v0i, const Vec4vf<M>& v1i, const Vec4vf<M>& vLi, const Vec4vf<M>& vRi, const Epilog& epilog) { const Vec3vf<M> ray_org(ray.org.x, ray.org.y, ray.org.z); const Vec3vf<M> ray_dir(ray.dir.x, ray.dir.y, ray.dir.z); const vfloat<M> ray_tnear(ray.tnear()); const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i); const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i); const Vec4vf<M> vL = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vLi); const Vec4vf<M> vR = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vRi); return __roundline_internal::intersectConeSphere<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar<Ray>(ray),v0,v1,vL,vR,epilog); } }; template<int M, int K> struct RoundLinearCurveIntersectorK { typedef CurvePrecalculationsK<K> Precalculations; struct ray_tfar { RayK<K>& ray; size_t k; __forceinline ray_tfar(RayK<K>& ray, size_t k) : ray(ray), k(k) {} __forceinline vfloat<M> operator() () const { return ray.tfar[k]; }; }; template<typename Epilog> static __forceinline bool intersect(const vbool<M>& valid_i, RayK<K>& ray, size_t k, IntersectContext* context, const LineSegments* geom, const Precalculations& pre, const Vec4vf<M>& v0i, const Vec4vf<M>& v1i, const Vec4vf<M>& vLi, const Vec4vf<M>& vRi, const Epilog& epilog) { const Vec3vf<M> ray_org(ray.org.x[k], ray.org.y[k], ray.org.z[k]); const Vec3vf<M> ray_dir(ray.dir.x[k], ray.dir.y[k], ray.dir.z[k]); const vfloat<M> ray_tnear = ray.tnear()[k]; const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i); const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i); const Vec4vf<M> vL = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vLi); const Vec4vf<M> vR = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vRi); return __roundline_internal::intersectConeSphere<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray,k),v0,v1,vL,vR,epilog); } }; } }