commit 43f696b08f98717039bfc0f0b9426540b9c0d1e1
Author: Robert Russell <robertrussell.72001@gmail.com>
Date: Tue, 9 Jul 2024 17:42:57 -0700
Initial commit
Save work so far
Diffstat:
| A | main.c | | | 345 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1 file changed, 345 insertions(+), 0 deletions(-)
diff --git a/main.c b/main.c
@@ -0,0 +1,345 @@
+#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;
+
+struct cycle2 {
+ usize len;
+ Vec2 *verts;
+};
+
+struct image2 {
+ usize w, h;
+ 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];
+}
+
+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; }
+
+Vec4 vec4_add(Vec4 u, Vec4 v) { Vec4 sum; vec_add(4, sum.data, u.data, v.data); return sum; }
+
+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],
+ };
+}
+
+void
+aff3_app(Vec2 *r, Aff3 a, Vec2 b) {
+ *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) {
+ *r = (Aff3){
+ x, 0.0f, 0.0f,
+ 0.0f, y, 0.0f,
+ };
+}
+
+void
+aff3_translate(Aff3 *r, f32 x, f32 y) {
+ *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);
+}
+
+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]; }
+
+void
+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]) },
+ };
+}
+
+void
+image2_bbox(Quad2 *bbox, Image2 image) {
+ *bbox = (Quad2){
+ { 0.0f, 0.0f },
+ { (f32)image.w, (f32)image.h },
+ };
+}
+
+int
+fill(Image2 image, Quad2 src_view, Cycle2 cycle, Vec4 color) {
+ ASSERT(cycle.len >= 3);
+
+ 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 *deltas = r_alloc(image.w * sizeof *deltas);
+ if (!deltas) return -1;
+
+ for (usize row = 0; row < image.h; row++) {
+ Quad2 row_bbox = {
+ { 0.0f, (f32)(image.h - row - 1) },
+ { (f32)image.w, (f32)(image.h - row) },
+ };
+ //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);
+
+ for (usize i = 0; i < cycle.len; i++) {
+ 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];
+
+ // 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;
+
+ f32 *e_bot, *e_top;
+ f32 *c_bot, *c_top;
+ if (e[0][1] < e[1][1]) {
+ e_bot = e[0], e_top = e[1];
+ c_bot = e[0], c_top = c[1];
+ } else {
+ e_bot = e[1], e_top = e[0];
+ c_bot = e[1], c_top = 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
+
+ // 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];
+ }
+
+ // TODO: This is a goofy way of calculuating these
+ f32 b, l;
+ if (eb[1] < rowb) {
+ b = rowb;
+ edge2_solve_x(&l, e, rowb);
+ } else {
+ b = eb[1];
+ l = eb[0];
+ }
+ f32 t, r;
+ if (et[1] > rowt) {
+ t = rowt;
+ edge2_solve_x(&r, e, rowt);
+ } else {
+ t = et[1];
+ r = et[0];
+ }
+
+ f32 dy = t - b;
+ f32 dx = r - l;
+
+ // ┌ ┐ ┘ └
+
+ // Case 1: The clipped edge is entirely to the right of the image.
+ // ╔═════════⋯═════════╗
+ // ║ ║
+ // ⋮ ⋮
+ // ║ ║
+ // * ║ ║
+ // / ║ ║
+ // / ║ ║
+ // / ║ ║
+ // / ║ ║
+ // * ║ ║
+ // ║ ║
+ // ⋮ ⋮
+ // ║ ║
+ // ╚═════════⋯═════════╝
+ if (rowr <= l) {
+ // 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 (r < rowl) {
+ deltas[0] += fabs(t - b);
+ 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);
+
+ // 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
+ // 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;
+
+ // TODO: Explain this.
+ if (col < image.w - 1)
+ deltas[col + 1] += h - a;
+
+ 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;
+ }
+
+ // ...
+
+ if (r < rowr) {
+ }
+ }
+ }
+ }
+
+ free(deltas);
+
+ return 0;
+}
+
+int
+main(void) {
+
+}