Ray tracer in c running slowly

Question

I was learning C the past few days and decided to make a mini project to get a better feel for it following the Ray Tracing In a Weekend Series. I got it working however I feel like it is pretty slow as the main image took 3 hours to run. I couldn't really find a good way to profile my code(On Mac M1) and I'm pretty sure that there any many things that are making the code super slow which I am unaware of. Here is the code:

main.c

#include "utils/vec3.h"
#include "utils/ray.h"
#include "utils/color.h"
#include "utils/sphere.h"
#include "utils/hittable.h"
#include "utils/hittable_list.h"
#include "utils/camera.h"
#include "utils/utility.h"
#include "utils/hittable_types.h"
#include "utils/material.h"
#include "utils/material_types.h"

#include <stdio.h>
#include <time.h>
#include <stdbool.h>
#include <math.h>

typedef struct pixel{
    int r;
    int g;
    int b;

} pixel;
color ray_color(const ray* r, const hittable_list *world, int depth){
    hit_record rec;

    if (depth <= 0) {
        return col(0, 0, 0);
    }
    if (hit(world, r, 0.0001, 10000, &rec)){
        ray scattered;
        color attenuation;
        if (scatter(rec.mat_ptr, r, &rec, &attenuation, &scattered)) {
            return mul(ray_color(&scattered, world, depth-1), attenuation);
        }
        return (color){0, 0, 0};
    }
    vec3 unit_direction = unit_vector(r->direction);
    double t = 0.5 * (unit_direction.y + 1.0);
    return add(scale(col(1.0, 1.0, 1.0), (1.0 - t)), scale(col(0.5, 0.7, 1.0), t));
}

void random_scene(hittable_list *world) {
    world->current_index = 0;
    material *ground_material = initialize_lambertian(col(0.5, 0.5, 0.5));
    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){p3(0.0, -1000, 0), 1000, ground_material}});
    for(int a = -11; a < 11; a++) {
        for(int b = -11; b < 11; b++) {
            double choose_mat = random_double();
            point3 center = p3(a + 0.9 * random_double(), 0.2, b + 0.9 * random_double());
            if (length(sub(center, p3(4, 0.2, 0)))> 0.9) {
                material *sphere_material;
                if (choose_mat < 0.8) {
                    color albedo = mul(random_vec(), random_vec());
                    sphere_material = initialize_lambertian(albedo);
                    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){center, 0.2, sphere_material}});
                }
                else if (choose_mat < 0.95){
                    color albedo = random_vec_mm(0.5, 0.1);
                    double fuzz = random_double_mm(0, 0.5);
                    sphere_material = initialize_metal(albedo, fuzz);
                    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){center, 0.2, sphere_material}});
                }
                else {
                    sphere_material = initialize_dielectric(1.5);
                    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){center, 0.2, sphere_material}});               
                }
            }
        }
    }
    
    material *material1 = initialize_dielectric(1.5);
    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){p3(0, 1, 0), 1.0, material1}});

    material *material2 = initialize_lambertian(col(0.4, 0.2, 0.1));
    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){p3(-4, 1, 0), 1.0, material2}});

    material *material3 = initialize_metal(col(0.7, 0.6, 0.5), 0.0);
    add_obj(world, (hittable) {.type = SPHERE, .s = (sphere){p3(4, 1, 0), 1.0, material3}});

}
int main() {

    // Image

    const float aspect_ratio = 3.0 / 2.0;
    const int image_width = 1200;
    const int image_height = (int) (image_width / aspect_ratio);
    const int samples_per_pixel = 32;
    const int max_depth = 8;

    // World
    hittable_list world;
    random_scene(&world);

    // Camera

    camera cam;

    point3 lookfrom = p3(13, 2, 3);
    point3 lookat = p3(0, 0, 0);
    vec3 vup = v3(0, 1, 0);
    double dist_to_focus = 10.0;
    double aperture = 0.1;
    setup_camera(&cam, lookfrom, lookat, vup, 20, aspect_ratio, aperture, dist_to_focus);

    clock_t begin = clock();
    printf("P3\n%d %d\n255\n", image_width, image_height);

    for(int j = image_height - 1; j >= 0; j--) {
        fprintf(stderr, "%d iterations remaining!\n", j);
        for(int i = 0; i < image_width; i++) {
            color pixel_color = col(0, 0, 0);
            for (int s = 0; s < samples_per_pixel; s++) {
                double u = (double) (((double)(i) + random_double()) / (image_width - 1));
                double v = (double) (((double)(j) + random_double()) / (image_height - 1));
                ray r = get_ray(&cam, u, v);
                iadd(&pixel_color, ray_color(&r, &world, max_depth));
            }
            write_color(pixel_color, samples_per_pixel);
        }
    }
    clock_t end = clock();
    fprintf(stderr, "Time spent: %f seconds!\n", (double)(end - begin) / CLOCKS_PER_SEC);
}

