Implement triangle interpolation based on Byricentric Coordinates

docs/byricentric.tex

@@ -0,0 +1,108 @@
+\documentclass{article}
+\usepackage{amsmath}
+\usepackage{tikz}
+\begin{document}
+\section{Barycentric Coordinates}
+
+\def\xa{0} \def\ya{0}
+\def\xb{1} \def\yb{2}
+\def\xc{3} \def\yc{-1}
+\def\xp{1} \def\yp{0.25}
+\begin{center}
+\begin{tikzpicture}
+  \draw[gray, thick] (\xa,\ya) -- (\xb,\yb);
+  \draw[gray, thick] (\xb,\yb) -- (\xc,\yc);
+  \draw[gray, thick] (\xc,\yc) -- (\xa,\ya);
+  \draw[gray, thick, dotted] (\xa,\ya) -- (\xp,\yp);
+  \draw[gray, thick, dotted] (\xb,\yb) -- (\xp,\yp);
+  \draw[gray, thick, dotted] (\xc,\yc) -- (\xp,\yp);
+  \filldraw[black] (\xa,\ya) circle (2pt) node[anchor=west]{$v_1$};
+  \filldraw[black] (\xb,\yb) circle (2pt) node[anchor=west]{$v_2$};
+  \filldraw[black] (\xc,\yc) circle (2pt) node[anchor=west]{$v_3$};
+  \filldraw[black] (\xp,\yp) circle (2pt) node[anchor=west]{$p$};
+\end{tikzpicture}
+\end{center}
+
+\begin{align}
+  U &= (u_1, u_2, u_3) \\
+  v_1 &= (x_1, y_1) \\
+  v_2 &= (x_2, y_2) \\
+  v_3 &= (x_3, y_3) \\
+  p   &= (x_p, y_p)
+\end{align}
+
+The Barycentric coordinates can be defined in terms of the following relationships:
+
+\begin{align}
+  \begin{cases}
+  & u_1 + u_2 + u_3 = 1 \\
+  & u_1x_1 + u_2x_2 + u_3x_3 = x_p \\
+  & u_1y_1 + u_2y_2 + u_3y_3 = y_p
+  \end{cases}
+\end{align}
+
+
+Let's reduce the amount of varibles in these equations:
+
+\begin{align}
+  & u_3 = 1 - u_1 - u_2 \\
+  & \begin{cases}
+    u_1(x_1 - x_3) + u_2(x_2 - x_3) &= x_p - x_3 \\
+    u_1(y_1 - y_3) + u_2(y_2 - y_3) &= y_p - y_3 \\
+  \end{cases}
+\end{align}
+
+Now we can turn the system of equations into matrix form:
+
+\begin{align}
+  & T =
+  \begin{bmatrix}
+    x_1 - x_3 & x_2 - x_3 \\
+    y_1 - y_3 & y_2 - y_3 \\
+  \end{bmatrix} \\
+  & U = \begin{bmatrix}
+    u1 \\ u2 \\
+  \end{bmatrix}\\
+  & R = \begin{bmatrix}
+    x_p - x_3 \\
+    y_p - y_3 \\
+  \end{bmatrix} \\
+  & D = 1 \\
+  & T \cdot U = R \cdot D
+\end{align}
+
+So the solution is
+
+\begin{align}
+  U = T^{-1} \cdot R
+\end{align}
+
+So the main effort goes towards finding $T^{-1}$
+
+\begin{align}
+  & T^{-1} = \frac{adj(T)}{det(T)} \\
+  & det(T) = (x_1 - x_3)(y_2 - y_3) - (x_2 - x_3)(y_1 - y_3) \\
+  & adj(T) = \begin{bmatrix}
+    y_2 - y_3 & x_3 - x_2 \\
+    y_3 - y_1 & x_1 - x_3 \\
+  \end{bmatrix} \\
+  & T^{-1} = \frac{D}{det(T)} \cdot \begin{bmatrix}
+    y_2 - y_3 & x_3 - x_2 \\
+    y_3 - y_1 & x_1 - x_3 \\
+  \end{bmatrix} \\
+  & T^{-1}\cdot R = \frac{D}{det(T)} \cdot \begin{bmatrix}
+    (y_2 - y_3)(x_p - x_3) + (x_3 - x_2)(y_p - y_3) \\
+    (y_3 - y_1)(x_p - x_3) + (x_1 - x_3)(y_p - y_3) \\
+  \end{bmatrix}
+\end{align}
+
+So, the final formula you need to find $(u_1, u_2, u_3)$ given points $v_1, v_2, v_3, p$ is
+
+\begin{align}
+  u_1 &= \frac{(y_2 - y_3)(x_p - x_3) + (x_3 - x_2)(y_p - y_3)}{(x_1 - x_3)(y_2 - y_3) - (x_2 - x_3)(y_1 - y_3)} \\
+  u_2 &= \frac{(y_3 - y_1)(x_p - x_3) + (x_1 - x_3)(y_p - y_3)}{(x_1 - x_3)(y_2 - y_3) - (x_2 - x_3)(y_1 - y_3)} \\
+  u_3 &= 1 - u_2 - u_1
+\end{align}
+
+\end{document}
+

examples/triangle.c

@@ -53,7 +53,7 @@ uint32_t *render(float dt)
         rotate_point(&x1, &y1);
         rotate_point(&x2, &y2);
         rotate_point(&x3, &y3);
