/*
 * MINIMAL 3D RENDERER
 * ===================
 * 
 * This renders a spinning cube in 3D. The core idea is transforming points
 * from 3D space onto a 2D screen, then filling triangles with pixels.
 * 
 * The rendering pipeline transforms each vertex through several spaces:
 * - Object space (the cube's original coordinates)
 * - World space (after moving/rotating the cube)  
 * - Camera space (relative to where the camera is looking)
 * - Screen space (final 2D pixel coordinates)
 * 
 * Backface culling skips triangles facing away from the camera, which both
 * optimizes performance and prevents drawing the hidden back faces of the cube.
 */

#include 
#include 
#include 
#include 

/* Basic 3D vector representing a point in space */
typedef struct { float x, y, z; } Vec3;

/* Triangle defined by three vertices and a color */
typedef struct {
    Vec3 v[3];
    unsigned int color;
} Triangle;

/* 
 * Cross product computes a vector perpendicular to two input vectors using
 * the determinant formula. Given vectors a=(ax,ay,az) and b=(bx,by,bz), the
 * cross product a×b produces a vector whose components are:
 * 
 * x = ay*bz - az*by
 * y = az*bx - ax*bz  
 * z = ax*by - ay*bx
 * 
 * Geometrically, the magnitude |a×b| equals the area of the parallelogram
 * formed by a and b. The direction follows the right-hand rule: if you curl
 * your right hand from a toward b, your thumb points in the direction of a×b.
 * 
 * For a triangle with vertices v0,v1,v2, computing (v1-v0)×(v2-v0) gives the
 * surface normal vector that points outward from the front face. This is
 * fundamental for lighting calculations and determining which side of a surface
 * we're looking at.
 */
Vec3 cross(Vec3 a, Vec3 b) {
    return (Vec3){a.y*b.z - a.z*b.y, a.z*b.x - a.x*b.z, a.x*b.y - a.y*b.x};
}

/* 
 * Dot product computes the scalar a·b = ax*bx + ay*by + az*bz. This measures
 * the projection of one vector onto another. Mathematically, a·b = |a||b|cos(θ)
 * where θ is the angle between the vectors and |a| denotes the magnitude
 * (length) of vector a computed as sqrt(ax²+ay²+az²).
 * 
 * From the cosine relationship, we can determine angles without computing them:
 * - If a·b > 0: angle is acute (< 90°), vectors point somewhat together
 * - If a·b = 0: angle is exactly 90°, vectors are perpendicular
 * - If a·b < 0: angle is obtuse (> 90°), vectors point somewhat apart
 * 
 * For backface culling, we compute normal·view_direction. If this is negative,
 * the surface normal points away from the camera (angle > 90°), meaning we're
 * looking at the back of the triangle. The mathematical property dot(a,b) < 0
 * precisely captures "pointing in opposite directions."
 */
float dot(Vec3 a, Vec3 b) {
    return a.x*b.x + a.y*b.y + a.z*b.z;
}

/* Vector subtraction gives the direction from b to a */
Vec3 sub(Vec3 a, Vec3 b) {
    return (Vec3){a.x-b.x, a.y-b.y, a.z-b.z};
}

/*
 * Rotation around the Y axis applies the 2D rotation matrix to the X-Z plane
 * while leaving Y unchanged. The standard 2D rotation matrix is:
 * 
 * [cos(θ)  -sin(θ)] [x]   [x*cos(θ) - z*sin(θ)]
 * [sin(θ)   cos(θ)] [z] = [x*sin(θ) + z*cos(θ)]
 * 
 * In 3D, this becomes:
 * x' = x*cos(θ) + z*sin(θ)
 * y' = y
 * z' = -x*sin(θ) + z*cos(θ)
 * 
 * The trigonometric functions encode circular motion. As θ increases from 0
 * to 2π, the point (x,z) traces a circle. At θ=0, cos(0)=1 and sin(0)=0, so
 * x'=x and z'=z (no rotation). At θ=π/2 (90°), cos(π/2)=0 and sin(π/2)=1,
 * so x'=z and z'=-x (90° rotation). The negative sign in the z' equation
 * determines the rotation direction (clockwise vs counter-clockwise when
 * viewed from above).
 */
