sphere normals with quadratic formula

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
     {