Sin descripción
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

GeometricTools.hlsl 12KB

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