Vec3 rotate_y(Vec3 v, float angle) {
    float c = cosf(angle), s = sinf(angle);
    return (Vec3){v.x*c + v.z*s, v.y, -v.x*s + v.z*c};
}

/*
 * Rotation around the X axis applies 2D rotation to the Y-Z plane while X
 * remains constant. The transformation matrix operates on the Y and Z
 * components:
 * 
 * x' = x
 * y' = y*cos(θ) - z*sin(θ)
 * z' = y*sin(θ) + z*cos(θ)
 * 
 * This rotates points around the X axis like a wheel spinning. Combined with
 * Y-axis rotation, we can achieve arbitrary 3D orientations through Euler
 * angles, though the order of rotations matters due to the non-commutative
 * property of matrix multiplication (R_x·R_y ≠ R_y·R_x).
 */
Vec3 rotate_x(Vec3 v, float angle) {
    float c = cosf(angle), s = sinf(angle);
    return (Vec3){v.x, v.y*c - v.z*s, v.y*s + v.z*c};
}

/*
 * Perspective projection transforms 3D points to 2D screen coordinates using
 * the pinhole camera model. A point (x,y,z) in camera space projects onto a
 * 2D image plane at distance f (focal length) from the origin.
 * 
 * By similar triangles, a point at distance z projects to position:
 * x_screen = f * (x/z)
 * y_screen = f * (y/z)
 * 
 * The division by z creates the perspective effect: larger z (farther away)
 * produces smaller screen coordinates, making distant objects appear smaller.
 * This is realistic (how cameras and eyes work) but means the cube appears
 * distorted - the front face looks larger than the back face because it's
 * closer to the camera.
 * 
 * The focal length f controls the field of view - larger f zooms in (telephoto
 * lens), smaller f zooms out (wide-angle lens). At f=600 pixels with an 800px
 * wide screen, we get approximately a 53° horizontal field of view, which is
 * moderate (not wide-angle, not telephoto).
 * 
 * We add 8 to z before dividing to move the cube away from the camera (which
 * sits at z=0). We then add 400 and 300 to x and y to translate the origin
 * from the top-left corner to the center of our 800×600 screen. Without this
 * translation, the cube would render in the top-left corner since screen
 * coordinates typically have (0,0) there rather than at the center.
 * 
 * Note: For a geometrically perfect cube (parallel edges, equal face sizes),
 * you would need orthographic projection where x_screen = x and y_screen = y
 * without any division by z. Orthographic projection is used in CAD and
 * technical drawings, but it lacks depth perception since all objects appear
 * the same size regardless of distance.
 */
Vec3 project(Vec3 v, float focal_length) {
    float scale = focal_length / (v.z + 8.0f);
    return (Vec3){v.x * scale + 400, v.y * scale + 300, v.z};
}

/*
 * The edge function computes the signed area of the triangle formed by points
 * a, b, and p using the 2D cross product (the z-component of a 3D cross product
 * where z=0). Mathematically:
 * 
 * edge(a,b,p) = (p-a) × (b-a) = (px-ax)(by-ay) - (py-ay)(bx-ax)
 * 
 * This is the z-component of the cross product and equals twice the signed area
 * of triangle (a,b,p). The sign indicates which side of line segment ab the
 * point p lies on:
 * - Positive: p is to the left of ab (counter-clockwise)
 * - Zero: p is exactly on line ab
 * - Negative: p is to the right of ab (clockwise)
 * 
 * For a triangle with vertices v0,v1,v2 in counter-clockwise order, a point p
 * is inside the triangle if and only if:
 * edge(v0,v1,p) ≥ 0 AND edge(v1,v2,p) ≥ 0 AND edge(v2,v0,p) ≥ 0
 * 
 * This works because being on the "left" side of all three edges means being
 * inside the triangle. The three edge function values also give us barycentric
 * coordinates (w0,w1,w2) which can interpolate vertex attributes across the
 * triangle's surface, though we don't use that feature here.
 */