camera.h

#ifndef CAMERA_H
#define CAMERA_H

#include "vec3.h"
#include "ray.h"
typedef struct {
    point3 lookfrom;
    point3 lookat;
    vec3 vup;
    double vfov;
    float aspect_ratio;
    double aperture;
    double focus_dist;
    float viewport_height;
    float viewport_width;
    double len_radius;
    vec3 u, v, w;

    point3 origin;
    vec3 horizontal;
    vec3 vertical;
    vec3 lower_left_corner;
} camera;

void setup_camera(camera *cam, point3 lookfrom_t, point3 lookat_t, vec3 vup_t, double vfov_t, double aspect_ratio_t, double aperture_t, double focus_dist_t) {
    cam->lookfrom = lookfrom_t;
    cam->lookat = lookat_t;
    cam->vup = vup_t;
    cam->vfov = vfov_t;
    cam->aspect_ratio = aspect_ratio_t;
    cam->aperture = aperture_t;
    cam->focus_dist = focus_dist_t;
    double theta = degrees_to_radians(cam->vfov);
    double h = tan(theta / 2);
    cam->viewport_height = 2.0 * h;
    cam->viewport_width = cam->aspect_ratio * cam->viewport_height;

    cam->w = unit_vector(sub(cam->lookfrom, cam->lookat));
    cam->u = unit_vector(cross(cam->vup, cam->w));
    cam->v = cross(cam->w, cam->u);

    cam->origin = cam->lookfrom;
    cam->horizontal = scale(cam->u, cam->viewport_width * cam->focus_dist);
    cam->vertical = scale(cam->v, cam->viewport_height * cam->focus_dist);
    cam->lower_left_corner = sub(sub(cam->origin, mdiv(cam->horizontal, 2)), mdiv(cam->vertical, 2));
    cam->lower_left_corner = sub(cam->lower_left_corner, scale(cam->w, cam->focus_dist));
    cam->len_radius = cam->aperture / 2;
}

ray get_ray(camera *cam, double s, double t) {
    vec3 rd = scale(random_in_unit_disk(), cam->len_radius);
    vec3 offset = add(scale(cam->u, rd.x), scale(cam->v, rd.y));
    vec3 dir = add(add(cam->lower_left_corner, scale(cam->horizontal, s)), scale(cam->vertical, t));
    ray r = {.origin = add(cam->origin, offset), .direction = sub(sub(dir, cam->origin), offset)};
    return r;

}
#endif

color.h

#ifndef COLOR_H
#define COLOR_H

#include <stdio.h>
#include <math.h>
#include "vec3.h"

static inline double clamp(double x, double min, double max) {
    if (x < min) return min;
    if (x > max) return max;
    return x;
}
void write_color(color pixel_color, int samples_per_pixel){

    double r = pixel_color.x;
    double g = pixel_color.y;
    double b = pixel_color.z;

    double scale = (1.0 / samples_per_pixel);
    r = sqrt(scale * r);
    g = sqrt(scale * g);
    b = sqrt(scale * b);

    printf("%d %d %d\n", (int) (256 * clamp(r, 0.0, 0.999)), (int) (256 * clamp(g, 0.0, 0.999)), (int) (256 * clamp(b, 0.0, 0.999)));
}

#endif

hittable.c

#include "hittable.h"


void set_face_normal(hit_record *rec, const ray *r, vec3 outward_normal) {
    rec->front_face = (dot(r->direction, outward_normal) < 0);
    rec->normal = rec->front_face ? outward_normal : negate(outward_normal);
}

hittable.h

#ifndef HITTABLE_H
#define HITTABLE_H

#include "ray.h"
#include "material.h"
#include <stdbool.h>

struct material;
typedef struct hit_record{
    point3 p;
    vec3 normal;
    struct material* mat_ptr;
    double t;
    bool front_face;

} hit_record;

void set_face_normal(hit_record *rec, const ray *r, vec3 outward_normal);


#endif

hittable_list.c

