Implement endpoint-based elliptical arc drawing

Allow drawing elliptical arcs by specifying end points instead of
center point, following some common vector graphic endpoint
parameterization format.

Ref:
https://www.w3.org/TR/SVG/implnote.html

apps/multi.c

@@ -220,6 +220,55 @@ static void apps_jelly_start(twin_screen_t *screen, int x, int y, int w, int h)
     twin_window_show(window);
 }
 
+static void draw_flower(twin_path_t *path, twin_fixed_t radius)
+{
+    const twin_angle_t angle_72 = 819;
+    const twin_angle_t angle_36 = 409;
+    twin_fixed_t p_x = twin_fixed_mul(radius, twin_cos(-angle_36));
+    twin_fixed_t p_y = twin_fixed_mul(radius, twin_sin(-angle_36));
+    twin_path_move(path, p_x, p_y);
+
+    for (twin_angle_t a = angle_36; a <= TWIN_ANGLE_360; a += angle_72) {
+        twin_fixed_t c_x = twin_fixed_mul(radius, twin_cos(a));
+        twin_fixed_t c_y = twin_fixed_mul(radius, twin_sin(a));
+        twin_fixed_t rx = radius;
+        twin_fixed_t ry = radius * 3;
+        twin_path_arc_ellipse(path, 1, 1, rx, ry, p_x, p_y, c_x, c_y,
+                              a - angle_36);
+        p_x = c_x;
+        p_y = c_y;
+    }
+
+    twin_path_close(path);
+}
+
+static void apps_flower_start(twin_screen_t *screen, int x, int y, int w, int h)
+{
+    twin_window_t *window = twin_window_create(
+        screen, TWIN_ARGB32, TwinWindowApplication, x, y, w, h);
+    twin_pixmap_t *pixmap = window->pixmap;
+    twin_path_t *stroke = twin_path_create();
+    twin_path_t *path = twin_path_create();
+    twin_path_translate(path, D(200), D(200));
+    twin_path_scale(path, D(10), D(10));
+    twin_path_translate(stroke, D(200), D(200));
+    twin_fill(pixmap, 0xffffffff, TWIN_SOURCE, 0, 0, w, h);
+    twin_window_set_name(window, "Flower");
+    twin_path_move(stroke, D(-200), D(0));
+    twin_path_draw(stroke, D(200), D(0));
+    twin_path_move(stroke, D(0), D(200));
+    twin_path_draw(stroke, D(0), D(-200));
+    twin_path_set_cap_style(stroke, TwinCapProjecting);
+    twin_paint_stroke(pixmap, 0xffcc9999, stroke, D(10));
+    draw_flower(path, D(3));
+    // twin_path_arc_circle(path, 0, 1, D(3), D(3), D(0), D(0), D(3));
+    // twin_path_draw(path, 0, 0);
+    // twin_path_close(path);
+    twin_paint_path(pixmap, 0xffe2d2d2, path);
+    twin_path_destroy(stroke);
+    twin_window_show(window);
+}
+
 void apps_multi_start(twin_screen_t *screen,
                       const char *name,
                       int x,
@@ -233,4 +282,5 @@ void apps_multi_start(twin_screen_t *screen,
     apps_quickbrown_start(screen, x += 20, y += 20, w, h);
     apps_ascii_start(screen, x += 20, y += 20, w, h);
     apps_jelly_start(screen, x += 20, y += 20, w / 2, h);
+    apps_flower_start(screen, x += 20, y += 20, w, h);
 }

include/twin.h

@@ -800,6 +800,7 @@ void twin_path_ellipse(twin_path_t *path,
                        twin_fixed_t y,
                        twin_fixed_t x_radius,
                        twin_fixed_t y_radius);
+
 void twin_path_arc(twin_path_t *path,
                    twin_fixed_t x,
                    twin_fixed_t y,
@@ -808,6 +809,26 @@ void twin_path_arc(twin_path_t *path,
                    twin_angle_t start,
                    twin_angle_t extent);
 
+void twin_path_arc_ellipse(twin_path_t *path,
+                           bool large_arc,
+                           bool sweep,
+                           twin_fixed_t radius_x,
+                           twin_fixed_t radius_y,
+                           twin_fixed_t cur_x,
+                           twin_fixed_t cur_y,
+                           twin_fixed_t target_x,
+                           twin_fixed_t target_y,
+                           twin_angle_t rotation);
+
+void twin_path_arc_circle(twin_path_t *path,
+                          bool large_arc,
+                          bool sweep,
+                          twin_fixed_t radius,
+                          twin_fixed_t cur_x,
+                          twin_fixed_t cur_y,
+                          twin_fixed_t target_x,
+                          twin_fixed_t target_y);
+
 void twin_path_rectangle(twin_path_t *path,
                          twin_fixed_t x,
                          twin_fixed_t y,
@@ -1104,6 +1125,10 @@ twin_fixed_t twin_tan(twin_angle_t a);
 
 void twin_sincos(twin_angle_t a, twin_fixed_t *sin, twin_fixed_t *cos);
 
+twin_angle_t twin_atan2(twin_fixed_t y, twin_fixed_t x);
+
+twin_angle_t twin_acos(twin_fixed_t x);
+
 /*
  * widget.c
  */

src/path.c

@@ -236,6 +236,146 @@ void twin_path_arc(twin_path_t *path,
     twin_path_set_matrix(path, save);
 }
 
+static twin_angle_t vector_angle(twin_fixed_t ux,
+                                 twin_fixed_t uy,
+                                 twin_fixed_t vx,
+                                 twin_fixed_t vy)
+{
+    twin_fixed_t dot = twin_fixed_mul(ux, vx) + twin_fixed_mul(uy, vy);
+
+    twin_fixed_t ua =
+        twin_fixed_sqrt(twin_fixed_mul(ux, ux) + twin_fixed_mul(uy, uy));
+    twin_fixed_t va =
+        twin_fixed_sqrt(twin_fixed_mul(vx, vx) + twin_fixed_mul(vy, vy));
+
+    /* cos(theta) = dot / (|u| * |v|) */
+    twin_fixed_t cos_theta = twin_fixed_div(dot, twin_fixed_mul(ua, va));
+    twin_fixed_t cross = twin_fixed_mul(ux, vy) - twin_fixed_mul(uy, vx);
+    twin_angle_t angle = twin_acos(cos_theta);
+    return (cross < 0) ? -angle : angle;
+}
+
+typedef struct {
+    twin_fixed_t cx, cy;
+    twin_angle_t start, extent;
+} twin_ellipse_para_t;
+
+static twin_ellipse_para_t get_center_parameters(twin_fixed_t x1,
+                                                 twin_fixed_t y1,
+                                                 twin_fixed_t x2,
+                                                 twin_fixed_t y2,
+                                                 bool fa,
+                                                 bool fs,
+                                                 twin_fixed_t rx,
+                                                 twin_fixed_t ry,
+                                                 twin_fixed_t phi)
+{
+    twin_fixed_t sin_phi = twin_sin(phi);
+    twin_fixed_t cos_phi = twin_cos(phi);
+
+    /* Simplify through translation/rotation */
+    twin_fixed_t x =
+        twin_fixed_mul(cos_phi, twin_fixed_mul(x1 - x2, TWIN_FIXED_HALF)) +
+        twin_fixed_mul(sin_phi, twin_fixed_mul(y1 - y2, TWIN_FIXED_HALF));
+
+    twin_fixed_t y =
+        twin_fixed_mul(-sin_phi, twin_fixed_mul(x1 - x2, TWIN_FIXED_HALF)) +
+        twin_fixed_mul(cos_phi, twin_fixed_mul(y1 - y2, TWIN_FIXED_HALF));
+
+    twin_fixed_t px = twin_fixed_mul(x, x);
+    twin_fixed_t py = twin_fixed_mul(y, y);
+    twin_fixed_t prx = twin_fixed_mul(rx, rx);
+    twin_fixed_t pry = twin_fixed_mul(ry, ry);
+
+    /* Correct out-of-range radii */
+    twin_fixed_t L = twin_fixed_div(px, prx) + twin_fixed_div(py, pry);
+
+    if (L > TWIN_FIXED_ONE) {
+        twin_fixed_t sqrt_L = twin_fixed_sqrt(L);
+        rx = twin_fixed_mul(sqrt_L, twin_fixed_abs(rx));
+        ry = twin_fixed_mul(sqrt_L, twin_fixed_abs(ry));
+    } else {
+        rx = twin_fixed_abs(rx);
+        ry = twin_fixed_abs(ry);
+    }
+
+    /* Compute center */
+    twin_fixed_t sign = (fa != fs) ? -1 : 1;
+    twin_fixed_t numerator = twin_fixed_mul(prx, pry) -
+                             twin_fixed_mul(prx, py) - twin_fixed_mul(pry, px);
+
+    twin_fixed_t denominator =
+        twin_fixed_mul(prx, py) + twin_fixed_mul(pry, px);
+    twin_fixed_t M =
+        sign * twin_fixed_sqrt(twin_fixed_div(numerator, denominator));
+    twin_fixed_t _cx =
+        twin_fixed_mul(M, twin_fixed_div(twin_fixed_mul(rx, y), ry));
+    twin_fixed_t _cy =
+        twin_fixed_mul(M, twin_fixed_div(twin_fixed_mul(-ry, x), rx));
+
+    twin_ellipse_para_t ret;
+    ret.cx = twin_fixed_mul(cos_phi, _cx) - twin_fixed_mul(sin_phi, _cy) +
+             twin_fixed_mul(x1 + x2, TWIN_FIXED_HALF);
+
+    ret.cy = twin_fixed_mul(sin_phi, _cx) + twin_fixed_mul(cos_phi, _cy) +
+             twin_fixed_mul(y1 + y2, TWIN_FIXED_HALF);
+
+    /* Compute θ and dθ */
+    ret.start = vector_angle(TWIN_FIXED_ONE, 0, twin_fixed_div(x - _cx, rx),
+                             twin_fixed_div(y - _cy, ry));
+    twin_angle_t extent = vector_angle(
+        twin_fixed_div(x - _cx, rx), twin_fixed_div(y - _cy, ry),
+        twin_fixed_div(-x - _cx, rx), twin_fixed_div(-y - _cy, ry));
+
+    if (fs && extent > TWIN_ANGLE_0)
+        extent -= TWIN_ANGLE_360;
+    if (!fs && extent < TWIN_ANGLE_0)
+        extent += TWIN_ANGLE_360;
+    ret.start %= TWIN_ANGLE_360;
+    extent %= TWIN_ANGLE_360;
+
+    ret.extent = extent;
+    return ret;
+}
+
+void twin_path_arc_ellipse(twin_path_t *path,
+                           bool large_arc,
+                           bool sweep,
+                           twin_fixed_t radius_x,
+                           twin_fixed_t radius_y,
+                           twin_fixed_t cur_x,
+                           twin_fixed_t cur_y,
+                           twin_fixed_t target_x,
+                           twin_fixed_t target_y,
+                           twin_angle_t rotation)
+{
+    twin_ellipse_para_t para;
+    para = get_center_parameters(cur_x, cur_y, target_x, target_y, large_arc,
+                                 sweep, radius_x, radius_y, rotation);
+    twin_matrix_t save = twin_path_current_matrix(path);
+
+    twin_path_translate(path, para.cx, para.cy);
+    twin_path_rotate(path, rotation);
+    twin_path_translate(path, -para.cx, -para.cy);
+    twin_path_arc(path, para.cx, para.cy, radius_x, radius_y, para.start,
+                  para.extent);
+
+    twin_path_set_matrix(path, save);
+}
+
+void twin_path_arc_circle(twin_path_t *path,
+                          bool large_arc,
+                          bool sweep,
+                          twin_fixed_t radius,
+                          twin_fixed_t cur_x,
+                          twin_fixed_t cur_y,
+                          twin_fixed_t target_x,
+                          twin_fixed_t target_y)
+{
+    twin_path_arc_ellipse(path, large_arc, sweep, radius, radius, cur_x, cur_y,
+                          target_x, target_y, TWIN_ANGLE_0);
+}
+
 void twin_path_rectangle(twin_path_t *path,
                          twin_fixed_t x,
                          twin_fixed_t y,

src/trig.c

@@ -97,3 +97,102 @@ void twin_sincos(twin_angle_t a, twin_fixed_t *sin, twin_fixed_t *cos)
         *cos = cos_val;
     }
 }
+
+static const twin_angle_t atan_table[] = {
+    0x0200, /* arctan(2^0)  = 45° -> 512     */
+    0x0130, /* arctan(2^-1) = 26.565° -> 303 */
+    0x009B, /* arctan(2^-2) = 14.036° -> 155 */
+    0x004F, /* arctan(2^-3) = 7.125° -> 79   */
+    0x0027, /* arctan(2^-4) = 3.576° -> 39   */
+    0x0014, /* arctan(2^-5) = 1.790° -> 20   */
+    0x000A, /* arctan(2^-6) = 0.895° -> 10   */
+    0x0005, /* arctan(2^-7) = 0.448° -> 5    */
+    0x0002, /* arctan(2^-8) = 0.224° -> 2    */
+    0x0001, /* arctan(2^-9) = 0.112° -> 1    */
+    0x0001, /* arctan(2^-10) = 0.056° -> 1   */
+    0x0000, /* arctan(2^-11) = 0.028° -> 0   */
+};
+
+static twin_angle_t twin_atan2_first_quadrant(twin_fixed_t y, twin_fixed_t x)
+{
+    if (x == 0 && y == 0)
+        return 0;
+    if (x == 0)
+        return TWIN_ANGLE_90;
+    if (y == 0)
+        return TWIN_ANGLE_0;
+    twin_fixed_t current_x = x;
+    twin_fixed_t current_y = y;
+    twin_angle_t angle = 0;
+    /* CORDIC iteration */
+    for (int i = 0; i < 12; i++) {
+        twin_fixed_t temp_x, temp_y;
+        if (current_y > 0) {
+            temp_x = current_x + (current_y >> i);
+            temp_y = current_y - (current_x >> i);
+            angle += atan_table[i];
+        } else {
+            temp_x = current_x - (current_y >> i);
+            temp_y = current_y + (current_x >> i);
+            angle -= atan_table[i];
+        }
+        current_x = temp_x;
+        current_y = temp_y;
+    }
+    return angle;
+}
+
+twin_angle_t twin_atan2(twin_fixed_t y, twin_fixed_t x)
+{
+    if (x == 0 && y == 0)
+        return TWIN_ANGLE_0;
+    if (x == 0)
+        return (y > 0) ? TWIN_ANGLE_90 : TWIN_ANGLE_270;
+    if (y == 0)
+        return (x > 0) ? TWIN_ANGLE_0 : TWIN_ANGLE_180;
+    int quadrant;
+    twin_fixed_t abs_x = x;
+    twin_fixed_t abs_y = y;
+    if (x >= 0 && y >= 0) {
+        quadrant = 1;
+    } else if (x < 0 && y >= 0) {
+        quadrant = 2;
+        abs_x = -x;
+    } else if (x < 0 && y < 0) {
+        quadrant = 3;
+        abs_x = -x;
+        abs_y = -y;
+    } else {
+        quadrant = 4;
+        abs_y = -y;
+    }
+    twin_angle_t angle = twin_atan2_first_quadrant(abs_y, abs_x);
+    switch (quadrant) {
+    case 1:
+        return angle;
+    case 2:
+        return 2048 - angle;
+    case 3:
+        return 2048 + angle;
+    case 4:
+        return 4096 - angle;
+    default:
+        return 0;
+    }
+}
+
+twin_angle_t twin_acos(twin_fixed_t x)
+{
+    if (x <= -TWIN_FIXED_ONE)
+        return TWIN_ANGLE_180;
+    if (x >= TWIN_FIXED_ONE)
+        return TWIN_ANGLE_0;
+    twin_fixed_t y = twin_fixed_sqrt(TWIN_FIXED_ONE - twin_fixed_mul(x, x));
+    twin_angle_t angle;
+    if (x >= 0) {
+        angle = twin_atan2_first_quadrant(y, x);
+    } else {
+        angle = TWIN_ANGLE_180 - twin_atan2_first_quadrant(y, -x);
+    }
+    return angle;
+}