498c21b9706cc6da04b919e4bbd3a2a0e9dfb0d8
[blender.git] / intern / cycles / util / util_math_intersect.h
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
24 ccl_device bool ray_sphere_intersect(
25         float3 ray_P, float3 ray_D, float ray_t,
26         float3 sphere_P, float sphere_radius,
27         float3 *isect_P, float *isect_t)
28 {
29         const float3 d = sphere_P - ray_P;
30         const float radiussq = sphere_radius*sphere_radius;
31         const float tsq = dot(d, d);
32
33         if(tsq > radiussq) {
34                 /* Ray origin outside sphere. */
35                 const float tp = dot(d, ray_D);
36                 if(tp < 0.0f) {
37                         /* Ray  points away from sphere. */
38                         return false;
39                 }
40                 const float dsq = tsq - tp*tp;  /* pythagoras */
41                 if(dsq > radiussq)  {
42                         /* Closest point on ray outside sphere. */
43                         return false;
44                 }
45                 const float t = tp - sqrtf(radiussq - dsq);  /* pythagoras */
46                 if(t < ray_t) {
47                         *isect_t = t;
48                         *isect_P = ray_P + ray_D*t;
49                         return true;
50                 }
51         }
52         return false;
53 }
54
55 ccl_device bool ray_aligned_disk_intersect(
56         float3 ray_P, float3 ray_D, float ray_t,
57         float3 disk_P, float disk_radius,
58         float3 *isect_P, float *isect_t)
59 {
60         /* Aligned disk normal. */
61         float disk_t;
62         const float3 disk_N = normalize_len(ray_P - disk_P, &disk_t);
63         const float div = dot(ray_D, disk_N);
64         if(UNLIKELY(div == 0.0f)) {
65                 return false;
66         }
67         /* Compute t to intersection point. */
68         const float t = -disk_t/div;
69         if(t < 0.0f || t > ray_t) {
70                 return false;
71         }
72         /* Test if within radius. */
73         float3 P = ray_P + ray_D*t;
74         if(len_squared(P - disk_P) > disk_radius*disk_radius) {
75                 return false;
76         }
77         *isect_P = P;
78         *isect_t = t;
79         return true;
80 }
81
82 #if defined(__KERNEL_CUDA__) && __CUDA_ARCH__ < 300
83 ccl_device_inline
84 #else
85 ccl_device_forceinline
86 #endif
87 bool ray_triangle_intersect(
88         float3 ray_P, float3 ray_dir, float ray_t,
89 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
90         const ssef *ssef_verts,
91 #else
92         const float3 tri_a, const float3 tri_b, const float3 tri_c,
93 #endif
94         float *isect_u, float *isect_v, float *isect_t)
95 {
96 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
97         typedef ssef float3;
98         const float3 tri_a(ssef_verts[0]);
99         const float3 tri_b(ssef_verts[1]);
100         const float3 tri_c(ssef_verts[2]);
101         const float3 P(ray_P);
102         const float3 dir(ray_dir);
103 #else
104 #  define dot3(a, b) dot(a, b)
105         const float3 P = ray_P;
106         const float3 dir = ray_dir;
107 #endif
108
109         /* Calculate vertices relative to ray origin. */
110         const float3 v0 = tri_c - P;
111         const float3 v1 = tri_a - P;
112         const float3 v2 = tri_b - P;
113
114         /* Calculate triangle edges. */
115         const float3 e0 = v2 - v0;
116         const float3 e1 = v0 - v1;
117         const float3 e2 = v1 - v2;
118
119         /* Perform edge tests. */
120 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
121         const float3 crossU = cross(v2 + v0, e0);
122         const float3 crossV = cross(v0 + v1, e1);
123         const float3 crossW = cross(v1 + v2, e2);
124
125         ssef crossX(crossU);
126         ssef crossY(crossV);
127         ssef crossZ(crossW);
128         ssef zero = _mm_setzero_ps();
129         _MM_TRANSPOSE4_PS(crossX, crossY, crossZ, zero);
130
131         const ssef dirX(ray_dir.x);
132         const ssef dirY(ray_dir.y);
133         const ssef dirZ(ray_dir.z);
134
135         ssef UVWW = madd(crossX, dirX, madd(crossY, dirY, crossZ * dirZ));
136 #else  /* __KERNEL_SSE2__ */
137         const float U = dot(cross(v2 + v0, e0), ray_dir);
138         const float V = dot(cross(v0 + v1, e1), ray_dir);
139         const float W = dot(cross(v1 + v2, e2), ray_dir);
140 #endif  /* __KERNEL_SSE2__ */
141
142 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
143         int uvw_sign = movemask(UVWW) & 0x7;
144         if (uvw_sign != 0)
145         {
146                 if (uvw_sign != 0x7)
147                         return false;
148         }
149 #else
150         const float minUVW = min(U, min(V, W));
151         const float maxUVW = max(U, max(V, W));
152
153         if(minUVW < 0.0f && maxUVW > 0.0f) {
154                 return false;
155         }
156 #endif
157
158
159         /* Calculate geometry normal and denominator. */
160         const float3 Ng1 = cross(e1, e0);
161         //const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0);
162         const float3 Ng = Ng1 + Ng1;
163         const float den = dot3(Ng, dir);
164         /* Avoid division by 0. */
165         if(UNLIKELY(den == 0.0f)) {
166                 return false;
167         }
168
169         /* Perform depth test. */
170         const float T = dot3(v0, Ng);
171         const int sign_den = (__float_as_int(den) & 0x80000000);
172         const float sign_T = xor_signmask(T, sign_den);
173         if((sign_T < 0.0f) ||
174            (sign_T > ray_t * xor_signmask(den, sign_den)))
175         {
176                 return false;
177         }
178
179         const float inv_den = 1.0f / den;
180 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
181         UVWW *= inv_den;
182         _mm_store_ss(isect_u, UVWW);
183         _mm_store_ss(isect_v, shuffle<1,1,3,3>(UVWW));
184 #else
185         *isect_u = U * inv_den;
186         *isect_v = V * inv_den;
187 #endif
188         *isect_t = T * inv_den;
189         return true;
190
191 #undef dot3
192 }
193
194 ccl_device bool ray_quad_intersect(float3 ray_P, float3 ray_D,
195                                    float ray_mint, float ray_maxt,
196                                    float3 quad_P,
197                                    float3 quad_u, float3 quad_v, float3 quad_n,
198                                    float3 *isect_P, float *isect_t,
199                                    float *isect_u, float *isect_v)
200 {
201         /* Perform intersection test. */
202         float t = -(dot(ray_P, quad_n) - dot(quad_P, quad_n)) / dot(ray_D, quad_n);
203         if(t < ray_mint || t > ray_maxt) {
204                 return false;
205         }
206         const float3 hit = ray_P + t*ray_D;
207         const float3 inplane = hit - quad_P;
208         const float u = dot(inplane, quad_u) / dot(quad_u, quad_u) + 0.5f;
209         if(u < 0.0f || u > 1.0f) {
210                 return false;
211         }
212         const float v = dot(inplane, quad_v) / dot(quad_v, quad_v) + 0.5f;
213         if(v < 0.0f || v > 1.0f) {
214                 return false;
215         }
216         /* Store the result. */
217         /* TODO(sergey): Check whether we can avoid some checks here. */
218         if(isect_P != NULL) *isect_P = hit;
219         if(isect_t != NULL) *isect_t = t;
220         if(isect_u != NULL) *isect_u = u;
221         if(isect_v != NULL) *isect_v = v;
222         return true;
223 }
224
225 CCL_NAMESPACE_END
226
227 #endif /* __UTIL_MATH_INTERSECT_H__ */