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 | |
67 | namespace 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 | |