#include "hittable.h"
#include "sphere.h"
#include "hittable_types.h"
#include "hittable_list.h"
#include <string.h>

bool hit_scene(const hittable obj, const ray* r, double t_min, double t_max, hit_record *rec) {
    switch(obj.type) {
        case SPHERE:
            return hit_sphere(obj.s, r, t_min, t_max, rec);
    }
    return false;
}
static inline void clear(hittable_list *world){
    memset(world->objects, 0, sizeof(world->objects));
}

void add_obj(hittable_list *world, hittable object) {
    world->objects[world->current_index++] = object;
}

bool hit(const hittable_list *world, const ray *r, double t_min, double t_max, hit_record *rec){
    hit_record temp_rec;
    bool hit_anything = false;
    double closest_so_far = t_max;

    for(int i = 0; i < world->current_index; i++) {
        //if(world->objects[i].hit(world->objects[i], r, t_min, closest_so_far, &temp_rec)){
        if(hit_scene(world->objects[i], r, t_min, closest_so_far, &temp_rec)){
            hit_anything = true;
            closest_so_far = temp_rec.t;
            *rec = temp_rec;
        }
    }
    return hit_anything;
}

hittable_list.h

#ifndef HITTABLE_LIST_H
#define HITTABLE_LIST_H

#include "hittable.h"
#include "sphere.h"
#include "hittable_types.h"

#include <string.h>

typedef struct {
    int type;
    union {
        sphere s;
    };
} hittable;

typedef struct {
    int current_index;
    hittable objects[500];
} hittable_list;


bool hit_scene(const hittable obj, const ray* r, double t_min, double t_max, hit_record *rec);

static inline void clear(hittable_list *world);
void add_obj(hittable_list *world, hittable object);
bool hit(const hittable_list *world, const ray *r, double t_min, double t_max, hit_record *rec);
#endif

hittable_types.h

#ifndef HITTABLE_TYPES_H
#define HITTABLE_TYPES_H
enum {SPHERE = 0} types;

#endif

material.c

#include "material.h"

static double reflectance(double cosine, double ref_idx) {
    double r0 = (1 - ref_idx) / (1 + ref_idx);
    r0 = r0 * r0;
    return r0 + (1-r0)*pow((1-cosine), 5);
}
bool lambertian_scatter(lambertian lam, const ray* r_in, const struct hit_record* rec, color* attenuation, ray* scattered) {
    vec3 scatter_direction = add(rec->normal, random_unit_vector());
    if (near_zero(scatter_direction)) {
        scatter_direction = rec->normal;
    }
    *scattered = (ray){rec->p, scatter_direction};
    *attenuation = lam.albedo;
    return true;
}

bool metal_scatter(metal m, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered) {
    vec3 reflected = reflect(unit_vector(r_in->direction), rec->normal);
    *scattered = (ray){rec->p, add(reflected, scale(random_in_unit_sphere(), m.fuzz))};
    *attenuation = m.albedo;
    return (dot(scattered->direction, rec->normal) > 0);
}

bool dielectric_scatter(dielectric d, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered) {
    *attenuation = (color) {1.0, 1.0, 1.0};
    double refraction_ratio = rec->front_face ? (1.0 / d.ir) : d.ir;

    vec3 unit_direction = unit_vector(r_in->direction);
    double cos_theta = fmin(dot(negate(unit_direction), rec->normal), 1.0);
    double sin_theta = sqrt(1.0 - cos_theta * cos_theta);

    bool cannot_refract = refraction_ratio * sin_theta > 1.0;
    vec3 direction;

    if (cannot_refract || reflectance(cos_theta, refraction_ratio) > random_double()) {
        direction = reflect(unit_direction, rec->normal);
    }
    else {
        direction = refract(unit_direction, rec->normal, refraction_ratio);
    }
    *scattered = (ray) {rec->p, direction};
    return true;
}
bool scatter(material *mat, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered) {
    switch(mat->type) {
        case LAMBERTIAN:
            return lambertian_scatter(mat->lam, r_in, rec, attenuation, scattered);

        case METAL:
            return metal_scatter(mat->m, r_in, rec, attenuation, scattered);

        case DIELECTRIC:
            return dielectric_scatter(mat->d, r_in, rec, attenuation, scattered);

    }
    return false;
}

