Merge pull request #399 from recp/ray

Some missing ray functions

CREDITS

@@ -82,3 +82,17 @@ http://github.com/microsoft/DirectXMath
 17. Pick Matrix
 
 glu project -> project.c
+
+18. Ray sphere intersection
+
+RAY TRACING GEMS
+HIGH-QUALITY AND REAL-TIME RENDERING WITH DXR AND OTHER APIS
+
+CHAPTER 7
+Precision Improvements for Ray/Sphere Intersection
+Eric Haines (1), Johannes Günther (2), and Tomas Akenine-Möller (1)
+  (1) NVIDIA
+  (2) Intel
+
+Wyman, C., and Haines, E. Getting Started with RTX Ray Tracing.
+https://github.com/NVIDIAGameWorks/GettingStartedWithRTXRayTracing

Makefile.am

@@ -216,7 +216,8 @@ cglm_struct_HEADERS = include/cglm/struct/mat4.h \
                       include/cglm/struct/project.h \
                       include/cglm/struct/sphere.h \
                       include/cglm/struct/color.h \
-                      include/cglm/struct/curve.h
+                      include/cglm/struct/curve.h \
+                      include/cglm/struct/ray.h
 
 cglm_struct_clipspacedir=$(includedir)/cglm/struct/clipspace
 cglm_struct_clipspace_HEADERS = include/cglm/struct/clipspace/persp_lh_no.h \

docs/source/ray.rst

@@ -13,6 +13,8 @@ Table of contents (click to go):
 Functions:
 
 1. :c:func:`glm_ray_triangle`
+#. :c:func:`glm_ray_sphere`
+#. :c:func:`glm_ray_at`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -29,3 +31,37 @@ Functions documentation
       | *[in]*       **v2**            third vertex of triangle
       | *[in, out]*  **d**             float pointer to save distance to intersection
       | *[out]*      **intersection**  whether there is intersection
+
+.. c:function:: bool glm_ray_sphere(vec3 origin, vec3 dir, vec4 s, float * __restrict t1, float * __restrict t2)
+
+    ray sphere intersection
+
+    - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+    - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+    - t1 < 0, t2 < 0: no intersection ahead of the ray
+    - the caller can check if the intersection points (t1 and t2) fall within a
+      specific range (for example, tmin < t1, t2 < tmax) to determine if the
+      intersections are within a desired segment of the ray
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     normalized ray direction
+      | *[in]*   **s**       sphere  [center.x, center.y, center.z, radii]
+      | *[out]*  **t1**      near point1 (closer to origin)
+      | *[out]*  **t2**      far point2 (farther from origin)
+
+    Return:
+      | whether there is intersection
+
+.. c:function:: bool  glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point)
+
+    point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+
+    Parameters:
+      | *[in]*   **origin**  ray origin
+      | *[in]*   **dir**     ray direction
+      | *[out]*  **t**       parameter
+      | *[out]*  **point**   point at t
+
+    Return:
+      | point at t

docs/source/vec2.rst

@@ -53,6 +53,8 @@ Functions:
 #. :c:func:`glm_vec2_clamp`
 #. :c:func:`glm_vec2_lerp`
 #. :c:func:`glm_vec2_make`
+#. :c:func:`glm_vec2_reflect`
+#. :c:func:`glm_vec2_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -394,3 +396,25 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec2_reflect(vec2 I, vec2 N, vec2 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

docs/source/vec3.rst

@@ -80,6 +80,9 @@ Functions:
 #. :c:func:`glm_vec3_clamp`
 #. :c:func:`glm_vec3_lerp`
 #. :c:func:`glm_vec3_make`