float edge_func(Vec3 a, Vec3 b, Vec3 p) {
    return (p.x - a.x)*(b.y - a.y) - (p.y - a.y)*(b.x - a.x);
}

/*
 * Rasterization converts a triangle from three vertex coordinates into a set
 * of pixels. The algorithm tests each pixel in a bounding rectangle to determine
 * if it lies inside the triangle.
 * 
 * The bounding box is computed by taking the minimum and maximum x and y
 * coordinates among the three vertices, giving us a tight rectangular region
 * that fully contains the triangle. This is an optimization - instead of testing
 * all 480,000 pixels in an 800×600 image, we only test pixels that could
 * possibly be inside.
 * 
 * We clamp the bounding box to screen dimensions [0,width)×[0,height) to avoid
 * accessing memory outside the framebuffer array, which would cause crashes or
 * undefined behavior.
 * 
 * For each pixel (x,y) in the bounding box, we test the point (x+0.5, y+0.5)
 * rather than (x,y) because pixel centers are at half-integer coordinates. A
 * pixel at array index (x,y) covers the square [x,x+1)×[y,y+1), so its center
 * is at (x+0.5, y+0.5). Testing at the center gives more accurate coverage.
 * 
 * We compute three edge function values w0, w1, w2 for the three edges of the
 * triangle. By the mathematical property of the edge function, if all three
 * values are non-negative (≥0), the point lies inside the triangle (including
 * on the edges). We then write the triangle's color to the framebuffer at
 * index [y*width + x], which maps 2D coordinates to the 1D array.
 */
void rasterize(Triangle t, unsigned int* pixels, int width, int height) {
    int min_x = fminf(fminf(t.v[0].x, t.v[1].x), t.v[2].x);
    int max_x = fmaxf(fmaxf(t.v[0].x, t.v[1].x), t.v[2].x);
    int min_y = fminf(fminf(t.v[0].y, t.v[1].y), t.v[2].y);
    int max_y = fmaxf(fmaxf(t.v[0].y, t.v[1].y), t.v[2].y);
    
    if (min_x < 0) min_x = 0;
    if (max_x >= width) max_x = width - 1;
    if (min_y < 0) min_y = 0;
    if (max_y >= height) max_y = height - 1;
    
    for (int y = min_y; y <= max_y; y++) {
        for (int x = min_x; x <= max_x; x++) {
            Vec3 p = {x + 0.5f, y + 0.5f, 0};
            float w0 = edge_func(t.v[1], t.v[2], p);
            float w1 = edge_func(t.v[2], t.v[0], p);
            float w2 = edge_func(t.v[0], t.v[1], p);
            
            if (w0 >= 0 && w1 >= 0 && w2 >= 0) {
                pixels[y * width + x] = t.color;
            }
        }
    }
}

/*
 * Backface culling determines whether a triangle faces toward or away from the
 * camera using the mathematical relationship between the surface normal and
 * view direction.
 * 
 * We compute the surface normal n as the cross product of two edges: n = e1×e2
 * where e1 = v1-v0 and e2 = v2-v0. For counter-clockwise vertex ordering, this
 * normal points outward from the front face of the triangle due to the right-
 * hand rule property of cross products.
 * 
 * The view direction vector v points from the surface toward the camera. In our
 * simplified camera model where the camera looks down the -Z axis, v = (0,0,-1).
 * 
 * The dot product n·v tells us the angle between the normal and view direction:
 * - If n·v > 0: angle < 90°, normal points toward camera (front face)
 * - If n·v = 0: angle = 90°, viewing edge-on
 * - If n·v < 0: angle > 90°, normal points away from camera (back face)
 * 
 * Mathematically, n·v < 0 ⟺ cos(θ) < 0 ⟺ θ > 90°, which means the surface
 * is facing away. We return true (should cull) in this case. This optimization
 * skips rendering roughly half the triangles in a closed mesh, since we can
 * never see the back faces of an opaque object.
 */
