raster

main.c (14230B)

---

 #include <math.h>
 #include <stdio.h>
 #include <string.h>
 #include <rcx/all.h>
 
 typedef f32 Vec2[2];
 typedef f32 Vec3[3];
 typedef f32 Vec4[4];
 typedef f32 Mat3[3][3];
 typedef f32 Aff3[2][3];
 typedef Vec2 Edge2[2];
 typedef Vec2 Quad2[2];
 typedef struct cycle2 Cycle2;
 typedef struct image2 Image2;
 typedef struct linear_gradient LinearGradient;
 typedef void Shader(Vec4 *pixel, Vec2 pos, f32 alpha, void *arg);
 
 struct cycle2 {
 	usize len;
 	Vec2 *verts;
 };
 
 // TODO: More format options (not only f32 RGBA)
 struct image2 {
 	usize w, h;
 	Vec4 *pixels; // Row major order, starting in upper left
 };
 
 struct linear_gradient {
 	Vec2 pos0;
 	Vec4 color0;
 	Vec2 pos1;
 	Vec4 color1;
 };
 
 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_sub(Vec2 *r, Vec2 u, Vec2 v);
 void vec3_sub(Vec3 *r, Vec3 u, Vec3 v);
 void vec4_sub(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);
 f32 vec2_dot(Vec2 u, Vec2 v);
 f32 vec3_dot(Vec3 u, Vec3 v);
 f32 vec4_dot(Vec4 u, Vec4 v);
 f32 vec2_hypot(Vec2 v);
 f32 vec2_dist(Vec2 u, Vec2 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, Shader shader, void *shader_arg);
 
 f32
 clampf(f32 x, f32 min, f32 max) {
 	x = x < min ? min : x;
 	x = x > max ? max : x;
 	return x;
 }
 
 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]; }
 
 void vec2_sub(Vec2 *r, Vec2 u, Vec2 v) { for (usize i = 0; i < 2; i++) (*r)[i] = u[i] - v[i]; }
 void vec3_sub(Vec3 *r, Vec3 u, Vec3 v) { for (usize i = 0; i < 3; i++) (*r)[i] = u[i] - v[i]; }
 void vec4_sub(Vec4 *r, Vec4 u, Vec4 v) { for (usize i = 0; i < 4; i++) (*r)[i] = u[i] - v[i]; }
 
 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]; }
 
 f32 vec2_dot(Vec2 u, Vec2 v) { return u[0] * v[0] + u[1] * v[1]; }
 f32 vec3_dot(Vec3 u, Vec3 v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2]; }
 f32 vec4_dot(Vec4 u, Vec4 v) { return u[0] * v[0] + u[1] * v[1] + u[2] * v[2] + u[3] * v[3]; }
 
 f32 vec2_hypot(Vec2 v) { return hypotf(v[0], v[1]); }
 f32 vec2_dist(Vec2 u, Vec2 v) { return hypotf(v[0] - u[0], v[1] - u[1]); }
 
 void
 aff3_mul(Aff3 *r, Aff3 a, Aff3 b) {
 	SET(*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],
 		}
 	}));
 }
 
 void
 aff3_app(Vec2 *r, Aff3 a, Vec2 b) {
 	SET(*r, ((Vec2){
 		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],
 	}));
 }
 
 void
 aff3_scale(Aff3 *r, f32 x, f32 y) {
 	SET(*r, ((Aff3){
 		{    x, 0.0f, 0.0f },
 		{ 0.0f,    y, 0.0f },
 	}));
 }
 
 void
 aff3_translate(Aff3 *r, f32 x, f32 y) {
 	SET(*r, ((Aff3){
 		{ 1.0f, 0.0f, x },
 		{ 0.0f, 1.0f, y },
 	}));
 }
 
 // Find the affine transformation that maps src to dst.
 void
 aff3_view(Aff3 *r, Quad2 dst, Quad2 src) {
 	f32 dstw = quad2_width(dst);
 	f32 dsth = quad2_height(dst);
 	f32 srcw = quad2_width(src);
 	f32 srch = quad2_height(src);
 
 	ASSERT(dstw > 0 && dsth > 0 && srcw > 0 && srch > 0);
 
 	Aff3 t0, s, t1;
 	aff3_translate(&t0, -src[0][0], -src[0][0]);
 	aff3_scale(&s, dstw / srcw, dsth / srch);
 	aff3_translate(&t1, dst[0][0], dst[0][1]);
 	aff3_mul(r, s, t0);
 	aff3_mul(r, t1, *r);
 }
 
 f32 quad2_width(Quad2 q) { return q[1][0] - q[0][0]; }
 f32 quad2_height(Quad2 q) { return q[1][1] - q[0][1]; }
 
 void
 quad2_transform(Quad2 *r, Aff3 t, Quad2 q) {
 	Vec2 v0, v1;
 	aff3_app(&v0, t, q[0]);
 	aff3_app(&v1, t, q[1]);
 	SET(*r, ((Quad2){
 		{ MIN(v0[0], v1[0]), MIN(v0[1], v1[1]) },
 		{ MAX(v0[0], v1[0]), MAX(v0[1], v1[1]) },
 	}));
 }
 
 void
 image2_bbox(Quad2 *r, Image2 image) {
 	SET(*r, ((Quad2){
 		{ 0.0f, 0.0f },
 		{ (f32)image.w, (f32)image.h },
 	}));
 }
 
 // TODO: Support filling collections of cycles (like Postscript), not only one
 int
 fill(Image2 image, Quad2 src_view, Cycle2 cycle, Shader shader, void *shader_arg) {
 	if (cycle.len < 3) return -1;
 
 	Quad2 dst_view;
 	image2_bbox(&dst_view, image);
 	Aff3 view_matrix;
 	aff3_view(&view_matrix, dst_view, src_view);
 
 	// TODO: Copy cycle, apply view_matrix to it (to avoid recalculuations),
 	// and sort by top vertex (to avoid repeatedly considering irrelevant
 	// edges). Benchmark the difference.
 
 	f32 *mem = r_alloc((2 * image.w + 1) * sizeof *mem);
 	if (!mem) return -1;
 	f32 *values = mem;
 	f32 *deltas = mem + image.w;
 
 	for (usize i = 0; i < image.h; i++) {
 		Quad2 row_bbox = {
 			{ 0.0f, (f32)(image.h - i - 1) },
 			{ (f32)image.w, (f32)(image.h - i) },
 		};
 
 		memset(values, 0, image.w * sizeof *values);
 		memset(deltas, 0, (image.w + 1) * sizeof *deltas);
 
 		for (usize k = 0; k < cycle.len; k++) {
 			Edge2 e;
 			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).
 			Edge2 c;
 
 			// 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]) {
 				eb = e[0], et = e[1];
 				cb = c[0], ct = c[1];
 			} else {
 				eb = e[1], et = e[0];
 				cb = c[1], ct = c[0];
 			}
 
 			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
 
 			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;
 
 			// 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 {
 				cb[0] = eb[0];
 				cb[1] = eb[1];
 			}
 			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 {
 				ct[0] = et[0];
 				ct[1] = et[1];
 			}
 
 			// 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];
 
 			// TODO: Diagrams for each case. ┌ ┐ ┘ └
 
 			// Case 1: The clipped edge is entirely to the right of the image.
 			if (row_bbox[1][0] <= cl[0]) {
 				// No pixels are right of the edge, so we do nothing.
 				continue;
 			}
 
 			// Case 2: The clipped edge is entirely to the left of the image.
 			//             ╔═════════⋯═════════╗
 			//             ⋮                   ⋮          
 			//         *   ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║╮
 			//        /    ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║├ Filled on previous iters
 			//       /     ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║╯
 			//      /      ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║← Filled on current iter
 			//     /       ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║╮
 			//    /        ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║├ Filled on following iters
 			//   *         ║▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒▒║╯
 			//             ⋮                   ⋮          
 			//             ╚═════════⋯═════════╝
 			// TODO: Misleading picture: the edge is clipped to the row.
 			if (cr[0] <= row_bbox[0][0]) {
 				deltas[0] += c[0][1] - c[1][1]; // Note the sign.
 				continue;
 			}
 
 			// 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.
 			if (coll == colr) {
 				// Compute the area of the trapezoid covering the pixel that
 				// intersects the edge.
 				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;
 
 				usize j = coll;
 				ASSERT(j < image.w);
 				values[j] += area;
 				deltas[j + 1] += h;
 
 				continue;
 			}
 
 			// Case 4: 
 			{
 				Vec2 p;
 				p[0] = MAX(coll + 1.0f, 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[1] - cl[1];
 
 				// 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;
 					}
 				}
 
 				// Account for the pentagon in the right-most pixel.
 				if (cr[0] < row_bbox[1][0]) {
 					Vec2 q;
 					q[0] = colr;
 					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 y = ((f32)(image.h - i) - 0.5f) / (f32)image.h;
 		f32 s = 0.0f;
 		for (usize j = 0; j < image.w; j++) {
 			s += deltas[j];
 			f32 alpha = clampf(values[j] + s, 0.0f, 1.0f);
 			Vec4 *pixel = &image.pixels[i * image.h + j];
 			f32 x = ((f32)j + 0.5f) / (f32)image.w;
 			shader(pixel, (Vec2){x, y}, alpha, shader_arg);
 		}
 	}
 
 	free(mem);
 
 	return 0;
 }
 
 void
 shader_over(Vec4 *pixel, Vec2 pos, f32 alpha, void *arg) {
 	Vec4 *color = arg;
 
 	Vec4 c;
 	vec4_mul(&c, alpha, *color);
 	vec4_mul(pixel, 1.0f - c[3], *pixel);
 	vec4_add(pixel, c, *pixel);
 }
 
 void
 shader_linear_gradient(Vec4 *pixel, Vec2 pos, f32 alpha, void *arg) {
 	LinearGradient *lg = arg;
 
 	// Project pos onto the gradient's axis.
 	// TODO: This is wrong somehow.
 	Vec2 v; vec2_sub(&v, lg->pos1, lg->pos0);
 	Vec2 w; vec2_sub(&w, pos, lg->pos0);
 	f32 t = vec2_dot(v, w) / vec2_hypot(v);
 
 	Vec4 color0; vec4_mul(&color0, 1.0f - t, lg->color0);
 	Vec4 color1; vec4_mul(&color1, t, lg->color1);
 	Vec4 color; vec4_add(&color, color0, color1);
 
 	vec4_mul(&color, alpha, color);
 	vec4_mul(pixel, 1.0f - color[3], *pixel);
 	vec4_add(pixel, color, *pixel);
 }
 
 void
 write_farbfeld(FILE *f, Image2 image) {
 	char header[16];
 	memcpy(header, "farbfeld", 8);
 	r_writeb32(header + 8, image.w);
 	r_writeb32(header + 12, image.h);
 	if (!fwrite(header, sizeof header, 1, f))
 		r_fatalf("write_farbfeld: failed to write header");
 
 	for (usize i = 0; i < image.h; i++) {
 		for (usize j = 0; j < image.w; j++) {
 			Vec4 *p = &image.pixels[i * image.w + j];
 
 			// TODO: Clamp pixel values?
 			u16 buf[4];
 			for (usize k = 0; k < 3; k++) {
 				f32 a = (*p)[k] / (*p)[3]; // Un-premultiply alpha
 				buf[k] = r_htob16(a * (f32)U16_MAX);
 			}
 			buf[3] = r_htob16((*p)[3] * (f32)U16_MAX);
 
 			if (!fwrite(buf, sizeof buf, 1, f))
 				r_fatalf("write_farbfeld: failed to write pixel data");
 		}
 	}
 }
 
 void
 fill_under(Image2 image, f32 f(f32), usize n, Shader shader, void *shader_arg) {
 	Cycle2 cycle = {
 		.len = n + 3,
 		.verts = r_ealloc((n + 3) * sizeof(Vec2)),
 	};
 
 	SET(cycle.verts[0], ((Vec2){0.0f, 0.0f}));
 	SET(cycle.verts[1], ((Vec2){1.0f, 0.0f}));
 	for (usize i = 0; i <= n; i++) {
 		f32 x = 1.0f - (f32)i / (f32)n;
 		SET(cycle.verts[i+2], ((Vec2){x, f(x)}));
 	}
 
 	Quad2 cycle_view = {
 		{ 0.0f, 0.0f },
 		{ 1.0f, 1.0f },
 	};
 
 	if (fill(image, cycle_view, cycle, shader, shader_arg) < 0)
 		r_fatalf("fill error");
 
 	free(cycle.verts);
 }
 
 static f32 pi = acosf(-1.0f);
 f32 waves(f32 x) { return sinf(50.0f * pi * x) * 0.01f + 0.4f; }
 f32 mountains(f32 x) { return sinf(3.0f * pi * x) * 0.2f + sinf(100.0f * pi * x) * 0.003f + 0.4f; }
 
 int
 main(void) {
 	Image2 image = {
 		.w = 1000,
 		.h = 1000,
 		.pixels = r_ealloc(1000 * 1000 * sizeof(Vec4)),
 	};
 
 	Vec4 bg = { 0.6f, 0.6f, 0.65f, 1.0f };
 	for (usize i = 0; i < image.h; i++) {
 		for (usize j = 0; j < image.w; j++)
 			SET(image.pixels[i * image.h + j], bg);
 	}
 
 	fill_under(image, waves, 1000, shader_linear_gradient,
 		&(LinearGradient){
 			.pos0 = { 0.0f, 0.0f },
 			.color0 = { 0.0f, 0.0f, 0.0f, 1.0f },
 			.pos1 = { 0.0f, 0.5f },
 			.color1 = { 0.3f, 0.3f, 0.9f, 1.0f },
 		});
 
 	fill_under(image, mountains, 1000, shader_linear_gradient,
 		&(LinearGradient){
 			.pos0 = { 0.0f, 0.0f },
 			.color0 = { 0.05f, 0.15f, 0.05f, 1.0f },
 			.pos1 = { 0.0f, 1.0f },
 			.color1 = { 0.15f, 0.5f, 0.15f, 1.0f },
 		});
 
 	write_farbfeld(stdout, image);
 }