+#. :c:func:`glm_vec3_faceforward`
+#. :c:func:`glm_vec3_reflect`
+#. :c:func:`glm_vec3_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -512,3 +515,35 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest)
+
+    A vector pointing in the same direction as another
+
+    Parameters:
+      | *[in]*  **N**     vector to orient
+      | *[in]*  **I**     incident vector
+      | *[in]*  **Nref**  reference vector
+      | *[out]* **dest**  destination: oriented vector, pointing away from the surface.
+
+.. c:function:: void glm_vec3_reflect(vec3 I, vec3 N, vec3 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

docs/source/vec4.rst

@@ -60,6 +60,8 @@ Functions:
 #. :c:func:`glm_vec4_lerp`
 #. :c:func:`glm_vec4_cubic`
 #. :c:func:`glm_vec4_make`
+#. :c:func:`glm_vec4_reflect`
+#. :c:func:`glm_vec4_refract`
 
 Functions documentation
 ~~~~~~~~~~~~~~~~~~~~~~~
@@ -424,3 +426,29 @@ Functions documentation
     Parameters:
       | *[in]*  **src**  pointer to an array of floats
       | *[out]* **dest** destination vector
+
+.. c:function:: void glm_vec4_reflect(vec4 I, vec4 N, vec4 dest) 
+
+    Reflection vector using an incident ray and a surface normal
+
+    Parameters:
+      | *[in]*  **I**     incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[out]* **dest**  destination: reflection result
+
+.. c:function:: void glm_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest)
+
+    Refraction vector using entering ray, surface normal and refraction index
+
+    If the angle between the entering ray I and the surface normal N is too 
+    great for a given refraction index, the return value is zero
+
+    this implementation does not explicitly preserve the 'w' component of the
+    incident vector 'I' in the output 'dest', users requiring the preservation of
+    the 'w' component should manually adjust 'dest' after calling this function.
+
+    Parameters:
+      | *[in]*  **I**     *❗️ normalized ❗️* incident vector
+      | *[in]*  **N**     *❗️ normalized ❗️* normal vector
+      | *[in]*  **eta**   ratio of indices of refraction ( η )
+      | *[out]* **dest**  destination: refraction result

include/cglm/call/ray.h

@@ -20,7 +20,19 @@ glmc_ray_triangle(vec3   origin,
                   vec3   v1,
                   vec3   v2,
                   float *d);
-
+
+CGLM_EXPORT
+bool
+glmc_ray_sphere(vec3 origin,
+                vec3 dir,
+                vec4 s,
+                float * __restrict t1,
+                float * __restrict t2);
+
+CGLM_EXPORT
+void
+glmc_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec2.h

@@ -197,6 +197,14 @@ CGLM_EXPORT
 void
 glmc_vec2_make(const float * __restrict src, vec2 dest);
 
+CGLM_EXPORT
+void
+glmc_vec2_reflect(vec2 I, vec2 N, vec2 dest);
+
+CGLM_EXPORT
+void
+glmc_vec2_refract(vec2 I, vec2 N, float eta, vec2 dest);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec3.h

@@ -334,6 +334,18 @@ CGLM_EXPORT
 void
 glmc_vec3_make(const float * __restrict src, vec3 dest);
 
+CGLM_EXPORT
+void
+glmc_vec3_faceforward(vec3 N, vec3 I, vec3 Nref, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_reflect(vec3 I, vec3 N, vec3 dest);
+
+CGLM_EXPORT
+void
+glmc_vec3_refract(vec3 I, vec3 N, float eta, vec3 dest);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/call/vec4.h

@@ -311,6 +311,14 @@ CGLM_EXPORT
 void
 glmc_vec4_make(const float * __restrict src, vec4 dest);
 
+CGLM_EXPORT
+void
+glmc_vec4_reflect(vec4 I, vec4 N, vec4 dest);
+
+CGLM_EXPORT
+void
+glmc_vec4_refract(vec4 I, vec4 N, float eta, vec4 dest);
+
 #ifdef __cplusplus
 }
 #endif

include/cglm/ray.h

@@ -7,12 +7,18 @@
 
 /*
  Functions:
-   CGLM_INLINE bool glm_line_triangle_intersect(vec3   origin,
-                                                vec3   direction,
-                                                vec3   v0,
-                                                vec3   v1,
-                                                vec3   v2,
-                                                float *d);
+   CGLM_INLINE bool glm_ray_triangle(vec3   origin,
+                                     vec3   direction,
+                                     vec3   v0,
+                                     vec3   v1,
+                                     vec3   v2,
+                                     float *d);
+ CGLM_INLINE bool glm_ray_sphere(vec3 origin,
+                                 vec3 dir,
+                                 vec4 s,
+                                 float * __restrict t1,
+                                 float * __restrict t2)
+ CGLM_INLINE void glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point);
 */
 
 #ifndef cglm_ray_h
