godot/thirdparty/embree/kernels/geometry/roundline_intersector.h
2022-11-25 13:09:04 +01:00

716 lines
31 KiB
C++

// 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 interior 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 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)
*/
__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);
}
};
}
}