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 ccl_device_forceinline bool ray_triangle_intersect(
  83         float3 ray_P, float3 ray_dir, float ray_t,
  84 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
  85         const ssef *ssef_verts,
  86 #else
  87         const float3 tri_a, const float3 tri_b, const float3 tri_c,
  88 #endif
  89         float *isect_u, float *isect_v, float *isect_t)
  90 {
  91 #if defined(__KERNEL_SSE2__) && defined(__KERNEL_SSE__)
  92         typedef ssef float3;
  93         const float3 tri_a(ssef_verts[0]);
  94         const float3 tri_b(ssef_verts[1]);
  95         const float3 tri_c(ssef_verts[2]);
  96         const float3 P(ray_P);
  97         const float3 dir(ray_dir);
  98 #else
  99 #  define dot3(a, b) dot(a, b)
 100         const float3 P = ray_P;
 101         const float3 dir = ray_dir;
 102 #endif
 103
 104         /* Calculate vertices relative to ray origin. */
 105         const float3 v0 = tri_c - P;
 106         const float3 v1 = tri_a - P;
 107         const float3 v2 = tri_b - P;
 108
 109         /* Calculate triangle edges. */
 110         const float3 e0 = v2 - v0;
 111         const float3 e1 = v0 - v1;
 112         const float3 e2 = v1 - v2;
 113
 114         /* Perform edge tests. */
 115 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
 116         const float3 crossU = cross(v2 + v0, e0);
 117         const float3 crossV = cross(v0 + v1, e1);
 118         const float3 crossW = cross(v1 + v2, e2);
 119
 120         ssef crossX(crossU);
 121         ssef crossY(crossV);
 122         ssef crossZ(crossW);
 123         ssef zero = _mm_setzero_ps();
 124         _MM_TRANSPOSE4_PS(crossX, crossY, crossZ, zero);
 125
 126         const ssef dirX(ray_dir.x);
 127         const ssef dirY(ray_dir.y);
 128         const ssef dirZ(ray_dir.z);
 129
 130         ssef UVWW = madd(crossX, dirX, madd(crossY, dirY, crossZ * dirZ));
 131 #else  /* __KERNEL_SSE2__ */
 132         const float U = dot(cross(v2 + v0, e0), ray_dir);
 133         const float V = dot(cross(v0 + v1, e1), ray_dir);
 134         const float W = dot(cross(v1 + v2, e2), ray_dir);
 135 #endif  /* __KERNEL_SSE2__ */
 136
 137 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
 138         int uvw_sign = movemask(UVWW) & 0x7;
 139         if (uvw_sign != 0)
 140         {
 141                 if (uvw_sign != 0x7)
 142                         return false;
 143         }
 144 #else
 145         const float minUVW = min(U, min(V, W));
 146         const float maxUVW = max(U, max(V, W));
 147
 148         if(minUVW < 0.0f && maxUVW > 0.0f) {
 149                 return false;
 150         }
 151 #endif
 152
 153
 154         /* Calculate geometry normal and denominator. */
 155         const float3 Ng1 = cross(e1, e0);
 156         //const Vec3vfM Ng1 = stable_triangle_normal(e2,e1,e0);
 157         const float3 Ng = Ng1 + Ng1;
 158         const float den = dot3(Ng, dir);
 159         /* Avoid division by 0. */
 160         if(UNLIKELY(den == 0.0f)) {
 161                 return false;
 162         }
 163
 164         /* Perform depth test. */
 165         const float T = dot3(v0, Ng);
 166         const int sign_den = (__float_as_int(den) & 0x80000000);
 167         const float sign_T = xor_signmask(T, sign_den);
 168         if((sign_T < 0.0f) ||
 169            (sign_T > ray_t * xor_signmask(den, sign_den)))
 170         {
 171                 return false;
 172         }
 173
 174         const float inv_den = 1.0f / den;
 175 #if defined(__KERNEL_SSE2__)  && defined (__KERNEL_SSE__)
 176         UVWW *= inv_den;
 177         _mm_store_ss(isect_u, UVWW);
 178         _mm_store_ss(isect_v, shuffle<1,1,3,3>(UVWW));
 179 #else
 180         *isect_u = U * inv_den;
 181         *isect_v = V * inv_den;
 182 #endif
 183         *isect_t = T * inv_den;
 184         return true;
 185
 186 #undef dot3
 187 }
 188
 189 ccl_device bool ray_quad_intersect(float3 ray_P, float3 ray_D,
 190                                    float ray_mint, float ray_maxt,
 191                                    float3 quad_P,
 192                                    float3 quad_u, float3 quad_v, float3 quad_n,
 193                                    float3 *isect_P, float *isect_t,
 194                                    float *isect_u, float *isect_v)
 195 {
 196         /* Perform intersection test. */
 197         float t = -(dot(ray_P, quad_n) - dot(quad_P, quad_n)) / dot(ray_D, quad_n);
 198         if(t < ray_mint || t > ray_maxt) {
 199                 return false;
 200         }
 201         const float3 hit = ray_P + t*ray_D;
 202         const float3 inplane = hit - quad_P;
 203         const float u = dot(inplane, quad_u) / dot(quad_u, quad_u) + 0.5f;
 204         if(u < 0.0f || u > 1.0f) {
 205                 return false;
 206         }
 207         const float v = dot(inplane, quad_v) / dot(quad_v, quad_v) + 0.5f;
 208         if(v < 0.0f || v > 1.0f) {
 209                 return false;
 210         }
 211         /* Store the result. */
 212         /* TODO(sergey): Check whether we can avoid some checks here. */
 213         if(isect_P != NULL) *isect_P = hit;
 214         if(isect_t != NULL) *isect_t = t;
 215         if(isect_u != NULL) *isect_u = u;
 216         if(isect_v != NULL) *isect_v = v;
 217         return true;
 218 }
 219
 220 CCL_NAMESPACE_END
 221
 222 #endif /* __UTIL_MATH_INTERSECT_H__ */