raster

2D vector graphics rasterizer
git clone git://git.rr3.xyz/raster
Log | Files | Refs

commit 77b8f4b38b90b00ece693bb82dcbbfbfd07ef223
parent 43f696b08f98717039bfc0f0b9426540b9c0d1e1
Author: Robert Russell <robertrussell.72001@gmail.com>
Date:   Tue,  9 Jul 2024 22:56:57 -0700

Successfully raster a triangle

Diffstat:
AMakefile | 2++
Mmain.c | 415++++++++++++++++++++++++++++++++++++++++++++++---------------------------------
2 files changed, 243 insertions(+), 174 deletions(-)

diff --git a/Makefile b/Makefile @@ -0,0 +1,2 @@ +main: main.c + gcc -Wall -o main main.c -lm -lrcx diff --git a/main.c b/main.c @@ -1,3 +1,4 @@ +#include <math.h> #include <stdio.h> #include <string.h> #include <rcx/all.h> @@ -22,52 +23,80 @@ struct image2 { Vec4 *pixels; // Row major order, starting in upper left }; -void -vec_add(usize n, f32 *sum, f32 *u, f32 *v) { - for (usize i = 0; i < n; i++) - sum[i] = u[i] + v[i]; +f32 clampf(f32 x, f32 min, f32 max); +void vec2_add(Vec2 *r, Vec2 u, Vec2 v); +void vec3_add(Vec3 *r, Vec3 u, Vec3 v); +void vec4_add(Vec4 *r, Vec4 u, Vec4 v); +void vec2_mul(Vec2 *r, f32 a, Vec2 v); +void vec3_mul(Vec3 *r, f32 a, Vec3 v); +void vec4_mul(Vec4 *r, f32 a, Vec4 v); +void aff3_mul(Aff3 *r, Aff3 a, Aff3 b); +void aff3_app(Vec2 *r, Aff3 a, Vec2 b); +void aff3_scale(Aff3 *r, f32 x, f32 y); +void aff3_translate(Aff3 *r, f32 x, f32 y); +void aff3_view(Aff3 *r, Quad2 dst, Quad2 src); +f32 quad2_width(Quad2 q); +f32 quad2_height(Quad2 q); +void quad2_transform(Quad2 *r, Aff3 t, Quad2 q); +void image2_bbox(Quad2 *r, Image2 image); +int fill(Image2 image, Quad2 src_view, Cycle2 cycle, Vec4 color); + +f32 +clampf(f32 x, f32 min, f32 max) { + x = x < min ? min : x; + x = x > max ? max : x; + return x; } -Vec2 vec2_add(Vec2 u, Vec2 v) { Vec2 sum; vec_add(2, sum.data, u.data, v.data); return sum; } - -Vec3 vec3_add(Vec3 u, Vec3 v) { Vec3 sum; vec_add(3, sum.data, u.data, v.data); return sum; } +void vec2_add(Vec2 *r, Vec2 u, Vec2 v) { for (usize i = 0; i < 2; i++) (*r)[i] = u[i] + v[i]; } +void vec3_add(Vec3 *r, Vec3 u, Vec3 v) { for (usize i = 0; i < 3; i++) (*r)[i] = u[i] + v[i]; } +void vec4_add(Vec4 *r, Vec4 u, Vec4 v) { for (usize i = 0; i < 4; i++) (*r)[i] = u[i] + v[i]; } -Vec4 vec4_add(Vec4 u, Vec4 v) { Vec4 sum; vec_add(4, sum.data, u.data, v.data); return sum; } +void vec2_mul(Vec2 *r, f32 a, Vec2 v) { for (usize i = 0; i < 2; i++) (*r)[i] = a * v[i]; } +void vec3_mul(Vec3 *r, f32 a, Vec3 v) { for (usize i = 0; i < 3; i++) (*r)[i] = a * v[i]; } +void vec4_mul(Vec4 *r, f32 a, Vec4 v) { for (usize i = 0; i < 4; i++) (*r)[i] = a * v[i]; } void aff3_mul(Aff3 *r, Aff3 a, Aff3 b) { - *r = (Aff3){ - a[0][0] * b[0][0] + a[0][1] * b[1][0], - a[0][0] * b[0][1] + a[0][1] * b[1][1], - a[0][0] * b[0][2] + a[0][1] * b[1][2] + a[0][2], - a[1][0] * b[0][0] + a[1][1] * b[1][0], - a[1][0] * b[0][1] + a[1][1] * b[1][1], - a[1][0] * b[0][2] + a[1][1] * b[1][2] + a[1][2], + Aff3 rr = { + { + a[0][0] * b[0][0] + a[0][1] * b[1][0], + a[0][0] * b[0][1] + a[0][1] * b[1][1], + a[0][0] * b[0][2] + a[0][1] * b[1][2] + a[0][2], + }, { + a[1][0] * b[0][0] + a[1][1] * b[1][0], + a[1][0] * b[0][1] + a[1][1] * b[1][1], + a[1][0] * b[0][2] + a[1][1] * b[1][2] + a[1][2], + } }; + memcpy(r, rr, sizeof *r); } void aff3_app(Vec2 *r, Aff3 a, Vec2 b) { - *r = (Vec2){ + Vec2 rr = { a[0][0] * b[0] + a[0][1] * b[1] + a[0][2], a[1][0] * b[0] + a[1][1] * b[1] + a[1][2], }; + memcpy(r, rr, sizeof *r); } void aff3_scale(Aff3 *r, f32 x, f32 y) { - *r = (Aff3){ - x, 0.0f, 0.0f, - 0.0f, y, 0.0f, + Aff3 rr = { + { x, 0.0f, 0.0f }, + { 0.0f, y, 0.0f }, }; + memcpy(r, rr, sizeof *r); } void aff3_translate(Aff3 *r, f32 x, f32 y) { - *r = (Aff3){ - 1.0f, 0.0f, x, - 0.0f, 1.0f, y, + Aff3 rr = { + { 1.0f, 0.0f, x }, + { 0.0f, 1.0f, y }, }; + memcpy(r, rr, sizeof *r); } // Find the affine transformation that maps src to dst. @@ -88,50 +117,6 @@ aff3_view(Aff3 *r, Quad2 dst, Quad2 src) { aff3_mul(r, t1, *r); } -void -edge2_bbox(Quad2 *bbox, Edge2 e) { - *bbox = (Quad2){ - { MIN(e[0][0], e[1][0]), MIN(e[0][1], e[1][1]) }, - { MAX(e[0][0], e[1][0]), MAX(e[0][1], e[1][1]) }, - }; -} - -int -edge2_solve_x(f32 *x, Edge2 e, f32 y) { - // We parameterize the edge as - // (top - bot) * t + bot, t in [0..1] - // where bot[1] <= top[1]. - Vec2 bot, top; - if (e[0][1] < e[1][1]) { - SET(bot, e[0][1]); - SET(top, e[1][1]); - } else { - SET(bot, e[1][1]); - SET(top, e[0][1]); - } - - // If the edge doesn't pass through y, then there is no solution. - if (y < bot[1] || y > top[1]) return -1; - - // If bot[1] == top[1] == y, then there isn't a unique solution, unless the - // edge is a point (x, y), in which case the only solution is x. - if (bot[1] == top[1]) { - if (bot[0] == top[0]) { - *x = bot[0]; - return 0; - } - return -1; - } - - // We now have bot[1] <= y <= top[1], and one of the inequalities is - // strict, so the denominator of the following is nonzero. - f32 t = (y - bot[1]) / (top[1] - bot[1]); - - *x = (top[0] - bot[0]) * t + bot[0]; - - return 0; -} - f32 quad2_width(Quad2 q) { return q[1][0] - q[0][0]; } f32 quad2_height(Quad2 q) { return q[1][1] - q[0][1]; } @@ -140,23 +125,25 @@ quad2_transform(Quad2 *r, Aff3 t, Quad2 q) { Vec2 v0, v1; aff3_app(&v0, t, q[0]); aff3_app(&v1, t, q[1]); - *r = (Quad2){ - (Vec2){ MIN(v0[0], v1[0]), MIN(v0[1], v1[1]) }, - (Vec2){ MAX(v0[0], v1[0]), MAX(v0[1], v1[1]) }, + Quad2 rr = { + { MIN(v0[0], v1[0]), MIN(v0[1], v1[1]) }, + { MAX(v0[0], v1[0]), MAX(v0[1], v1[1]) }, }; + memcpy(r, rr, sizeof *r); } void -image2_bbox(Quad2 *bbox, Image2 image) { - *bbox = (Quad2){ +image2_bbox(Quad2 *r, Image2 image) { + Quad2 rr = { { 0.0f, 0.0f }, { (f32)image.w, (f32)image.h }, }; + memcpy(r, rr, sizeof *r); } int fill(Image2 image, Quad2 src_view, Cycle2 cycle, Vec4 color) { - ASSERT(cycle.len >= 3); + if (cycle.len < 3) return -1; Quad2 dst_view; image2_bbox(&dst_view, image); @@ -167,103 +154,81 @@ fill(Image2 image, Quad2 src_view, Cycle2 cycle, Vec4 color) { // and sort by top vertex (to avoid repeatedly considering irrelevant // edges). Benchmark the difference. - f32 *deltas = r_alloc(image.w * sizeof *deltas); - if (!deltas) return -1; + f32 *mem = r_alloc((2 * image.w + 1) * sizeof *mem); + if (!mem) return -1; + f32 *values = mem; + f32 *deltas = mem + image.w; - for (usize row = 0; row < image.h; row++) { + for (usize i = 0; i < image.h; i++) { Quad2 row_bbox = { - { 0.0f, (f32)(image.h - row - 1) }, - { (f32)image.w, (f32)(image.h - row) }, + { 0.0f, (f32)(image.h - i - 1) }, + { (f32)image.w, (f32)(image.h - i) }, }; - //f32 rowt = (f32)(image.h - row); - //f32 rowb = rowt - 1.0f; - //f32 rowl = 0.0f; - //f32 rowr = (f32)image.w; - memset(deltas, 0, sizeof *deltas); + memset(values, 0, image.w * sizeof *deltas); + memset(deltas, 0, (image.w + 1) * sizeof *deltas); - for (usize i = 0; i < cycle.len; i++) { + for (usize k = 0; k < cycle.len; k++) { Edge2 e; - aff3_app(&e[0], view_matrix, cycle.verts[i]); - aff3_app(&e[1], view_matrix, cycle.verts[(i+1)%cycle.len]); - - f32 dx = e[1][0] - e[0][0]; - f32 dy = e[1][1] - e[0][1]; + aff3_app(&e[0], view_matrix, cycle.verts[k]); + aff3_app(&e[1], view_matrix, cycle.verts[(k+1)%cycle.len]); // c, initialized below, is e vertically clipped to the row. We // can't horizontally clip e to the row, because, e.g., e could be // entirely to the left of the image while still contributing area - // (see Case 1 below). + // (see case 1 below). Edge2 c; - f32 *e_bot, *e_top; - f32 *c_bot, *c_top; + // Compute pointers to the bottom-most (b) and top-most (t) + // vertices of e and c. + f32 *eb, *et; + f32 *cb, *ct; if (e[0][1] < e[1][1]) { - e_bot = e[0], e_top = e[1]; - c_bot = e[0], c_top = c[1]; + eb = e[0], et = e[1]; + cb = c[0], ct = c[1]; } else { - e_bot = e[1], e_top = e[0]; - c_bot = e[1], c_top = c[0]; + eb = e[1], et = e[0]; + cb = c[1], ct = c[0]; } - if (e_top[1] <= row_bbox[0][1]) continue; // e is below the row - if (row_bbox[1][1] <= e_bot[1]) continue; // e is above the row + if (et[1] <= row_bbox[0][1]) continue; // e is below the row + if (row_bbox[1][1] <= eb[1]) continue; // e is above the row - // Initialize c. - if (e_bot[1] < row_bbox[0][1]) { - c_bot[0] = dx / dy * (row_bbox[0][1] - e_bot[1]) + e_bot[0]; - c_bot[1] = row_bbox[0][1]; - } else { - c_bot[0] = e_bot[0]; - c_bot[1] = e_bot[1]; - } - if (e_top[1] > row_bbox[1][1]) { - c_top[0] = dx / dy * (row_bbox[1][1] - e_top[1]) + e_top[0]; - c_top[1] = row_bbox[1][1]; - } else { - c_top[0] = e_top[0]; - c_top[1] = e_top[1]; - } + f32 e_dx = e[1][0] - e[0][0]; + f32 e_dy = e[1][1] - e[0][1]; + f32 dx_dy = e_dx / e_dy; + f32 dy_dx = e_dy / e_dx; - // TODO: This is a goofy way of calculuating these - f32 b, l; - if (eb[1] < rowb) { - b = rowb; - edge2_solve_x(&l, e, rowb); + // Initialize c (through cb and ct). + if (eb[1] < row_bbox[0][1]) { + // Intersect e with the bottom of the row. + cb[0] = dx_dy * (row_bbox[0][1] - eb[1]) + eb[0]; + cb[1] = row_bbox[0][1]; } else { - b = eb[1]; - l = eb[0]; + cb[0] = eb[0]; + cb[1] = eb[1]; } - f32 t, r; - if (et[1] > rowt) { - t = rowt; - edge2_solve_x(&r, e, rowt); + if (et[1] > row_bbox[1][1]) { + // Intersect e with the top of the row. + ct[0] = dx_dy * (row_bbox[1][1] - et[1]) + et[0]; + ct[1] = row_bbox[1][1]; } else { - t = et[1]; - r = et[0]; + ct[0] = et[0]; + ct[1] = et[1]; } - f32 dy = t - b; - f32 dx = r - l; + // Compute pointers to the left-most (l) and right-most (r) + // vertices of c. + f32 *cl, *cr; + if (e[0][0] < e[1][0]) + cl = c[0], cr = c[1]; + else + cl = c[1], cr = c[0]; // ┌ ┐ ┘ └ // Case 1: The clipped edge is entirely to the right of the image. - // ╔═════════⋯═════════╗ - // ║ ║ - // ⋮ ⋮ - // ║ ║ - // * ║ ║ - // / ║ ║ - // / ║ ║ - // / ║ ║ - // / ║ ║ - // * ║ ║ - // ║ ║ - // ⋮ ⋮ - // ║ ║ - // ╚═════════⋯═════════╝ - if (rowr <= l) { + if (row_bbox[1][0] <= cl[0]) { // No pixels are right of the edge, so we do nothing. continue; } @@ -281,65 +246,167 @@ fill(Image2 image, Quad2 src_view, Cycle2 cycle, Vec4 color) { // ⋮ ⋮ // ╚═════════⋯═════════╝ // TODO: Misleading picture: the edge is clipped to the row. - if (r < rowl) { - deltas[0] += fabs(t - b); + if (cr[0] <= row_bbox[0][0]) { + deltas[0] += c[0][1] - c[1][1]; // Note the sign. continue; } - // The indices of the first and last pixels intersecting the edge. - // These could be out of bounds. - isize col0 = floorf(l); - isize col1 = floorf(r); + // The indices of the left-most and right-most pixels intersecting + // the edge. These could be out of bounds. + f32 coll = floorf(cl[0]); + f32 colr = floorf(cr[0]); // Case 3: The clipped edge intersects exactly one pixel in the // current row. // TODO: Explain why we must separate this case. - // (Non infinite slope.) - if (col0 == col1) { - isize col = col0; - - // The failure of cases 1 and 2 implies col is in bounds. - ASSERT(0 <= col && col < image.w); - - // Compute the area a of the trapezoid covering the pixel that + if (coll == colr) { + // Compute the area of the trapezoid covering the pixel that // intersects the edge. - f32 w0 = (f32)(col + 1) - l; - f32 w1 = (f32)(col + 1) - r; - f32 h = dy; - f32 a = (w0 + w1) / 2.0f * h; - deltas[col] += a; + f32 w0 = (coll + 1.0f) - c[0][0]; + f32 w1 = (coll + 1.0f) - c[1][0]; + f32 h = c[0][1] - c[1][1]; // Note the sign. + f32 area = (w0 + w1) / 2.0f * h; - // TODO: Explain this. - if (col < image.w - 1) - deltas[col + 1] += h - a; + usize j = coll; + ASSERT(j < image.w); + values[j] += area; + deltas[j + 1] += h; continue; } - f32 x0 = l > rowl ? ceilf(l) : rowl; - // Case 4: { - f32 dy_dx = dy / dx; // dx != 0, but it could be small... - - if (rowl <= l) { - f32 iy = dy_dx * ( - l) + b; + Vec2 p; + p[0] = MAX(cl[0], row_bbox[0][0]); + p[1] = dy_dx * (p[0] - cl[0]) + cl[1]; + + // Assume e goes from SW to NE. Then we can multiply areas + // by this sign for correctness in the other three possible + // configurations. + // TODO: Explain with diagram. + // TODO: Use this convention in other cases too? + f32 sign = cl == c[0] ? -1.0f : +1.0f; + + // Account for the triangle in the left-most pixel. + if (row_bbox[0][0] <= cl[0]) { + f32 w = p[0] - cl[0]; + f32 h = p[1] - cl[1]; + f32 area = w * h / 2.0f; + + usize j = coll; + ASSERT(j < image.w); + values[j] += sign * area; } - // ... + f32 rect_area = p[0] - cl[0]; + + // Account for the trapezoids in the middle pixels. + f32 trap_coll = MAX(coll + 1.0f, 0.0f); + f32 trap_colr = MIN(colr - 1.0f, (f32)(image.w - 1)); + if (trap_coll <= trap_colr) { + usize jl = trap_coll; + usize jr = trap_colr; + for (usize j = jl; j <= jr; j++) { + values[j] += sign * (rect_area + dy_dx / 2.0f); + rect_area += dy_dx; + } + } - if (r < rowr) { + // Account for the pentagon in the right-most pixel. + if (cr[0] < row_bbox[1][0]) { + Vec2 q; + q[0] = cr[0]; + q[1] = dy_dx * (q[0] - cr[0]) + cr[1]; + + // Compute area of trapezoid on top of the rectangle. + // TODO: Diagram. + f32 w0 = (colr + 1.0f) - cr[0]; + f32 w1 = 1.0f; + f32 h = cr[1] - q[1]; + f32 trap_area = (w0 + w1) / 2.0f * h; + + usize j = colr; + ASSERT(j < image.w); + values[j] += sign * (rect_area + trap_area); + deltas[j + 1] += sign * (cr[1] - cl[1]); } } } + + f32 s = 0.0f; + for (usize j = 0; j < image.w; j++) { + s += deltas[j]; + f32 a = clampf(values[j] + s, 0.0f, 1.0f); + Vec4 *p = &image.pixels[i * image.h + j]; + // TODO: More blending options. + vec4_mul(p, 1.0f - a, *p); + vec4_add(p, color, *p); + } } - free(deltas); + free(mem); return 0; } +void +image2_write_farbfeld(FILE *f, Image2 image) { + usize w = image.w; + usize h = image.h; + + char header[16]; + memcpy(header, "farbfeld", 8); + r_writeb32(header + 8, w); + r_writeb32(header + 12, h); + if (!fwrite(header, sizeof header, 1, f)) + r_fatalf("image2_write_farbfeld: failed to write header"); + + for (usize i = 0; i < h; i++) { + for (usize j = 0; j < w; j++) { + Vec4 *p = &image.pixels[i * w + j]; + + u16 buf[4]; + for (usize k = 0; k < 3; k++) + buf[k] = r_htob16(((*p)[k] / (*p)[3]) * (f32)U16_MAX); + buf[3] = r_htob16((*p)[3] * (f32)U16_MAX); + + if (!fwrite(buf, sizeof buf, 1, f)) + r_fatalf("image2_write_farbfeld: failed to write pixel data"); + } + } +} + int main(void) { - + Vec4 white = { 0.0f, 0.0f, 0.0f, 1.0f }; + Vec4 black = { 1.0f, 1.0f, 1.0f, 1.0f }; + + Cycle2 cycle = { + .len = 3, + .verts = (Vec2[]){ + { 0.0f, 0.0f }, + { 1.0f, 0.0f }, + { 0.5f, 1.0f }, + }, + }; + Quad2 cycle_view = { + { 0.0f, 0.0f }, + { 1.0f, 1.0f }, + }; + + Image2 image = { + .w = 1000, + .h = 1000, + .pixels = r_ealloc(1000 * 1000 * sizeof(Vec4)), + }; + for (usize i = 0; i < image.h; i++) { + for (usize j = 0; j < image.w; j++) + memcpy(&image.pixels[i * image.h + j], white, sizeof white); + } + + if (fill(image, cycle_view, cycle, black) < 0) + r_fatalf("fill error"); + + image2_write_farbfeld(stdout, image); }