711 lines
30 KiB
C++
711 lines
30 KiB
C++
// Copyright 2009-2020 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]); }
|
|
|
|
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;
|
|
|
|
struct ray_tfar {
|
|
Ray& ray;
|
|
__forceinline ray_tfar(Ray& ray) : ray(ray) {}
|
|
__forceinline vfloat<M> operator() () const { return ray.tfar; };
|
|
};
|
|
|
|
template<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(context,geom,ray_org,v0i);
|
|
const Vec4vf<M> v1 = enlargeRadiusToMinWidth(context,geom,ray_org,v1i);
|
|
const Vec4vf<M> vL = enlargeRadiusToMinWidth(context,geom,ray_org,vLi);
|
|
const Vec4vf<M> vR = enlargeRadiusToMinWidth(context,geom,ray_org,vRi);
|
|
return __roundline_internal::intersectConeSphere(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(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(context,geom,ray_org,v0i);
|
|
const Vec4vf<M> v1 = enlargeRadiusToMinWidth(context,geom,ray_org,v1i);
|
|
const Vec4vf<M> vL = enlargeRadiusToMinWidth(context,geom,ray_org,vLi);
|
|
const Vec4vf<M> vR = enlargeRadiusToMinWidth(context,geom,ray_org,vRi);
|
|
return __roundline_internal::intersectConeSphere(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray,k),v0,v1,vL,vR,epilog);
|
|
}
|
|
};
|
|
}
|
|
}
|