A game about forced loneliness, made by TACStudios
tacstudios.tngl.sh
1#ifndef UNITY_GEOMETRICTOOLS_INCLUDED
2#define UNITY_GEOMETRICTOOLS_INCLUDED
3
4//-----------------------------------------------------------------------------
5// Transform functions
6//-----------------------------------------------------------------------------
7
8// Rotate around a pivot point and an axis
9float3 Rotate(float3 pivot, float3 position, float3 rotationAxis, float angle)
10{
11 rotationAxis = normalize(rotationAxis);
12 float3 cpa = pivot + rotationAxis * dot(rotationAxis, position - pivot);
13 return cpa + ((position - cpa) * cos(angle) + cross(rotationAxis, (position - cpa)) * sin(angle));
14}
15
16float3x3 RotationFromAxisAngle(float3 A, float sinAngle, float cosAngle)
17{
18 float c = cosAngle;
19 float s = sinAngle;
20
21 return float3x3(A.x * A.x * (1 - c) + c, A.x * A.y * (1 - c) - A.z * s, A.x * A.z * (1 - c) + A.y * s,
22 A.x * A.y * (1 - c) + A.z * s, A.y * A.y * (1 - c) + c, A.y * A.z * (1 - c) - A.x * s,
23 A.x * A.z * (1 - c) - A.y * s, A.y * A.z * (1 - c) + A.x * s, A.z * A.z * (1 - c) + c);
24}
25
26//-----------------------------------------------------------------------------
27// Solver
28//-----------------------------------------------------------------------------
29
30// Solves the quadratic equation of the form: a*t^2 + b*t + c = 0.
31// Returns 'false' if there are no real roots, 'true' otherwise.
32// Ensures that roots.x <= roots.y.
33bool SolveQuadraticEquation(float a, float b, float c, out float2 roots)
34{
35 float det = Sq(b) - 4.0 * a * c;
36
37 float sqrtDet = sqrt(det);
38 roots.x = (-b - sign(a) * sqrtDet) / (2.0 * a);
39 roots.y = (-b + sign(a) * sqrtDet) / (2.0 * a);
40
41 return (det >= 0.0);
42}
43
44//-----------------------------------------------------------------------------
45// Intersection functions
46//-----------------------------------------------------------------------------
47
48bool IntersectRayAABB(float3 rayOrigin, float3 rayDirection,
49 float3 boxMin, float3 boxMax,
50 float tMin, float tMax,
51 out float tEntr, out float tExit)
52{
53 // Could be precomputed. Clamp to avoid INF. clamp() is a single ALU on GCN.
54 // rcp(FLT_EPS) = 16,777,216, which is large enough for our purposes,
55 // yet doesn't cause a lot of numerical issues associated with FLT_MAX.
56 float3 rayDirInv = clamp(rcp(rayDirection), -rcp(FLT_EPS), rcp(FLT_EPS));
57
58 // Perform ray-slab intersection (component-wise).
59 float3 t0 = boxMin * rayDirInv - (rayOrigin * rayDirInv);
60 float3 t1 = boxMax * rayDirInv - (rayOrigin * rayDirInv);
61
62 // Find the closest/farthest distance (component-wise).
63 float3 tSlabEntr = min(t0, t1);
64 float3 tSlabExit = max(t0, t1);
65
66 // Find the farthest entry and the nearest exit.
67 tEntr = Max3(tSlabEntr.x, tSlabEntr.y, tSlabEntr.z);
68 tExit = Min3(tSlabExit.x, tSlabExit.y, tSlabExit.z);
69
70 // Clamp to the range.
71 tEntr = max(tEntr, tMin);
72 tExit = min(tExit, tMax);
73
74 return tEntr < tExit;
75}
76
77// This simplified version assume that we care about the result only when we are inside the box
78float IntersectRayAABBSimple(float3 start, float3 dir, float3 boxMin, float3 boxMax)
79{
80 float3 invDir = rcp(dir);
81
82 // Find the ray intersection with box plane
83 float3 rbmin = (boxMin - start) * invDir;
84 float3 rbmax = (boxMax - start) * invDir;
85
86 float3 rbminmax = float3((dir.x > 0.0) ? rbmax.x : rbmin.x, (dir.y > 0.0) ? rbmax.y : rbmin.y, (dir.z > 0.0) ? rbmax.z : rbmin.z);
87
88 return min(min(rbminmax.x, rbminmax.y), rbminmax.z);
89}
90
91// Assume Sphere is at the origin (i.e start = position - spherePosition)
92bool IntersectRaySphere(float3 start, float3 dir, float radius, out float2 intersections)
93{
94 float a = dot(dir, dir);
95 float b = dot(dir, start) * 2.0;
96 float c = dot(start, start) - radius * radius;
97
98 return SolveQuadraticEquation(a, b, c, intersections);
99}
100
101// This simplified version assume that we care about the result only when we are inside the sphere
102// Assume Sphere is at the origin (i.e start = position - spherePosition) and dir is normalized
103// Ref: http://http.developer.nvidia.com/GPUGems/gpugems_ch19.html
104float IntersectRaySphereSimple(float3 start, float3 dir, float radius)
105{
106 float b = dot(dir, start) * 2.0;
107 float c = dot(start, start) - radius * radius;
108 float discriminant = b * b - 4.0 * c;
109
110 return abs(sqrt(discriminant) - b) * 0.5;
111}
112
113float3 IntersectRayPlane(float3 rayOrigin, float3 rayDirection, float3 planeOrigin, float3 planeNormal)
114{
115 float dist = dot(planeNormal, planeOrigin - rayOrigin) / dot(planeNormal, rayDirection);
116 return rayOrigin + rayDirection * dist;
117}
118
119// Same as above but return intersection distance and true / false if the ray hit/miss
120bool IntersectRayPlane(float3 rayOrigin, float3 rayDirection, float3 planePosition, float3 planeNormal, out float t)
121{
122 bool res = false;
123 t = -1.0;
124
125 float denom = dot(planeNormal, rayDirection);
126 if (abs(denom) > 1e-5)
127 {
128 float3 d = planePosition - rayOrigin;
129 t = dot(d, planeNormal) / denom;
130 res = (t >= 0);
131 }
132
133 return res;
134}
135
136bool RayPlaneSegmentIntersect(in float3 rayOrigin, in float3 rayDirection, float3 planeNormal, float planeDistance, inout float3 intersectionPoint)
137{
138 float denom = dot(rayDirection, planeNormal);
139 float lambda = (denom != 0.0) ? (planeDistance - dot(rayOrigin, planeNormal)) / denom : -1.0;
140 if ((lambda >= 0.0) && (lambda <= 1.0))
141 {
142 intersectionPoint = rayOrigin + lambda * rayDirection;
143 return true;
144 }
145 else
146 return false;
147}
148
149// Can support cones with an elliptic base: pre-scale 'coneAxisX' and 'coneAxisY' by (h/r_x) and (h/r_y).
150// Returns parametric distances 'tEntr' and 'tExit' along the ray,
151// subject to constraints 'tMin' and 'tMax'.
152bool IntersectRayCone(float3 rayOrigin, float3 rayDirection,
153 float3 coneOrigin, float3 coneDirection,
154 float3 coneAxisX, float3 coneAxisY,
155 float tMin, float tMax,
156 out float tEntr, out float tExit)
157{
158 // Inverse transform the ray into a coordinate system with the cone at the origin facing along the Z axis.
159 float3x3 rotMat = float3x3(coneAxisX, coneAxisY, coneDirection);
160
161 float3 o = mul(rotMat, rayOrigin - coneOrigin);
162 float3 d = mul(rotMat, rayDirection);
163
164 // Cone equation (facing along Z): (h/r*x)^2 + (h/r*y)^2 - z^2 = 0.
165 // Cone axes are premultiplied with (h/r).
166 // Set up the quadratic equation: a*t^2 + b*t + c = 0.
167 float a = d.x * d.x + d.y * d.y - d.z * d.z;
168 float b = o.x * d.x + o.y * d.y - o.z * d.z;
169 float c = o.x * o.x + o.y * o.y - o.z * o.z;
170
171 float2 roots;
172
173 // Check whether we have at least 1 root.
174 bool hit = SolveQuadraticEquation(a, 2 * b, c, roots);
175
176 tEntr = roots.x;
177 tExit = roots.y;
178 float3 pEntr = o + tEntr * d;
179 float3 pExit = o + tExit * d;
180
181 // Clip the negative cone.
182 bool pEntrNeg = pEntr.z < 0;
183 bool pExitNeg = pExit.z < 0;
184 if (pEntrNeg && pExitNeg) { hit = false; }
185 if (pEntrNeg) { tEntr = tExit; tExit = tMax; }
186 if (pExitNeg) { tExit = tEntr; tEntr = tMin; }
187
188 // Clamp using the values passed into the function.
189 tEntr = clamp(tEntr, tMin, tMax);
190 tExit = clamp(tExit, tMin, tMax);
191
192 // Check for grazing intersections.
193 if (tEntr == tExit) { hit = false; }
194
195 return hit;
196}
197
198bool IntersectSphereAABB(float3 position, float radius, float3 aabbMin, float3 aabbMax)
199{
200 float x = max(aabbMin.x, min(position.x, aabbMax.x));
201 float y = max(aabbMin.y, min(position.y, aabbMax.y));
202 float z = max(aabbMin.z, min(position.z, aabbMax.z));
203 float distance2 = ((x - position.x) * (x - position.x) + (y - position.y) * (y - position.y) + (z - position.z) * (z - position.z));
204 return distance2 < radius * radius;
205}
206
207//-----------------------------------------------------------------------------
208// Miscellaneous functions
209//-----------------------------------------------------------------------------
210
211// Box is AABB
212float DistancePointBox(float3 position, float3 boxMin, float3 boxMax)
213{
214 return length(max(max(position - boxMax, boxMin - position), float3(0.0, 0.0, 0.0)));
215}
216
217float3 ProjectPointOnPlane(float3 position, float3 planePosition, float3 planeNormal)
218{
219 return position - (dot(position - planePosition, planeNormal) * planeNormal);
220}
221
222// Plane equation: {(a, b, c) = N, d = -dot(N, P)}.
223// Returns the distance from the plane to the point 'p' along the normal.
224// Positive -> in front (above), negative -> behind (below).
225float DistanceFromPlane(float3 p, float4 plane)
226{
227 return dot(float4(p, 1.0), plane);
228}
229
230// Returns 'true' if the triangle is outside of the frustum.
231// 'epsilon' is the (negative) distance to (outside of) the frustum below which we cull the triangle.
232bool CullTriangleFrustum(float3 p0, float3 p1, float3 p2, float epsilon, float4 frustumPlanes[6], int numPlanes)
233{
234 bool outside = false;
235
236 for (int i = 0; i < numPlanes; i++)
237 {
238 // If all 3 points are behind any of the planes, we cull.
239 outside = outside || Max3(DistanceFromPlane(p0, frustumPlanes[i]),
240 DistanceFromPlane(p1, frustumPlanes[i]),
241 DistanceFromPlane(p2, frustumPlanes[i])) < epsilon;
242 }
243
244 return outside;
245}
246
247// Returns 'true' if the edge of the triangle is outside of the frustum.
248// The edges are defined s.t. they are on the opposite side of the point with the given index.
249// 'epsilon' is the (negative) distance to (outside of) the frustum below which we cull the triangle.
250//output packing:
251// x,y,z - one component per triangle edge, true if outside, false otherwise
252// w - true if entire triangle is outside of at least 1 plane of the frustum, false otherwise
253bool4 CullFullTriangleAndEdgesFrustum(float3 p0, float3 p1, float3 p2, float epsilon, float4 frustumPlanes[6], int numPlanes)
254{
255 bool4 edgesOutsideXYZ_triangleOutsideW = false;
256
257 for (int i = 0; i < numPlanes; i++)
258 {
259 bool3 pointsOutside = bool3(DistanceFromPlane(p0, frustumPlanes[i]) < epsilon,
260 DistanceFromPlane(p1, frustumPlanes[i]) < epsilon,
261 DistanceFromPlane(p2, frustumPlanes[i]) < epsilon);
262
263 bool3 edgesOutside;
264 // If both points of the edge are behind any of the planes, we cull.
265 edgesOutside.x = pointsOutside.y && pointsOutside.z;
266 edgesOutside.y = pointsOutside.x && pointsOutside.z;
267 edgesOutside.z = pointsOutside.x && pointsOutside.y;
268
269 edgesOutsideXYZ_triangleOutsideW = bool4(edgesOutsideXYZ_triangleOutsideW.x || edgesOutside.x,
270 edgesOutsideXYZ_triangleOutsideW.y || edgesOutside.y,
271 edgesOutsideXYZ_triangleOutsideW.z || edgesOutside.z,
272 all(pointsOutside));
273 }
274
275 return edgesOutsideXYZ_triangleOutsideW;
276}
277
278// Returns 'true' if the edge of the triangle is outside of the frustum.
279// The edges are defined s.t. they are on the opposite side of the point with the given index.
280// 'epsilon' is the (negative) distance to (outside of) the frustum below which we cull the triangle.
281//output packing:
282// x,y,z - one component per triangle edge, true if outside, false otherwise
283bool3 CullTriangleEdgesFrustum(float3 p0, float3 p1, float3 p2, float epsilon, float4 frustumPlanes[6], int numPlanes)
284{
285 return CullFullTriangleAndEdgesFrustum(p0, p1, p2, epsilon, frustumPlanes, numPlanes).xyz;
286}
287
288bool CullTriangleBackFaceView(float3 p0, float3 p1, float3 p2, float epsilon, float3 V, float winding)
289{
290 float3 edge1 = p1 - p0;
291 float3 edge2 = p2 - p0;
292
293 float3 N = cross(edge1, edge2);
294 float NdotV = dot(N, V) * winding;
295
296 // Optimize:
297 // NdotV / (length(N) * length(V)) < Epsilon
298 // NdotV < Epsilon * length(N) * length(V)
299 // NdotV < Epsilon * sqrt(dot(N, N)) * sqrt(dot(V, V))
300 // NdotV < Epsilon * sqrt(dot(N, N) * dot(V, V))
301 return NdotV < epsilon * sqrt(dot(N, N) * dot(V, V));
302}
303
304// Returns 'true' if a triangle defined by 3 vertices is back-facing.
305// 'epsilon' is the (negative) value of dot(N, V) below which we cull the triangle.
306// 'winding' can be used to change the order: pass 1 for (p0 -> p1 -> p2), or -1 for (p0 -> p2 -> p1).
307bool CullTriangleBackFace(float3 p0, float3 p1, float3 p2, float epsilon, float3 viewPos, float winding)
308{
309 float3 V = viewPos - p0;
310 return CullTriangleBackFaceView(p0, p1, p2, epsilon, V, winding);
311}
312
313#endif // UNITY_GEOMETRICTOOLS_INCLUDED