material *initialize_lambertian(color albedo) {
    material *mat = (material *)malloc(sizeof(material));
    mat->type = LAMBERTIAN;
    mat->lam = (lambertian) {albedo};
    return mat;
}
material *initialize_metal(color albedo, double fuzz) {
    material *mat = (material *) malloc(sizeof(material));
    mat->type = METAL;
    mat->m = (metal) {albedo, fuzz};
    return mat;
}
material *initialize_dielectric(double ir) {
    material *mat = (material *) malloc(sizeof(material));
    mat->type = DIELECTRIC;
    mat->d = (dielectric) {ir};
    return mat;
}

material.h

#ifndef MATERIAL_H
#define MATERIAL_H
#include "vec3.h"
#include "ray.h"
#include "hittable.h"
#include "material_types.h"

#include <math.h>

struct hit_record;
typedef struct {
    color albedo;

} lambertian;

typedef struct {
    color albedo;
    double fuzz;
} metal;

typedef struct {
    double ir;
} dielectric;

typedef struct material{
    int type;
    union {
        lambertian lam;
        metal m;
        dielectric d;
    };
} material;

bool lambertian_scatter(lambertian lam, const ray* r_in, const struct hit_record* rec, color *attenuation, ray* scattered);
bool metal_scatter(metal m, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered);
bool dielectric_scatter(dielectric d, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered);
bool scatter(material *mat, const ray* r_in, const struct hit_record *rec, color *attenuation, ray* scattered);
material *initialize_lambertian(color albedo);
material *initialize_metal(color albedo, double fuzz);
material *initialize_dielectric(double ir);

#endif

material_types.h

#ifndef MATERIAL_TYPES_H
#define MATERIAL_TYPES_H
enum {LAMBERTIAN = 0, METAL = 1, DIELECTRIC = 2} mat_types;

#endif

ray.c

#include "ray.h"

point3 rayat(const ray *r, const double t) {
    return add(r->origin, scale(r->direction, t));
}

ray.h

#ifndef RAY_H
#define RAY_H

#include "vec3.h"

typedef struct{
    point3 origin;
    vec3 direction;
} ray;

point3 rayat(const ray *r, const double t);

#endif

sphere.c

#include "hittable.h"
#include "vec3.h"
#include "material.h"
#include "sphere.h"

bool hit_sphere(const sphere s, const ray* r, double t_min, double t_max, hit_record *rec) {
    vec3 oc = sub(r->origin, s.center);
    double a = length_squared(r->direction);
    double half_b = dot(oc, r->direction);
    double c = length_squared(oc) - s.radius * s.radius;
    double discriminant = half_b * half_b - a * c;
    if (discriminant < 0) {
        return false;
    }

    double sqrtd = sqrt(discriminant);
    double root = (-half_b - sqrtd) / a;
    if (root < t_min || t_max < root) {
        root = (-half_b + sqrtd) / a;
        if (root < t_min || t_max < root) {
            return false;
        }
    }


    rec->t = root;
    rec->p = rayat(r, rec->t);
    vec3 outward_normal = mdiv(sub(rec->p, s.center), s.radius);
    set_face_normal(rec, r, outward_normal);
    rec->mat_ptr = s.mat_ptr;

    return true;
}

sphere.h

#ifndef SPHERE_H
#define SPHERE_H

#include "hittable.h"
#include "vec3.h"
#include "material.h"

typedef struct s sphere;
typedef struct s{
    point3 center;
    double radius;
    material *mat_ptr;
    //bool (*hit)(sphere, const ray*, double, double, hit_record*);

} sphere;

bool hit_sphere(const sphere s, const ray* r, double t_min, double t_max, hit_record *rec);
#endif

utility.c

#include "utility.h"


double random_double()
{
    return (double)rand() / (double)RAND_MAX ;
}

double random_double_mm(double min, double max) {
    return min + (max-min)*random_double();
}

double degrees_to_radians(double degrees) {
    return degrees * (M_PI / 180.0);
}

utility.h

#ifndef UTILITY_H
#define UTILITY_H
#include <stdlib.h>
#include <math.h>

double random_double();

double random_double_mm(double min, double max);

double degrees_to_radians(double degrees);
#endif

vec3.h

#ifndef VEC3_H
#define VEC3_H

#include "utility.h"


#include <math.h>
#include <stdio.h>
#include <stdbool.h>

struct vec{
    double x;
    double y;
    double z;
};

typedef struct vec vec3;
typedef struct vec color;
typedef struct vec point3;

static inline struct vec v3(double x, double y, double z){
    vec3 vector = {x, y, z};
    return vector;
}

