1// Copyright 2009-2021 Intel Corporation 2// SPDX-License-Identifier: Apache-2.0 3 4#pragma once 5 6#include "rkcommon/math/box.ih" 7 8// Ray intersection structures ////////////////////////////////////////////// 9 10struct Hit 11{ 12 bool hit; 13 float t; 14 vec3f N; 15 float u; 16}; 17 18struct Intersections 19{ 20 Hit entry; 21 Hit exit; 22}; 23 24// Ray intersection helpers /////////////////////////////////////////////////// 25 26// robust ray-sphere intersection 27inline Intersections intersectSphere(const vec3f &rayOrg, 28 const vec3f &rayDir, 29 const uniform vec3f ¢er, 30 const uniform float radius) 31{ 32 Intersections isect; 33 isect.entry.hit = false; 34 isect.exit.hit = false; 35 isect.entry.t = inf; 36 isect.exit.t = -inf; 37 38 const vec3f d = rayDir; 39 const float rd2 = 1.0f / dot(d, d); // 1/a 40 const vec3f CO = center - rayOrg; 41 // transformation to avoid missing a small sphere which is far away: 42 // the standard c=CO^2-r^2 would quickly loose term r due to float arithmetic 43 const float projCO = dot(CO, d) * rd2; // in ray-space 44 const vec3f perp = CO - projCO * d; 45 const float l2 = dot(perp, perp); 46 const uniform float r2 = sqr(radius); 47 if (l2 > r2) 48 return isect; 49 float td = sqrt((r2 - l2) * rd2); 50 isect.entry.hit = true; 51 isect.exit.hit = true; 52 isect.entry.t = projCO - td; 53 isect.exit.t = projCO + td; 54 55 // above solutions are problematic if rays starts close to the sphere 56 // (due to catastrophic cancellation, because then |projCO| ~ td) 57 // the usual recommendation is to choose the one solution with same sign: 58 // const float t1 = projCO + floatbits(signbits(projCO)|intbits(td)); 59 // and compute the other solution via t1*t2=c/a: 60 // const float t2 = (dot(CO, CO) - r2) / t1 * rd2; 61 // this is more precise, but still problematic in particular for large 62 // spheres, because |CO| ~ r; slightly better alternative, but costly sqrt: 63 // const float f = sqrt(dot(CO, CO)); 64 // const float t2 = (f - radius) * (f + radius) / t1 * rd2; 65 // the only variant I found that has high enough precision to avoid 66 // self-intersections of 2ndary rays is to (re-)compute most terms (CO, dot, 67 // r2, t2) with doubles; large spheres are a rare usecase for OSPRay, thus we 68 // use instead as a workaround an additional, radius-dependent epsilon 69 70 // cannot easily be moved to postIntersect 71 // we need hit in object space, in postIntersect it is in world-space 72 isect.entry.N = -td * d - perp; 73 isect.exit.N = td * d - perp; 74 75 return isect; 76} 77 78inline Intersections intersectCylinder(const vec3f &rayOrg, 79 const vec3f &rayDir, 80 const uniform vec3f &v0, 81 const uniform vec3f &v1, 82 const uniform float radius) 83{ 84 Intersections isect; 85 isect.entry.hit = false; 86 isect.exit.hit = false; 87 isect.entry.t = inf; 88 isect.exit.t = -inf; 89 90 const uniform vec3f cZ = v1 - v0; 91 const vec3f q = rayOrg - v0; 92 93 const uniform float z2 = dot(cZ, cZ); 94 const float d = dot(cZ, rayDir); 95 const float c = dot(cZ, q); 96 97 const float A = z2 - sqr(d); 98 const float B = z2 * dot(q, rayDir) - c * d; 99 const float C = z2 * dot(q, q) - sqr(c) - sqr(radius) * z2; 100 101 float radical = B * B - A * C; 102 if (radical < 0.f) { 103 return isect; 104 } 105 106 radical = sqrt(radical); 107 108 const float tin = (-B - radical) / A; 109 const float tout = (-B + radical) / A; 110 111 // first hit 112 const float yin = c + tin * d; 113 if (yin > 0.f && yin < z2) { 114 // body hit 115 isect.entry.hit = true; 116 isect.entry.t = tin; 117 isect.entry.u = yin * rcp(z2); 118 isect.entry.N = (q + tin * rayDir - cZ * yin * rcp(z2)) * rcp(radius); 119 } else { 120 const float tcapin = (((yin < 0.f) ? 0.f : z2) - c) / d; 121 if (abs(B + A * tcapin) < radical) { 122 // cap hit 123 isect.entry.hit = false; 124 isect.entry.t = tin; 125 isect.entry.u = (yin < 0.f) ? 0.f : 1.f; 126 const float us = signbits(yin) ? -1. : 1.; 127 isect.entry.N = cZ * us / z2; 128 } 129 } 130 131 // second hit 132 const float yout = c + tout * d; 133 if (yout > 0.f && yout < z2) { 134 // body hit 135 isect.exit.hit = true; 136 isect.exit.t = tout; 137 isect.exit.u = yout * rcp(z2); 138 isect.exit.N = (q + tout * rayDir - cZ * yout * rcp(z2)) * rcp(radius); 139 } else { 140 const float tcapout = (((yout < 0.f) ? 0.f : z2) - c) / d; 141 if (abs(B + A * tcapout) < radical) { 142 // cap hit 143 isect.exit.hit = false; 144 isect.exit.t = tout; 145 isect.exit.u = (yout < 0.f) ? 0.f : 1.f; 146 const float us = signbits(yout) ? -1. : 1.; 147 isect.exit.N = cZ * us / z2; 148 } 149 } 150 151 return isect; 152} 153 154inline Intersections intersectCapsule(const vec3f &rayOrg, 155 const vec3f &rayDir, 156 const uniform vec3f &v0, 157 const uniform vec3f &v1, 158 const uniform float radius) 159{ 160 Intersections isect_pipe = intersectCylinder(rayOrg, rayDir, v0, v1, radius); 161 const Intersections isect_sph1 = intersectSphere(rayOrg, rayDir, v0, radius); 162 const Intersections isect_sph2 = intersectSphere(rayOrg, rayDir, v1, radius); 163 164 const float t_in = 165 min(min(isect_sph1.entry.t, isect_sph2.entry.t), isect_pipe.entry.t); 166 const float t_out = 167 max(max(isect_sph1.exit.t, isect_sph2.exit.t), isect_pipe.exit.t); 168 169 isect_pipe.entry.hit |= isect_sph1.entry.hit | isect_sph2.entry.hit; 170 isect_pipe.entry.t = t_in; 171 isect_pipe.exit.hit |= isect_sph1.exit.hit | isect_sph2.exit.hit; 172 isect_pipe.exit.t = t_out; 173 174 if (isect_sph1.entry.t == t_in) { 175 isect_pipe.entry.u = 0.f; 176 isect_pipe.entry.N = isect_sph1.entry.N; 177 } else if (isect_sph2.entry.t == t_in) { 178 isect_pipe.entry.u = 1.f; 179 isect_pipe.entry.N = isect_sph2.entry.N; 180 } 181 182 if (isect_sph1.exit.t == t_out) { 183 isect_pipe.exit.u = 0.f; 184 isect_pipe.exit.N = isect_sph1.exit.N; 185 } else if (isect_sph2.exit.t == t_out) { 186 isect_pipe.exit.u = 1.f; 187 isect_pipe.exit.N = isect_sph2.exit.N; 188 } 189 190 return isect_pipe; 191} 192 193inline Intersections intersectBox( 194 const vec3f &rayOrg, const vec3f &rayDir, const uniform box3f &box) 195{ 196 Intersections isect; 197 198 const vec3f mins = (box.lower - rayOrg) * rcp_safe(rayDir); 199 const vec3f maxs = (box.upper - rayOrg) * rcp_safe(rayDir); 200 const vec3f nears = min(mins, maxs); 201 const vec3f fars = max(mins, maxs); 202 203 isect.entry.t = reduce_max(nears); 204 if (isect.entry.t == nears.x) 205 isect.entry.N = make_vec3f(rayDir.x > 0.0f ? -1.0f : 1.0f, 0.0f, 0.0f); 206 else if (isect.entry.t == nears.y) 207 isect.entry.N = make_vec3f(0.0f, rayDir.y > 0.0f ? -1.0f : 1.0f, 0.0f); 208 else 209 isect.entry.N = make_vec3f(0.0f, 0.0f, rayDir.z > 0.0f ? -1.0f : 1.0f); 210 isect.exit.t = reduce_min(fars); 211 if (isect.exit.t == fars.x) 212 isect.exit.N = make_vec3f(rayDir.x > 0.0f ? 1.0f : -1.0f, 0.0f, 0.0f); 213 else if (isect.exit.t == fars.y) 214 isect.exit.N = make_vec3f(0.0f, rayDir.y > 0.0f ? 1.0f : -1.0f, 0.0f); 215 else 216 isect.exit.N = make_vec3f(0.0f, 0.0f, rayDir.z > 0.0f ? 1.0f : -1.0f); 217 isect.entry.hit = isect.entry.t < isect.exit.t; 218 isect.exit.hit = isect.entry.hit; 219 220 return isect; 221} 222 223inline Hit intersectPlane( 224 const vec3f &rayOrg, const vec3f &rayDir, const uniform vec4f &plane) 225{ 226 Hit hit; 227 hit.hit = false; 228 const uniform vec3f normal = make_vec3f(plane); 229 const float DdN = dot(rayDir, normal); 230 hit.hit = DdN != 0.0f; 231 hit.N = normal; 232 hit.t = (plane.w - dot(rayOrg, normal)) * rcpf(DdN); 233 234 return hit; 235} 236