1 /*
2 * Copyright 2011-2017 Blender Foundation
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16
17 #ifndef __UTIL_MATH_INTERSECT_H__
18 #define __UTIL_MATH_INTERSECT_H__
19
20 CCL_NAMESPACE_BEGIN
21
22 /* Ray Intersection */
23
ray_sphere_intersect(float3 ray_P,float3 ray_D,float ray_t,float3 sphere_P,float sphere_radius,float3 * isect_P,float * isect_t)24 ccl_device bool ray_sphere_intersect(float3 ray_P,
25 float3 ray_D,
26 float ray_t,
27 float3 sphere_P,
28 float sphere_radius,
29 float3 *isect_P,
30 float *isect_t)
31 {
32 const float3 d = sphere_P - ray_P;
33 const float radiussq = sphere_radius * sphere_radius;
34 const float tsq = dot(d, d);
35
36 if (tsq > radiussq) {
37 /* Ray origin outside sphere. */
38 const float tp = dot(d, ray_D);
39 if (tp < 0.0f) {
40 /* Ray points away from sphere. */
41 return false;
42 }
43 const float dsq = tsq - tp * tp; /* pythagoras */
44 if (dsq > radiussq) {
45 /* Closest point on ray outside sphere. */
46 return false;
47 }
48 const float t = tp - sqrtf(radiussq - dsq); /* pythagoras */
49 if (t < ray_t) {
50 *isect_t = t;
51 *isect_P = ray_P + ray_D * t;
52 return true;
53 }
54 }
55 return false;
56 }
57
ray_aligned_disk_intersect(float3 ray_P,float3 ray_D,float ray_t,float3 disk_P,float disk_radius,float3 * isect_P,float * isect_t)58 ccl_device bool ray_aligned_disk_intersect(float3 ray_P,
59 float3 ray_D,
60 float ray_t,
61 float3 disk_P,
62 float disk_radius,
63 float3 *isect_P,
64 float *isect_t)
65 {
66 /* Aligned disk normal. */
67 float disk_t;
68 const float3 disk_N = normalize_len(ray_P - disk_P, &disk_t);
69 const float div = dot(ray_D, disk_N);
70 if (UNLIKELY(div == 0.0f)) {
71 return false;
72 }
73 /* Compute t to intersection point. */
74 const float t = -disk_t / div;
75 if (t < 0.0f || t > ray_t) {
76 return false;
77 }
78 /* Test if within radius. */
79 float3 P = ray_P + ray_D * t;
80 if (len_squared(P - disk_P) > disk_radius * disk_radius) {
81 return false;
82 }
83 *isect_P = P;
84 *isect_t = t;
85 return true;
86 }
87
ray_triangle_intersect(float3 ray_P,float3 ray_dir,float ray_t,const ssef * ssef_verts,float * isect_u,float * isect_v,float * isect_t)88 ccl_device_forceinline bool ray_triangle_intersect(float3 ray_P,
89 float3 ray_dir,
90 float ray_t,
91 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
92 const ssef *ssef_verts,
93 #else
94 const float3 tri_a,
95 const float3 tri_b,
96 const float3 tri_c,
97 #endif
98 float *isect_u,
99 float *isect_v,
100 float *isect_t)
101 {
102 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
103 typedef ssef float3;
104 const float3 tri_a(ssef_verts[0]);
105 const float3 tri_b(ssef_verts[1]);
106 const float3 tri_c(ssef_verts[2]);
107 const float3 P(ray_P);
108 const float3 dir(ray_dir);
109 #else
110 # define dot3(a, b) dot(a, b)
111 const float3 P = ray_P;
112 const float3 dir = ray_dir;
113 #endif
114
115 /* Calculate vertices relative to ray origin. */
116 const float3 v0 = tri_c - P;
117 const float3 v1 = tri_a - P;
118 const float3 v2 = tri_b - P;
119
120 /* Calculate triangle edges. */
121 const float3 e0 = v2 - v0;
122 const float3 e1 = v0 - v1;
123 const float3 e2 = v1 - v2;
124
125 /* Perform edge tests. */
126 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
127 const float3 crossU = cross(v2 + v0, e0);
128 const float3 crossV = cross(v0 + v1, e1);
129 const float3 crossW = cross(v1 + v2, e2);
130
131 ssef crossX(crossU);
132 ssef crossY(crossV);
133 ssef crossZ(crossW);
134 ssef zero = _mm_setzero_ps();
135 _MM_TRANSPOSE4_PS(crossX, crossY, crossZ, zero);
136
137 const ssef dirX(ray_dir.x);
138 const ssef dirY(ray_dir.y);
139 const ssef dirZ(ray_dir.z);
140
141 ssef UVWW = madd(crossX, dirX, madd(crossY, dirY, crossZ * dirZ));
142 #else /* __KERNEL_SSE2__ */
143 const float U = dot(cross(v2 + v0, e0), ray_dir);
144 const float V = dot(cross(v0 + v1, e1), ray_dir);
145 const float W = dot(cross(v1 + v2, e2), ray_dir);
146 #endif /* __KERNEL_SSE2__ */
147
148 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
149 int uvw_sign = movemask(UVWW) & 0x7;
150 if (uvw_sign != 0) {
151 if (uvw_sign != 0x7) {
152 return false;
153 }
154 }
155 #else
156 const float minUVW = min(U, min(V, W));
157 const float maxUVW = max(U, max(V, W));
158
159 if (minUVW < 0.0f && maxUVW > 0.0f) {
160 return false;
161 }
162 #endif
163
164 /* Calculate geometry normal and denominator. */
165 const float3 Ng1 = cross(e1, e0);
166 // const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0);
167 const float3 Ng = Ng1 + Ng1;
168 const float den = dot3(Ng, dir);
169 /* Avoid division by 0. */
170 if (UNLIKELY(den == 0.0f)) {
171 return false;
172 }
173
174 /* Perform depth test. */
175 const float T = dot3(v0, Ng);
176 const int sign_den = (__float_as_int(den) & 0x80000000);
177 const float sign_T = xor_signmask(T, sign_den);
178 if ((sign_T < 0.0f) || (sign_T > ray_t * xor_signmask(den, sign_den))) {
179 return false;
180 }
181
182 const float inv_den = 1.0f / den;
183 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
184 UVWW *= inv_den;
185 _mm_store_ss(isect_u, UVWW);
186 _mm_store_ss(isect_v, shuffle<1, 1, 3, 3>(UVWW));
187 #else
188 *isect_u = U * inv_den;
189 *isect_v = V * inv_den;
190 #endif
191 *isect_t = T * inv_den;
192 return true;
193
194 #undef dot3
195 }
196
197 /* Tests for an intersection between a ray and a quad defined by
198 * its midpoint, normal and sides.
199 * If ellipse is true, hits outside the ellipse that's enclosed by the
200 * quad are rejected.
201 */
ray_quad_intersect(float3 ray_P,float3 ray_D,float ray_mint,float ray_maxt,float3 quad_P,float3 quad_u,float3 quad_v,float3 quad_n,float3 * isect_P,float * isect_t,float * isect_u,float * isect_v,bool ellipse)202 ccl_device bool ray_quad_intersect(float3 ray_P,
203 float3 ray_D,
204 float ray_mint,
205 float ray_maxt,
206 float3 quad_P,
207 float3 quad_u,
208 float3 quad_v,
209 float3 quad_n,
210 float3 *isect_P,
211 float *isect_t,
212 float *isect_u,
213 float *isect_v,
214 bool ellipse)
215 {
216 /* Perform intersection test. */
217 float t = -(dot(ray_P, quad_n) - dot(quad_P, quad_n)) / dot(ray_D, quad_n);
218 if (t < ray_mint || t > ray_maxt) {
219 return false;
220 }
221 const float3 hit = ray_P + t * ray_D;
222 const float3 inplane = hit - quad_P;
223 const float u = dot(inplane, quad_u) / dot(quad_u, quad_u);
224 if (u < -0.5f || u > 0.5f) {
225 return false;
226 }
227 const float v = dot(inplane, quad_v) / dot(quad_v, quad_v);
228 if (v < -0.5f || v > 0.5f) {
229 return false;
230 }
231 if (ellipse && (u * u + v * v > 0.25f)) {
232 return false;
233 }
234 /* Store the result. */
235 /* TODO(sergey): Check whether we can avoid some checks here. */
236 if (isect_P != NULL)
237 *isect_P = hit;
238 if (isect_t != NULL)
239 *isect_t = t;
240 if (isect_u != NULL)
241 *isect_u = u + 0.5f;
242 if (isect_v != NULL)
243 *isect_v = v + 0.5f;
244 return true;
245 }
246
247 CCL_NAMESPACE_END
248
249 #endif /* __UTIL_MATH_INTERSECT_H__ */
250