1// Copyright 2009-2021 Intel Corporation
2// SPDX-License-Identifier: Apache-2.0
3
4#pragma once
5
6#include "../common/ray.h"
7#include "curve_intersector_precalculations.h"
8
9
10/*
11
12 This file implements the intersection of a ray with a round linear
13 curve segment. We define the geometry of such a round linear curve
14 segment from point p0 with radius r0 to point p1 with radius r1
15 using the cone that touches spheres p0/r0 and p1/r1 tangentially
16 plus the sphere p1/r1. We denote the tangentially touching cone from
17 p0/r0 to p1/r1 with cone(p0,r0,p1,r1) and the cone plus the ending
18 sphere with cone_sphere(p0,r0,p1,r1).
19
20 For multiple connected round linear curve segments this construction
21 yield a proper shape when viewed from the outside. Using the
22 following CSG we can also handle the interior in most common cases:
23
24 round_linear_curve(pl,rl,p0,r0,p1,r1,pr,rr) =
25 cone_sphere(p0,r0,p1,r1) - cone(pl,rl,p0,r0) - cone(p1,r1,pr,rr)
26
27 Thus by subtracting the neighboring cone geometries, we cut away
28 parts of the center cone_sphere surface which lie inside the
29 combined curve. This approach works as long as geometry of the
30 current cone_sphere penetrates into direct neighbor segments only,
31 and not into segments further away.
32
33 To construct a cone that touches two spheres at p0 and p1 with r0
34 and r1, one has to increase the cone radius at r0 and r1 to obtain
35 larger radii w0 and w1, such that the infinite cone properly touches
36 the spheres. From the paper "Ray Tracing Generalized Tube
37 Primitives: Method and Applications"
38 (https://www.researchgate.net/publication/334378683_Ray_Tracing_Generalized_Tube_Primitives_Method_and_Applications)
39 one can derive the following equations for these increased
40 radii:
41
42 sr = 1.0f / sqrt(1-sqr(dr)/sqr(p1-p0))
43 w0 = sr*r0
44 w1 = sr*r1
45
46 Further, we want the cone to start where it touches the sphere at p0
47 and to end where it touches sphere at p1. Therefore, we need to
48 construct clipping locations y0 and y1 for the start and end of the
49 cone. These start and end clipping location of the cone can get
50 calculated as:
51
52 Y0 = - r0 * (r1-r0) / length(p1-p0)
53 Y1 = length(p1-p0) - r1 * (r1-r0) / length(p1-p0)
54
55 Where the cone starts a distance Y0 and ends a distance Y1 away of
56 point p0 along the cone center. The distance between Y1-Y0 can get
57 calculated as:
58
59 dY = length(p1-p0) - (r1-r0)^2 / length(p1-p0)
60
61 In the code below, Y will always be scaled by length(p1-p0) to
62 obtain y and you will find the terms r0*(r1-r0) and
63 (p1-p0)^2-(r1-r0)^2.
64
65 */
66
67namespace embree
68{
69 namespace isa
70 {
71 template<int M>
72 struct RoundLineIntersectorHitM
73 {
74 __forceinline RoundLineIntersectorHitM() {}
75
76 __forceinline RoundLineIntersectorHitM(const vfloat<M>& u, const vfloat<M>& v, const vfloat<M>& t, const Vec3vf<M>& Ng)
77 : vu(u), vv(v), vt(t), vNg(Ng) {}
78
79 __forceinline void finalize() {}
80
81 __forceinline Vec2f uv (const size_t i) const { return Vec2f(vu[i],vv[i]); }
82 __forceinline float t (const size_t i) const { return vt[i]; }
83 __forceinline Vec3fa Ng(const size_t i) const { return Vec3fa(vNg.x[i],vNg.y[i],vNg.z[i]); }
84
85 __forceinline Vec2vf<M> uv() const { return Vec2vf<M>(vu,vv); }
86 __forceinline vfloat<M> t () const { return vt; }
87 __forceinline Vec3vf<M> Ng() const { return vNg; }
88
89 public:
90 vfloat<M> vu;
91 vfloat<M> vv;
92 vfloat<M> vt;
93 Vec3vf<M> vNg;
94 };
95
96 namespace __roundline_internal
97 {
98 template<int M>
99 struct ConeGeometry
100 {
101 ConeGeometry (const Vec4vf<M>& a, const Vec4vf<M>& b)
102 : 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) {}
103
104 /*
105
106 This function tests if a point is accepted by first cone
107 clipping plane.
108
109 First, we need to project the point onto the line p0->p1:
110
111 Y = (p-p0)*(p1-p0)/length(p1-p0)
112
113 This value y is the distance to the projection point from
114 p0. The clip distances are calculated as:
115
116 Y0 = - r0 * (r1-r0) / length(p1-p0)
117 Y1 = length(p1-p0) - r1 * (r1-r0) / length(p1-p0)
118
119 Thus to test if the point p is accepted by the first
120 clipping plane we need to test Y > Y0 and to test if it
121 is accepted by the second clipping plane we need to test
122 Y < Y1.
123
124 By multiplying the calculations with length(p1-p0) these
125 calculation can get simplied to:
126
127 y = (p-p0)*(p1-p0)
128 y0 = - r0 * (r1-r0)
129 y1 = (p1-p0)^2 - r1 * (r1-r0)
130
131 and the test y > y0 and y < y1.
132
133 */
134
135 __forceinline vbool<M> isClippedByPlane (const vbool<M>& valid_i, const Vec3vf<M>& p) const
136 {
137 const Vec3vf<M> p0p = p - p0;
138 const vfloat<M> y = dot(p0p,dP);
139 const vfloat<M> cap0 = -r0dr;
140 const vbool<M> inside_cone = y > cap0;
141 return valid_i & (p0.x != vfloat<M>(inf)) & (p1.x != vfloat<M>(inf)) & inside_cone;
142 }
143
144 /*
145
146 This function tests whether a point lies inside the capped cone
147 tangential to its ending spheres.
148
149 Therefore one has to check if the point is inside the
150 region defined by the cone clipping planes, which is
151 performed similar as in the previous function.
152
153 To perform the inside cone test we need to project the
154 point onto the line p0->p1:
155
156 dP = p1-p0
157 Y = (p-p0)*dP/length(dP)
158
159 This value Y is the distance to the projection point from
160 p0. To obtain a parameter value u going from 0 to 1 along
161 the line p0->p1 we calculate:
162
163 U = Y/length(dP)
164
165 The radii to use at points p0 and p1 are:
166
167 w0 = sr * r0
168 w1 = sr * r1
169 dw = w1-w0
170
171 Using these radii and u one can directly test if the point
172 lies inside the cone using the formula dP*dP < wy*wy with:
173
174 wy = w0 + u*dw
175 py = p0 + u*dP - p
176
177 By multiplying the calculations with length(p1-p0) and
178 inserting the definition of w can obtain simpler equations:
179
180 y = (p-p0)*dP
181 ry = r0 + y/dP^2 * dr
182 wy = sr*ry
183 py = p0 + y/dP^2*dP - p
184 y0 = - r0 * dr
185 y1 = dP^2 - r1 * dr
186
187 Thus for the in-cone test we get:
188
189 py^2 < wy^2
190 <=> py^2 < sr^2 * ry^2
191 <=> py^2 * ( dP^2 - dr^2 ) < dP^2 * ry^2
192
193 This can further get simplified to:
194
195 (p0-p)^2 * (dP^2 - dr^2) - y^2 < dP^2 * r0^2 + 2.0f*r0*dr*y;
196
197 */
198
199 __forceinline vbool<M> isInsideCappedCone (const vbool<M>& valid_i, const Vec3vf<M>& p) const
200 {
201 const Vec3vf<M> p0p = p - p0;
202 const vfloat<M> y = dot(p0p,dP);
203 const vfloat<M> cap0 = -r0dr+vfloat<M>(ulp);
204 const vfloat<M> cap1 = -r1*dr + dPdP;
205
206 vbool<M> inside_cone = valid_i & (p0.x != vfloat<M>(inf)) & (p1.x != vfloat<M>(inf));
207 inside_cone &= y > cap0; // start clipping plane
208 inside_cone &= y < cap1; // end clipping plane
209 inside_cone &= sqr(p0p)*g - sqr(y) < dPdP * sqr_r0 + 2.0f*r0dr*y; // in cone test
210 return inside_cone;
211 }
212
213 protected:
214 Vec3vf<M> p0;
215 Vec3vf<M> p1;
216 Vec3vf<M> dP;
217 vfloat<M> dPdP;
218 vfloat<M> r0;
219 vfloat<M> sqr_r0;
220 vfloat<M> r1;
221 vfloat<M> dr;
222 vfloat<M> drdr;
223 vfloat<M> r0dr;
224 vfloat<M> g;
225 };
226
227 template<int M>
228 struct ConeGeometryIntersector : public ConeGeometry<M>
229 {
230 using ConeGeometry<M>::p0;
231 using ConeGeometry<M>::p1;
232 using ConeGeometry<M>::dP;
233 using ConeGeometry<M>::dPdP;
234 using ConeGeometry<M>::r0;
235 using ConeGeometry<M>::sqr_r0;
236 using ConeGeometry<M>::r1;
237 using ConeGeometry<M>::dr;
238 using ConeGeometry<M>::r0dr;
239 using ConeGeometry<M>::g;
240
241 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)
242 : 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) {}
243
244 /*
245
246 This function intersects a ray with a cone that touches a
247 start sphere p0/r0 and end sphere p1/r1.
248
249 To find this ray/cone intersections one could just
250 calculate radii w0 and w1 as described above and use a
251 standard ray/cone intersection routine with these
252 radii. However, it turns out that calculations can get
253 simplified when deriving a specialized ray/cone
254 intersection for this special case. We perform
255 calculations relative to the cone origin p0 and define:
256
257 O = ray_org - p0
258 dO = ray_dir
259 dP = p1-p0
260 dr = r1-r0
261 dw = w1-w0
262
263 For some t we can compute the potential hit point h = O + t*dO and
264 project it onto the cone vector dP to obtain u = (h*dP)/(dP*dP). In
265 case of an intersection, the squared distance from the hit point
266 projected onto the cone center line to the hit point should be equal
267 to the squared cone radius at u:
268
269 (u*dP - h)^2 = (w0 + u*dw)^2
270
271 Inserting the definition of h, u, w0, and dw into this formula, then
272 factoring out all terms, and sorting by t^2, t^1, and t^0 terms
273 yields a quadratic equation to solve.
274
275 Inserting u:
276 ( (h*dP)*dP/dP^2 - h )^2 = ( w0 + (h*dP)*dw/dP^2 )^2
277
278 Multiplying by dP^4:
279 ( (h*dP)*dP - h*dP^2 )^2 = ( w0*dP^2 + (h*dP)*dw )^2
280
281 Inserting w0 and dw:
282 ( (h*dP)*dP - h*dP^2 )^2 = ( r0*dP^2 + (h*dP)*dr )^2 / (1-dr^2/dP^2)
283 ( (h*dP)*dP - h*dP^2 )^2 *(dP^2 - dr^2) = dP^2 * ( r0*dP^2 + (h*dP)*dr )^2
284
285 Now one can insert the definition of h, factor out, and presort by t:
286 ( ((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
287 ( (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
288
289 Factoring out further and sorting by t^2, t^1 and t^0 yields:
290
291 0 = t^2 * [ ((dO*dP)*dP - dO-dP^2)^2 * (dP^2 - dr^2) - dP^2*(dO*dP)^2*dr^2 ]
292 + 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 ]
293 + t^0 * [ ( (O*dP)*dP - O*dP^2)^2 * (dP^2-dr^2) - dP^2*(r0*dP^2 + (O*dP)*dr)^2 ]
294
295 This can be simplified to:
296
297 0 = t^2 * [ (dP^2 - dr^2)*dO^2 - (dO*dP)^2 ]
298 + 2*t^1 * [ (dP^2 - dr^2)*(O*dO) - (dO*dP)*(O*dP + r0*dr) ]
299 + t^0 * [ (dP^2 - dr^2)*O^2 - (O*dP)^2 - r0^2*dP^2 - 2.0f*r0*dr*(O*dP) ]
300
301 Solving this quadratic equation yields the values for t at which the
302 ray intersects the cone.
303
304 */
305
306 __forceinline bool intersectCone(vbool<M>& valid, vfloat<M>& lower, vfloat<M>& upper)
307 {
308 /* return no hit by default */
309 lower = pos_inf;
310 upper = neg_inf;
311
312 /* compute quadratic equation A*t^2 + B*t + C = 0 */
313 const vfloat<M> OO = dot(O,O);
314 const vfloat<M> OdO = dot(dO,O);
315 const vfloat<M> A = g * dOdO - sqr(dOdP);
316 const vfloat<M> B = 2.0f * (g*OdO - dOdP*yp);
317 const vfloat<M> C = g*OO - sqr(OdP) - sqr_r0*dPdP - 2.0f*r0dr*OdP;
318
319 /* we miss the cone if determinant is smaller than zero */
320 const vfloat<M> D = B*B - 4.0f*A*C;
321 valid &= (D >= 0.0f & g > 0.0f); // if g <= 0 then the cone is inside a sphere end
322
323 /* When rays are parallel to the cone surface, then the
324 * ray may be inside or outside the cone. We just assume a
325 * miss in that case, which is fine as rays inside the
326 * cone would anyway hit the ending spheres in that
327 * case. */
328 valid &= abs(A) > min_rcp_input;
329 if (unlikely(none(valid))) {
330 return false;
331 }
332
333 /* compute distance to front and back hit */
334 const vfloat<M> Q = sqrt(D);
335 const vfloat<M> rcp_2A = rcp(2.0f*A);
336 t_cone_front = (-B-Q)*rcp_2A;
337 y_cone_front = yp + t_cone_front*dOdP;
338 lower = select( (y_cone_front > -(float)ulp) & (y_cone_front <= g) & (g > 0.0f), t_cone_front, vfloat<M>(pos_inf));
339#if !defined (EMBREE_BACKFACE_CULLING_CURVES)
340 t_cone_back = (-B+Q)*rcp_2A;
341 y_cone_back = yp + t_cone_back *dOdP;
342 upper = select( (y_cone_back > -(float)ulp) & (y_cone_back <= g) & (g > 0.0f), t_cone_back , vfloat<M>(neg_inf));
343#endif
344 return true;
345 }
346
347 /*
348 This function intersects the ray with the end sphere at
349 p1. We already clip away hits that are inside the
350 neighboring cone segment.
351
352 */
353
354 __forceinline void intersectEndSphere(vbool<M>& valid,
355 const ConeGeometry<M>& coneR,
356 vfloat<M>& lower, vfloat<M>& upper)
357 {
358 /* calculate front and back hit with end sphere */
359 const Vec3vf<M> O1 = org - p1;
360 const vfloat<M> O1dO = dot(O1,dO);
361 const vfloat<M> h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r1));
362 const vfloat<M> rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat<M>(neg_inf) );
363
364 /* clip away front hit if it is inside next cone segment */
365 t_sph1_front = (-O1dO - rhs1)*rcp_dOdO;
366 const Vec3vf<M> hit_front = org + t_sph1_front*dO;
367 vbool<M> valid_sph1_front = h2 >= 0.0f & yp + t_sph1_front*dOdP > g & !coneR.isClippedByPlane (valid, hit_front);
368 lower = select(valid_sph1_front, t_sph1_front, vfloat<M>(pos_inf));
369
370#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
371 /* clip away back hit if it is inside next cone segment */
372 t_sph1_back = (-O1dO + rhs1)*rcp_dOdO;
373 const Vec3vf<M> hit_back = org + t_sph1_back*dO;
374 vbool<M> valid_sph1_back = h2 >= 0.0f & yp + t_sph1_back*dOdP > g & !coneR.isClippedByPlane (valid, hit_back);
375 upper = select(valid_sph1_back, t_sph1_back, vfloat<M>(neg_inf));
376#else
377 upper = vfloat<M>(neg_inf);
378#endif
379 }
380
381 __forceinline void intersectBeginSphere(const vbool<M>& valid,
382 vfloat<M>& lower, vfloat<M>& upper)
383 {
384 /* calculate front and back hit with end sphere */
385 const Vec3vf<M> O1 = org - p0;
386 const vfloat<M> O1dO = dot(O1,dO);
387 const vfloat<M> h2 = sqr(O1dO) - dOdO*(sqr(O1) - sqr(r0));
388 const vfloat<M> rhs1 = select( h2 >= 0.0f, sqrt(h2), vfloat<M>(neg_inf) );
389
390 /* clip away front hit if it is inside next cone segment */
391 t_sph0_front = (-O1dO - rhs1)*rcp_dOdO;
392 vbool<M> valid_sph1_front = valid & h2 >= 0.0f & yp + t_sph0_front*dOdP < 0;
393 lower = select(valid_sph1_front, t_sph0_front, vfloat<M>(pos_inf));
394
395#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
396 /* clip away back hit if it is inside next cone segment */
397 t_sph0_back = (-O1dO + rhs1)*rcp_dOdO;
398 vbool<M> valid_sph1_back = valid & h2 >= 0.0f & yp + t_sph0_back*dOdP < 0;
399 upper = select(valid_sph1_back, t_sph0_back, vfloat<M>(neg_inf));
400#else
401 upper = vfloat<M>(neg_inf);
402#endif
403 }
404
405 /*
406
407 This function calculates the geometry normal of some cone hit.
408
409 For a given hit point h (relative to p0) with a cone
410 starting at p0 with radius w0 and ending at p1 with
411 radius w1 one normally calculates the geometry normal by
412 first calculating the parmetric u hit location along the
413 cone:
414
415 u = dot(h,dP)/dP^2
416
417 Using this value one can now directly calculate the
418 geometry normal by bending the connection vector (h-u*dP)
419 from hit to projected hit with some cone dependent value
420 dw/sqrt(dP^2) * normalize(dP):
421
422 Ng = normalize(h-u*dP) - dw/length(dP) * normalize(dP)
423
424 The length of the vector (h-u*dP) can also get calculated
425 by interpolating the radii as w0+u*dw which yields:
426
427 Ng = (h-u*dP)/(w0+u*dw) - dw/dP^2 * dP
428
429 Multiplying with (w0+u*dw) yield a scaled Ng':
430
431 Ng' = (h-u*dP) - (w0+u*dw)*dw/dP^2*dP
432
433 Inserting the definition of w0 and dw and refactoring
434 yield a further scaled Ng'':
435
436 Ng'' = (dP^2 - dr^2) (h-q) - (r0+u*dr)*dr*dP
437
438 Now inserting the definition of u gives and multiplying
439 with the denominator yields:
440
441 Ng''' = (dP^2-dr^2)*(dP^2*h-dot(h,dP)*dP) - (dP^2*r0+dot(h,dP)*dr)*dr*dP
442
443 Factoring out, cancelling terms, dividing by dP^2, and
444 factoring again yields finally:
445
446 Ng'''' = (dP^2-dr^2)*h - dP*(dot(h,dP) + r0*dr)
447
448 */
449
450 __forceinline Vec3vf<M> Ng_cone(const vbool<M>& front_hit) const
451 {
452#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
453 const vfloat<M> y = select(front_hit, y_cone_front, y_cone_back);
454 const vfloat<M> t = select(front_hit, t_cone_front, t_cone_back);
455 const Vec3vf<M> h = O + t*dO;
456 return g*h-dP*y;
457#else
458 const Vec3vf<M> h = O + t_cone_front*dO;
459 return g*h-dP*y_cone_front;
460#endif
461 }
462
463 /* compute geometry normal of sphere hit as the difference
464 * vector from hit point to sphere center */
465
466 __forceinline Vec3vf<M> Ng_sphere1(const vbool<M>& front_hit) const
467 {
468#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
469 const vfloat<M> t_sph1 = select(front_hit, t_sph1_front, t_sph1_back);
470 return org+t_sph1*dO-p1;
471#else
472 return org+t_sph1_front*dO-p1;
473#endif
474 }
475
476 __forceinline Vec3vf<M> Ng_sphere0(const vbool<M>& front_hit) const
477 {
478#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
479 const vfloat<M> t_sph0 = select(front_hit, t_sph0_front, t_sph0_back);
480 return org+t_sph0*dO-p0;
481#else
482 return org+t_sph0_front*dO-p0;
483#endif
484 }
485
486 /*
487 This function calculates the u coordinate of a
488 hit. Therefore we use the hit distance y (which is zero
489 at the first cone clipping plane) and divide by distance
490 g between the clipping planes.
491
492 */
493
494 __forceinline vfloat<M> u_cone(const vbool<M>& front_hit) const
495 {
496#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
497 const vfloat<M> y = select(front_hit, y_cone_front, y_cone_back);
498 return clamp(y*rcp(g));
499#else
500 return clamp(y_cone_front*rcp(g));
501#endif
502 }
503
504 private:
505 Vec3vf<M> org;
506 Vec3vf<M> O;
507 Vec3vf<M> dO;
508 vfloat<M> dOdO;
509 vfloat<M> rcp_dOdO;
510 vfloat<M> OdP;
511 vfloat<M> dOdP;
512
513 /* for ray/cone intersection */
514 private:
515 vfloat<M> yp;
516 vfloat<M> y_cone_front;
517 vfloat<M> t_cone_front;
518#if !defined (EMBREE_BACKFACE_CULLING_CURVES)
519 vfloat<M> y_cone_back;
520 vfloat<M> t_cone_back;
521#endif
522
523 /* for ray/sphere intersection */
524 private:
525 vfloat<M> t_sph1_front;
526 vfloat<M> t_sph0_front;
527#if !defined (EMBREE_BACKFACE_CULLING_CURVES)
528 vfloat<M> t_sph1_back;
529 vfloat<M> t_sph0_back;
530#endif
531 };
532
533
534 template<int M, typename Epilog, typename ray_tfar_func>
535 static __forceinline bool intersectConeSphere(const vbool<M>& valid_i,
536 const Vec3vf<M>& ray_org_in, const Vec3vf<M>& ray_dir,
537 const vfloat<M>& ray_tnear, const ray_tfar_func& ray_tfar,
538 const Vec4vf<M>& v0, const Vec4vf<M>& v1,
539 const Vec4vf<M>& vL, const Vec4vf<M>& vR,
540 const Epilog& epilog)
541 {
542 vbool<M> valid = valid_i;
543
544 /* move ray origin closer to make calculations numerically stable */
545 const vfloat<M> dOdO = sqr(ray_dir);
546 const vfloat<M> rcp_dOdO = rcp(dOdO);
547 const Vec3vf<M> center = vfloat<M>(0.5f)*(v0.xyz()+v1.xyz());
548 const vfloat<M> dt = dot(center-ray_org_in,ray_dir)*rcp_dOdO;
549 const Vec3vf<M> ray_org = ray_org_in + dt*ray_dir;
550
551 /* intersect with cone from v0 to v1 */
552 vfloat<M> t_cone_lower, t_cone_upper;
553 ConeGeometryIntersector<M> cone (ray_org, ray_dir, dOdO, rcp_dOdO, v0, v1);
554 vbool<M> validCone = valid;
555 cone.intersectCone(validCone, t_cone_lower, t_cone_upper);
556
557 valid &= (validCone | (cone.g <= 0.0f)); // if cone is entirely in sphere end - check sphere
558 if (unlikely(none(valid)))
559 return false;
560
561 /* cone hits inside the neighboring capped cones are inside the geometry and thus ignored */
562 const ConeGeometry<M> coneL (v0, vL);
563 const ConeGeometry<M> coneR (v1, vR);
564#if !defined(EMBREE_BACKFACE_CULLING_CURVES)
565 const Vec3vf<M> hit_lower = ray_org + t_cone_lower*ray_dir;
566 const Vec3vf<M> hit_upper = ray_org + t_cone_upper*ray_dir;
567 t_cone_lower = select (!coneL.isInsideCappedCone (validCone, hit_lower) & !coneR.isInsideCappedCone (validCone, hit_lower), t_cone_lower, vfloat<M>(pos_inf));
568 t_cone_upper = select (!coneL.isInsideCappedCone (validCone, hit_upper) & !coneR.isInsideCappedCone (validCone, hit_upper), t_cone_upper, vfloat<M>(neg_inf));
569#endif
570
571 /* intersect ending sphere */
572 vfloat<M> t_sph1_lower, t_sph1_upper;
573 vfloat<M> t_sph0_lower = vfloat<M>(pos_inf);
574 vfloat<M> t_sph0_upper = vfloat<M>(neg_inf);
575 cone.intersectEndSphere(valid, coneR, t_sph1_lower, t_sph1_upper);
576
577 const vbool<M> isBeginPoint = valid & (vL[0] == vfloat<M>(pos_inf));
578 if (unlikely(any(isBeginPoint))) {
579 cone.intersectBeginSphere (isBeginPoint, t_sph0_lower, t_sph0_upper);
580 }
581
582 /* CSG union of cone and end sphere */
583 vfloat<M> t_sph_lower = min(t_sph0_lower, t_sph1_lower);
584 vfloat<M> t_cone_sphere_lower = min(t_cone_lower, t_sph_lower);
585#if !defined (EMBREE_BACKFACE_CULLING_CURVES)
586 vfloat<M> t_sph_upper = max(t_sph0_upper, t_sph1_upper);
587 vfloat<M> t_cone_sphere_upper = max(t_cone_upper, t_sph_upper);
588
589 /* filter out hits that are not in tnear/tfar range */
590 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);
591 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);
592
593 /* check if there is a first hit */
594 const vbool<M> valid_first = valid_lower | valid_upper;
595 if (unlikely(none(valid_first)))
596 return false;
597
598 /* construct first hit */
599 const vfloat<M> t_first = select(valid_lower, t_cone_sphere_lower, t_cone_sphere_upper);
600 const vbool<M> cone_hit_first = t_first == t_cone_lower | t_first == t_cone_upper;
601 const vbool<M> sph0_hit_first = t_first == t_sph0_lower | t_first == t_sph0_upper;
602 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)));
603 const vfloat<M> u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat<M>(zero), vfloat<M>(one)));
604
605 /* invoke intersection filter for first hit */
606 RoundLineIntersectorHitM<M> hit(u_first,zero,dt+t_first,Ng_first);
607 const bool is_hit_first = epilog(valid_first, hit);
608
609 /* check for possible second hits before potentially accepted hit */
610 const vfloat<M> t_second = t_cone_sphere_upper;
611 const vbool<M> valid_second = valid_lower & valid_upper & (dt+t_cone_sphere_upper <= ray_tfar());
612 if (unlikely(none(valid_second)))
613 return is_hit_first;
614
615 /* invoke intersection filter for second hit */
616 const vbool<M> cone_hit_second = t_second == t_cone_lower | t_second == t_cone_upper;
617 const vbool<M> sph0_hit_second = t_second == t_sph0_lower | t_second == t_sph0_upper;
618 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)));
619 const vfloat<M> u_second = select(cone_hit_second, cone.u_cone(false), select (sph0_hit_second, vfloat<M>(zero), vfloat<M>(one)));
620
621 hit = RoundLineIntersectorHitM<M>(u_second,zero,dt+t_second,Ng_second);
622 const bool is_hit_second = epilog(valid_second, hit);
623
624 return is_hit_first | is_hit_second;
625#else
626 /* filter out hits that are not in tnear/tfar range */
627 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);
628
629 /* check if there is a valid hit */
630 if (unlikely(none(valid_lower)))
631 return false;
632
633 /* construct first hit */
634 const vbool<M> cone_hit_first = t_cone_sphere_lower == t_cone_lower | t_cone_sphere_lower == t_cone_upper;
635 const vbool<M> sph0_hit_first = t_cone_sphere_lower == t_sph0_lower | t_cone_sphere_lower == t_sph0_upper;
636 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)));
637 const vfloat<M> u_first = select(cone_hit_first, cone.u_cone(valid_lower), select (sph0_hit_first, vfloat<M>(zero), vfloat<M>(one)));
638
639 /* invoke intersection filter for first hit */
640 RoundLineIntersectorHitM<M> hit(u_first,zero,dt+t_cone_sphere_lower,Ng_first);
641 const bool is_hit_first = epilog(valid_lower, hit);
642
643 return is_hit_first;
644#endif
645 }
646
647 } // end namespace __roundline_internal
648
649 template<int M>
650 struct RoundLinearCurveIntersector1
651 {
652 typedef CurvePrecalculations1 Precalculations;
653
654 template<typename Ray>
655 struct ray_tfar {
656 Ray& ray;
657 __forceinline ray_tfar(Ray& ray) : ray(ray) {}
658 __forceinline vfloat<M> operator() () const { return ray.tfar; };
659 };
660
661 template<typename Ray, typename Epilog>
662 static __forceinline bool intersect(const vbool<M>& valid_i,
663 Ray& ray,
664 IntersectContext* context,
665 const LineSegments* geom,
666 const Precalculations& pre,
667 const Vec4vf<M>& v0i, const Vec4vf<M>& v1i,
668 const Vec4vf<M>& vLi, const Vec4vf<M>& vRi,
669 const Epilog& epilog)
670 {
671 const Vec3vf<M> ray_org(ray.org.x, ray.org.y, ray.org.z);
672 const Vec3vf<M> ray_dir(ray.dir.x, ray.dir.y, ray.dir.z);
673 const vfloat<M> ray_tnear(ray.tnear());
674 const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i);
675 const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i);
676 const Vec4vf<M> vL = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vLi);
677 const Vec4vf<M> vR = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vRi);
678 return __roundline_internal::intersectConeSphere<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar<Ray>(ray),v0,v1,vL,vR,epilog);
679 }
680 };
681
682 template<int M, int K>
683 struct RoundLinearCurveIntersectorK
684 {
685 typedef CurvePrecalculationsK<K> Precalculations;
686
687 struct ray_tfar {
688 RayK<K>& ray;
689 size_t k;
690 __forceinline ray_tfar(RayK<K>& ray, size_t k) : ray(ray), k(k) {}
691 __forceinline vfloat<M> operator() () const { return ray.tfar[k]; };
692 };
693
694 template<typename Epilog>
695 static __forceinline bool intersect(const vbool<M>& valid_i,
696 RayK<K>& ray, size_t k,
697 IntersectContext* context,
698 const LineSegments* geom,
699 const Precalculations& pre,
700 const Vec4vf<M>& v0i, const Vec4vf<M>& v1i,
701 const Vec4vf<M>& vLi, const Vec4vf<M>& vRi,
702 const Epilog& epilog)
703 {
704 const Vec3vf<M> ray_org(ray.org.x[k], ray.org.y[k], ray.org.z[k]);
705 const Vec3vf<M> ray_dir(ray.dir.x[k], ray.dir.y[k], ray.dir.z[k]);
706 const vfloat<M> ray_tnear = ray.tnear()[k];
707 const Vec4vf<M> v0 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v0i);
708 const Vec4vf<M> v1 = enlargeRadiusToMinWidth<M>(context,geom,ray_org,v1i);
709 const Vec4vf<M> vL = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vLi);
710 const Vec4vf<M> vR = enlargeRadiusToMinWidth<M>(context,geom,ray_org,vRi);
711 return __roundline_internal::intersectConeSphere<M>(valid_i,ray_org,ray_dir,ray_tnear,ray_tfar(ray,k),v0,v1,vL,vR,epilog);
712 }
713 };
714 }
715}
716