static inline struct vec col(double x, double y, double z){
    color vector = {x, y, z};
    return vector;
}
static inline struct vec p3(double x, double y, double z){
    point3 vector = {x, y, z};
    return vector;
}
static inline struct vec negate(const struct vec v){
    return v3(-v.x, -v.y, -v.z);
}

static inline void iadd(struct vec *v, const struct vec u){
    v->x += u.x;
    v->y += u.y;
    v->z += u.z;
}

static inline struct vec add(const struct vec v, const struct vec u){
    return v3(v.x + u.x, v.y + u.y, v.z + u.z);
}

static inline struct vec sub(const struct vec v, const struct vec u){
    return v3(v.x - u.x, v.y - u.y, v.z - u.z);
}

static inline void iscale(struct vec*v, const double t){
    v->x *= t;
    v->y *= t;
    v->z *= t;

}

static inline struct vec scale(const struct vec v, const double t) {
    return v3(v.x * t, v.y * t, v.z *t);
}

static inline void imul(struct vec *v, const struct vec u) {
    v->x *= u.x;
    v->y *= u.y;
    v->z *= u.z;
}

static inline struct vec mul(const struct vec v, const struct vec u){
    return v3(v.x * u.x, v.y * u.y, v.z * u.z);
}


static inline void idiv(struct vec *v, const double t) {
    v->x /= t;
    v->y /= t;
    v->z /= t;
}

static inline struct vec mdiv(const struct vec v, const double t){
    return v3(v.x * (1 / t), v.y * (1 / t), v.z * (1 / t));
}

static inline double length_squared(const struct vec v) {
    return v.x * v.x + v.y * v.y + v.z * v.z;
}

static inline double length(const struct vec v) {
    return sqrt(length_squared(v));
}

//UTILITY

static inline void print_vec(const struct vec v){
    fprintf(stderr, "(%f, %f, %f)\n", v.x, v.y, v.z);
}

static inline double dot(const struct vec v, const struct vec u) {
    return v.x * u.x + v.y * u.y + v.z * u.z;
}

static inline struct vec cross(const struct vec v, const struct vec u) {
    return v3(v.y * u.z - v.z * u.y,
                v.z * u.x - v.x * u.z,
                v.x * u.y - v.y * u.x);
}

static inline struct vec unit_vector(struct vec v){
    return mdiv(v, length(v));
}

static inline struct vec random_vec() {
        return v3(random_double(), random_double(), random_double());
}

static inline struct vec random_vec_mm(double min, double max) {
        return v3(random_double_mm(min, max), random_double_mm(min, max), random_double_mm(min, max));
}

static inline struct vec random_in_unit_sphere() {
    while(true) {
        struct vec p = random_vec_mm(-1, 1);
        if(length_squared(p) >= 1) continue;
        return p;
    }

}
static inline struct vec random_unit_vector() {
    return unit_vector(random_vec());
}

static inline struct vec random_in_unit_disk() {
    while(true) {
        struct vec p = v3(random_double_mm(-1, 1), random_double_mm(-1, 1), 0);
        if(length_squared(p) >= 1) continue;
        return p;
    }
}

static inline bool near_zero(struct vec v) {
    double s = 1e-8;
    return (fabs(v.x) < s) && (fabs(v.y) < s) && (fabs(v.z) < s);
}

static inline struct vec reflect(struct vec v, struct vec n) {
    return sub(v, scale(n, dot(v, n) * 2));
}

static inline struct vec refract(const struct vec uv, const struct vec n, double etai_over_etat) {
    double cost_theta = fmin(dot(negate(uv), n), 1);
    struct vec r_out_perp = scale(add(uv, scale(n, cost_theta)), etai_over_etat);
    struct vec r_out_parallel = scale(n, -sqrt(fabs(1 - length_squared(r_out_perp))));
    return add(r_out_perp, r_out_parallel);
}
#endif

The whole code can be found in this repo if that is easier for anyone: https://github.com/Lumijek/c-raytracer

Answer 1

Compile with optimizations enabled

I strongly suspect you have been compiling your code without enabling the compiler's optimization flags. Consider compiling with:

-Ofast -ffast-math -flto

This enables lots of optimizations, including some "unsafe" ones that deal with floating point edge cases that generally should not be an issue for a raytracer. -flto is especially important, as it allows functions defined in one .c file to be inlined in another, which avoids function call overhead and possibly allows even more optimizations.

