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 &center,
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