-        olivec_triangle(oc, x1, y1, x2, y2, x3, y3, 0xFF2020AA);
+        olivec_triangle3(oc, x1, y1, x2, y2, x3, y3, 0xFF2020AA, 0xFF20AA20, 0xFFAA2020);
     }
 
     // Circle

olive.c

@@ -285,6 +285,107 @@ OLIVECDEF void olivec_line(Olivec_Canvas oc, int x1, int y1, int x2, int y2, uin
     }
 }
 
+uint32_t mix_colors3(uint32_t c1, uint32_t c2, uint32_t c3,
+                     int t1, int t2, int t3, int den)
+{
+    int64_t r1 = OLIVEC_RED(c1);
+    int64_t g1 = OLIVEC_GREEN(c1);
+    int64_t b1 = OLIVEC_BLUE(c1);
+    int64_t a1 = OLIVEC_ALPHA(c1);
+
+    int64_t r2 = OLIVEC_RED(c2);
+    int64_t g2 = OLIVEC_GREEN(c2);
+    int64_t b2 = OLIVEC_BLUE(c2);
+    int64_t a2 = OLIVEC_ALPHA(c2);
+
+    int64_t r3 = OLIVEC_RED(c3);
+    int64_t g3 = OLIVEC_GREEN(c3);
+    int64_t b3 = OLIVEC_BLUE(c3);
+    int64_t a3 = OLIVEC_ALPHA(c3);
+
+    int64_t r4 = (r1*t1 + r2*t2 + r3*t3)/den;
+    int64_t g4 = (g1*t1 + g2*t2 + g3*t3)/den;
+    int64_t b4 = (b1*t1 + b2*t2 + b3*t3)/den;
+    int64_t a4 = (a1*t1 + a2*t2 + a3*t3)/den;
+
+    return OLIVEC_RGBA(r4, g4, b4, a4);
+}
+
+void barycentric(int x1, int y1, int x2, int y2, int x3, int y3,
+                 int xp, int yp,
+                 int *u1, int *u2, int *det)
+{
+    *det = ((x1 - x3)*(y2 - y3) - (x2 - x3)*(y1 - y3));
+    *u1  = ((y2 - y3)*(xp - x3) + (x3 - x2)*(yp - y3));
+    *u2  = ((y3 - y1)*(xp - x3) + (x1 - x3)*(yp - y3));
+    // u3 = det - u1 - u2
+}
+
+OLIVECDEF void olivec_triangle3(Olivec_Canvas oc, int x1, int y1, int x2, int y2, int x3, int y3,
+                                uint32_t c1, uint32_t c2, uint32_t c3)
+{
+    if (y1 > y2) {
+        OLIVEC_SWAP(int, x1, x2);
+        OLIVEC_SWAP(int, y1, y2);
+        OLIVEC_SWAP(int, c1, c2);
+    }
+
+    if (y2 > y3) {
+        OLIVEC_SWAP(int, x2, x3);
+        OLIVEC_SWAP(int, y2, y3);
+        OLIVEC_SWAP(int, c2, c3);
+    }
+
+    if (y1 > y2) {
+        OLIVEC_SWAP(int, x1, x2);
+        OLIVEC_SWAP(int, y1, y2);
+        OLIVEC_SWAP(int, c1, c2);
+    }
+
+    int dx12 = x2 - x1;
+    int dy12 = y2 - y1;
+    int dx13 = x3 - x1;
+    int dy13 = y3 - y1;
+
+    for (int y = y1; y <= y2; ++y) {
+        // TODO: move boundary checks outside of loops in olivec_fill_triangle
+        if (0 <= y && (size_t) y < oc.height) {
+            int s1 = dy12 != 0 ? (y - y1)*dx12/dy12 + x1 : x1;
+            int s2 = dy13 != 0 ? (y - y1)*dx13/dy13 + x1 : x1;
+            if (s1 > s2) OLIVEC_SWAP(int, s1, s2);
+            for (int x = s1; x <= s2; ++x) {
+                if (0 <= x && (size_t) x < oc.width) {
+                    int u1, u2, det;
+                    barycentric(x1, y1, x2, y2, x3, y3, x, y, &u1, &u2, &det);
+                    uint32_t color = mix_colors3(c1, c2, c3, u1, u2, det - u1 - u2, det);
+                    olivec_blend_color(&OLIVEC_PIXEL(oc, x, y), color);
+                }
+            }
+        }
+    }
+
+    int dx32 = x2 - x3;
+    int dy32 = y2 - y3;
+    int dx31 = x1 - x3;
+    int dy31 = y1 - y3;
+
+    for (int y = y2; y <= y3; ++y) {
+        if (0 <= y && (size_t) y < oc.height) {
+            int s1 = dy32 != 0 ? (y - y3)*dx32/dy32 + x3 : x3;
+            int s2 = dy31 != 0 ? (y - y3)*dx31/dy31 + x3 : x3;
+            if (s1 > s2) OLIVEC_SWAP(int, s1, s2);
+            for (int x = s1; x <= s2; ++x) {
+                if (0 <= x && (size_t) x < oc.width) {
+                    int u1, u2, det;
+                    barycentric(x1, y1, x2, y2, x3, y3, x, y, &u1, &u2, &det);
+                    uint32_t color = mix_colors3(c1, c2, c3, u1, u2, det - u1 - u2, det);
+                    olivec_blend_color(&OLIVEC_PIXEL(oc, x, y), color);
+                }
+            }
+        }
+    }
+}
+
 // TODO: AA for triangle
 OLIVECDEF void olivec_triangle(Olivec_Canvas oc, int x1, int y1, int x2, int y2, int x3, int y3, uint32_t color)
 {