On some platforms -march=native might also help, but this unfortunately is not supported by the default compiler on macOS on an M1.

With this, it only takes 51 seconds for your program to run on my M1.

I also tried changing all doubles to floats, which doesn't seem to affect the quality of the output noticably, but this has very little effect on the performance. It might start to matter if you have more objects in your scene though, as it would allow more of them to fit in the CPU cache.

Parallelize your code

You can trivially parallelize your code by starting one thread per CPU core, and spread rendering individual lines of the image over multiple threads. The only thing to worry about is ensuring the output is in the correct order. I suggest allocating enough memory for the whole image, and after all threads have finished, just let the main thread convert it to PPM format.

On a M1, you have 8 cores, but beware that some are faster than others, so just splitting the work in 8 equal pieces is too naive to make optimal use of all cores.

How fast should it run?

It might be instructive to look at how much work your program has to do, and then check how much time that would take on a given processor. You are rendering a 1200 x 800 pixel image using 32 samples per pixel, and there are about 500 objects in the scene that you test against for every ray. So it takes about 15 billion operations, each of which still takes lots of CPU cycles to perform. Consider a single M1 core running at 3.2 GHz, then if it would only take one CPU cycle on average to test a ray against an object, your program would take about 5 seconds to run. With optimizations enabled, your code seems to take about 11 CPU cycles on average. That's pretty great!

I think you cannot squeeze much more out of it without changing the algorithms and datastructures. In particular, testing against every object is wasteful, and doesn't scale to larger scenes. A typical solution is to group objects together, so you can test if a ray intersects the bounding box of a group, and if not you don't have to test against the individual objects, otherwise you have to. You can nest groups in each other in various ways. A possible data structure you could use to organize your objects for efficient traversal is an octree.

Incorrect enum declarations

When you declare an enum like so:

enum {SPHERE = 0} types;

You are actually declaring a variable types which has an anonymous enum type. If multiple .c files #include the file with this declaration, then it can create a conflict when linking. It seems Xcode allows it, but GCC doesn't like it. However, I think you meant to write this instead, which avoids the issue:

enum types {SPHERE = 0};

You can also just omit types altogether, as you never use it anywhere else. The same goes for the other enum.

Don't malloc() materials

It is a bit weird to see lots of objects being created and put into statically allocated arrays, but only for three materials you are using malloc(). I would either malloc() most objects dynamically, and avoid any statically sized arrays, or I would avoid malloc() altogether, and just store the material in sphere by value.

Unnecessary casts

In C, you should not have to cast the pointer returned from malloc() to another type. This is only necessary for C++, and there you should not be using malloc() anyway.

Inside an initializer list you don't have to "cast" when you want to initialize substructures. So you can write:

add_obj(world, (hittable) {
    .type = SPHERE,
    .s = {p3(0.0, -1000, 0), 1000, &ground_material} // no (sphere) necessary
});

It's also not necessary to cast things to double that are already double, like:

double u = (double) (((double)(i) + random_double()) / (image_width - 1));

This can be written as:

double u = (i + random_double()) / (image_width - 1);

It can sometimes be a bit tricky to see when you need to cast and when not, so adding casts "just to be safe" might seem like a good strategy, but it also adds noise to the code.

Code style

Overall the code style is OK, there are just some minor formatting issues, like spaces missing here and there, random empty lines in some places, sometimes writing type* foo, other times type *foo. Consider running a code formatting tool like ClangFormat on your code, this will format everything for you. It doesn't really matter what style you use, as long as it is consistent.

Answer 2

I feel like it is pretty slow ...

Profile

Profile code and get a real assessment of where code is spending its time to get a short-list of candidate optimization re-writes.

float vs. double

Consider float instead of double for color math. I very much doubt the reduced precision of float will render a visible degradation over double, yet potentially 2x+ faster.

Be sure to use float objects, constants and functions like float r;, 1.0f / samples_per_pixel, sqrtf() and ((float) M_PI).

Use integer math

write_color() could do the x > max test of clamp() with int math. x < min is never true as r,g,b are all non-negative. Something like:

//                       r*256*256 
uint32fast_t r1 = ((uint32fast_t)r << 16) >> samples_per_pixel_shift /* 5 */;
unsigned r2 = sqrtf(r1);
if (r2 >= 256) r2 = 255;
printf("%d", r2);

I'd even consider re-writing to use no FP math and roll my own unsigned isqrt32(uint32_t).