@@ -31,7 +37,6 @@
  * @param[in, out] d         distance to intersection
  * @return whether there is intersection
  */
-
 CGLM_INLINE
 bool
 glm_ray_triangle(vec3   origin,
@@ -74,4 +79,93 @@ glm_ray_triangle(vec3   origin,
   return dist > epsilon;
 }
 
+/*!
+ * @brief ray sphere intersection
+ *
+ * - t1 > 0, t2 > 0: ray intersects the sphere at t1 and t2 both ahead of the origin
+ * - t1 < 0, t2 > 0: ray starts inside the sphere, exits at t2
+ * - t1 < 0, t2 < 0: no intersection ahead of the ray
+ * - the caller can check if the intersection points (t1 and t2) fall within a
+ *   specific range (for example, tmin < t1, t2 < tmax) to determine if the
+ *   intersections are within a desired segment of the ray
+ *
+ * @param[in]  origin ray origin
+ * @param[out] dir    normalized ray direction
+ * @param[in]  s      sphere  [center.x, center.y, center.z, radii]
+ * @param[in]  t1     near point1 (closer to origin)
+ * @param[in]  t2     far point2 (farther from origin)
+ */
+CGLM_INLINE
+bool 
+glm_ray_sphere(vec3 origin,
+               vec3 dir,
+               vec4 s,
+               float * __restrict t1,
+               float * __restrict t2) {
+  vec3  dp;
+  float r2, ddp, dpp, dscr, q, tmp, _t1, _t2;
+
+  /* ensure dir is normalized */
+  glm_vec3_sub(s, origin, dp);
+
+  ddp = glm_vec3_dot(dir, dp);
+  dpp = glm_vec3_norm2(dp);
+
+  /* compute the remedy term for numerical stability */
+  glm_vec3_mulsubs(dir, ddp, dp); /* dp: remedy term */
+
+  r2   = s[3] * s[3];
+  dscr = r2 - glm_vec3_norm2(dp);
+
+  if (dscr < 0.0f) {
+    /* no intersection */
+    return false;
+  }
+
+  dscr = sqrtf(dscr);
+  q    = (ddp >= 0.0f) ? (ddp + dscr) : (ddp - dscr);
+
+  /*
+     include Press, William H., Saul A. Teukolsky,
+     William T. Vetterling, and Brian P. Flannery,
+     "Numerical Recipes in C," Cambridge University Press, 1992.
+   */
+  _t1 = q;
+  _t2 = (dpp - r2) / q;
+
+  /* adjust t1 and t2 to ensure t1 is the closer intersection */
+  if (_t1 > _t2) {
+    tmp = _t1;
+    _t1 = _t2;
+    _t2 = tmp;
+  }
+
+  *t1 = _t1;
+  *t2 = _t2;
+
+  /* check if the closest intersection (t1) is behind the ray's origin */
+  if (_t1 < 0.0f && _t2 < 0.0f) {
+    /* both intersections are behind the ray, no visible intersection */
+    return false;
+  }
+
+  return true;
+}
+
+/*!
+ * @brief point using t by 𝐏(𝑡)=𝐀+𝑡𝐛
+ *
+ * @param[in]  orig  origin of ray
+ * @param[in]  dir   direction of ray
+ * @param[in]  t     parameter
+ * @param[out] point point at t
+ */
+CGLM_INLINE
+void
+glm_ray_at(vec3 orig, vec3 dir, float t, vec3 point) {
+  vec3 dst;
+  glm_vec3_scale(dir, t, dst);
+  glm_vec3_add(orig, dst, point);
+}
+
 #endif

include/cglm/struct.h

@@ -41,6 +41,7 @@ extern "C" {
 #include "struct/sphere.h"
 #include "struct/curve.h"
 #include "struct/affine2d.h"
+#include "struct/ray.h"
 
 #ifdef __cplusplus
 }