From 315c1fee82a9f556a0fa3e574b7b1698a556ab71 Mon Sep 17 00:00:00 2001
From: Anton Ljungdahl <anton.ljungdahl@fysik.su.se>
Date: Wed, 23 Apr 2025 19:34:09 +0200
Subject: [PATCH] sphere normals with quadratic formula

---
 src/main.cu   |  44 ++++++++++++++++++++++++++++++++------------
 timeBuild.ctm | Bin 964 -> 1124 bytes
 2 files changed, 32 insertions(+), 12 deletions(-)

diff --git a/src/main.cu b/src/main.cu
index 1373299..0e9c1bc 100644
--- a/src/main.cu
+++ b/src/main.cu
@@ -150,21 +150,21 @@ __host__ __device__ function Vec3F32 scale_V3F32(F32 s, Vec3F32 v)
   return out;
 }
 
-__host__ __device__ function F32 dot_V3F32(Vec3F32 a, Vec3F32 b)
+__device__ function F32 dot_V3F32(Vec3F32 a, Vec3F32 b)
 {
   return a.x*b.x + a.y*b.y + a.z*b.z;
 }
 
 
-__device__ function Vec3F32 ray_point_F32(F32 t, RayF32 *ray)
+__device__ function Vec3F32 ray_point_F32(F32 t, RayF32 ray)
 {
-  Vec3F32 out = add_V3F32(ray->origin, scale_V3F32(t, ray->direction));
+  Vec3F32 out = add_V3F32(ray.origin, scale_V3F32(t, ray.direction));
   return out;
 }
 
 __device__ function F32 mag_V3F32(Vec3F32 a)
 {
-  return a.x*a.x + a.y*a.y + a.z*a.z;
+  return dot_V3F32(a, a);
 }
 
 __device__ function F32 norm_V3F32(Vec3F32 a)
@@ -226,27 +226,40 @@ __host__ function void write_buffer_to_ppm(Vec3F32 *buffer, U32 image_width, U32
   fclose(file);
 }
 
-__device__ function U32 hit_sphere(Vec3F32 center, F32 radius, RayF32 r)
+__device__ function F32 hit_sphere(Vec3F32 center, F32 radius, RayF32 r)
 {
   // We take the quadratic formula -b +- sqrt(b*b-4ac) / 2a,
   // and we calculate only the sqrt part. If there is a hit with the sphere we either
   // have two solutions (positive sqrt), one solution (zero sqrt) or no solution (negative sqrt).
   // If we have no solution we have no hit on the sphere centered at center, with the given radius.
-  // No hit -> return 0
-  // Otherwise we do have a hit, return 1
 
   // Compare lines with RTIOW
   // (C-Q)
   Vec3F32 oc = sub_V3F32(center, r.origin);
-  // a = d.d
+  // a = D.D
   F32 a = dot_V3F32(r.direction, r.direction);
-  // b = -2d . (C-Q)
+  // b = -2D . (C-Q)
   F32 b = dot_V3F32(scale_V3F32(-2.0f, r.direction), oc);
   // c = (C-Q) . (C-Q) - r*r
   F32 c = dot_V3F32(oc, oc) - radius*radius;
 
   F32 discriminant = b*b - 4*a*c;
-  U32 out = discriminant >= 0.0f ? 1 : 0;
+
+  // We are actually solving for the parameter t in the expression of a point P(t) that 
+  // intersects the sphere. This is the quadratic problem we get by solving for t in 
+  // (C - P(t)) . (C - P(t)) = r*r, r being the radius and P(t) = tD+Q,
+  // where D is the direction of the ray and Q the origin of the ray.
+  F32 out = 0.0f;
+  if(discriminant < 0.0f)
+  {
+    out = -1.0f;
+  }
+  else
+  {
+    // t = (-b += sqrt(b*b-4ac))/2a, and here we take the smallest solution to get the point 
+    // on the sphere closest to the ray origin.
+    out = (-b - __fsqrt_rn(discriminant))/(2*a); 
+  }
 
   return out;
 }
@@ -284,9 +297,16 @@ __global__ function void cuda_main(Vec3F32 *pixelbuffer, U32 *idxbuffer)
 
     Vec3F32 sphere_center = vec3F32(0.0f, 0.0f, -1.0f);
     F32 sphere_radius = 0.5f;
-    if(hit_sphere(sphere_center, sphere_radius, r)) 
+
+    // t is the parameter of the (closest) sphere-ray intersection point P(t) = tD+Q,
+    // where Q is the ray origin and D the ray direction.
+    F32 t = hit_sphere(sphere_center, sphere_radius, r);
+    if(t > 0.0f)
     {
-      pixel_color = vec3F32(1.0f, 0.0f, 0.0f);
+      Vec3F32 intersection_point = ray_point_F32(t, r);
+      Vec3F32 N = sub_V3F32(intersection_point, sphere_center);
+      N = scale_V3F32(1.0f/sphere_radius, N);
+      pixel_color = scale_V3F32(0.5f, add_V3F32(N, vec3F32(1.0f, 1.0f, 1.0f)));
     } 
     else
     {
diff --git a/timeBuild.ctm b/timeBuild.ctm
index 957b33d17fcd1a737fb2eb01fc03d18301e7087c..20697c4b1c19a07a1849d6bd59a92dc4b798b5ed 100644
GIT binary patch
delta 169
zcmX@Y{)A)05$5{L=U&xCTW&KlGB7Z-urM$<2~>B0_{>251$G98jI``^F#ZD|-~7Sj
z7cljULFzx>cn_0r0?MCQbvdLM$-E^eqDx`?DM0lprJ0jJe31G491MyKuNAX@!W6Uu
P6&N{f{RdN!3^Wh`U4}FZ

delta 7
OcmaFDafE%t5oQ1l{sS)n