int should_cull(Vec3 v0, Vec3 v1, Vec3 v2) {
    Vec3 edge1 = sub(v1, v0);
    Vec3 edge2 = sub(v2, v0);
    Vec3 normal = cross(edge1, edge2);
    Vec3 view_dir = {0, 0, -1};
    return dot(normal, view_dir) < 0;
}

/*
 * Save the framebuffer as a PPM image file. PPM is an uncompressed format
 * that's trivial to write: a simple header followed by raw RGB bytes. Most
 * image viewers can open PPM files directly.
 */
void save_ppm(const char* filename, unsigned int* pixels, int w, int h) {
    FILE* f = fopen(filename, "wb");
    fprintf(f, "P6\n%d %d\n255\n", w, h);
    for (int i = 0; i < w * h; i++) {
        unsigned char rgb[3] = {
            (pixels[i] >> 16) & 0xFF,
            (pixels[i] >> 8) & 0xFF,
            pixels[i] & 0xFF
        };
        fwrite(rgb, 1, 3, f);
    }
    fclose(f);
}

int main(void) {
    const int W = 800, H = 600;
    unsigned int* pixels = calloc(W * H, sizeof(unsigned int));
    
    /* 
     * A cube has 8 vertices at the corners. These define the geometry in
     * object space, centered at the origin with sides of length 2.
     */
    Vec3 cube_verts[8] = {
        {-1,-1,-1}, {1,-1,-1}, {1,1,-1}, {-1,1,-1},
        {-1,-1,1}, {1,-1,1}, {1,1,1}, {-1,1,1}
    };
    
    /*
     * The 12 triangles that make up the 6 faces of the cube. Each face uses
     * two triangles, and each triangle references three vertices by index.
     * The order matters for backface culling: vertices are listed counter-
     * clockwise when viewing the front face.
     */
    int cube_faces[12][3] = {
        {0,1,2}, {0,2,3}, {4,7,6}, {4,6,5}, {0,4,5}, {0,5,1},
        {1,5,6}, {1,6,2}, {2,6,7}, {2,7,3}, {3,7,4}, {3,4,0}
    };
    
    /* Face colors alternate between purple, cyan, and yellow */
    unsigned int colors[6] = {
        0xFF00FF, 0xFF00FF, 0x00FFFF, 0x00FFFF, 0xFFFF00, 0xFFFF00
    };
    
    /* Render 60 frames of animation */
    for (int frame = 0; frame < 60; frame++) {
        memset(pixels, 0, W * H * sizeof(unsigned int));
        
        float angle = frame * 0.05f;
        
        /* Transform and project all vertices once */
        Vec3 transformed[8];
        for (int i = 0; i < 8; i++) {
            Vec3 v = cube_verts[i];
            v = rotate_y(v, angle);
            v = rotate_x(v, angle * 0.7f);
            transformed[i] = project(v, 600);
        }
        
        /* Render each triangle if it's front-facing */
        for (int i = 0; i < 12; i++) {
            int i0 = cube_faces[i][0];
            int i1 = cube_faces[i][1];
            int i2 = cube_faces[i][2];
            
            Vec3 v0 = transformed[i0];
            Vec3 v1 = transformed[i1];
            Vec3 v2 = transformed[i2];
            
            /* Skip triangles facing away from camera */
            if (should_cull(v0, v1, v2)) continue;
            
            Triangle tri = {{v0, v1, v2}, colors[i/2]};
            rasterize(tri, pixels, W, H);
        }
        
        char filename[32];
        sprintf(filename, "frame_%03d.ppm", frame);
        save_ppm(filename, pixels, W, H);
        printf("Rendered frame %d/60\n", frame + 1);
    }
    
    free(pixels);
    printf("Done! Convert to video: ffmpeg -i frame_%%03d.ppm output.mp4\n");
    return 0;
}