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,54 @@ 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,
+                        int number_of_petals)
+{
+    const twin_angle_t angle_shift = TWIN_ANGLE_360 / number_of_petals;
+    const twin_angle_t angle_start = angle_shift / 2;
+    twin_fixed_t p_x = twin_fixed_mul(radius, twin_cos(-angle_start));
+    twin_fixed_t p_y = twin_fixed_mul(radius, twin_sin(-angle_start));
+    twin_path_move(path, p_x, p_y);
+
+    for (twin_angle_t a = angle_start; a <= TWIN_ANGLE_360; a += angle_shift) {
+        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_start);
+        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), 5);
+    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 +281,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

@@ -670,7 +670,7 @@ void twin_clear_file(twin_file_t *file);
  */
 
 #define twin_fixed_mul(a, b) ((twin_fixed_t) (((int64_t) (a) * (b)) >> 16))
-#define twin_fixed_div(a, b) ((twin_fixed_t) ((((int64_t) (a)) << 16) / b))
+#define twin_fixed_div(a, b) ((twin_fixed_t) ((((int64_t) (a)) << 16) / (b)))
 
 twin_fixed_t twin_fixed_sqrt(twin_fixed_t a);
 
@@ -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,152 @@ 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) = (u ⋅ v) / (|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_param_t;
+
+static twin_ellipse_param_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;
+    double px_d = twin_fixed_to_double(px);
+    double py_d = twin_fixed_to_double(py);
+    double prx_d = twin_fixed_to_double(prx);
+    double pry_d = twin_fixed_to_double(pry);
+
+    twin_fixed_t A =
+        twin_double_to_fixed(pry_d / (py_d + (pry_d / prx_d) * px_d));
+    twin_fixed_t B =
+        twin_double_to_fixed(prx_d / (prx_d + (px_d / py_d) * pry_d));
+    twin_fixed_t C =
+        twin_double_to_fixed(pry_d / (pry_d + (py_d / px_d) * prx_d));
+
+    twin_fixed_t pM = A - B - C;
+    twin_fixed_t M = sign * twin_fixed_sqrt(pM);
+    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_param_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_param_t param;
+    param = 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, param.cx, param.cy);
+    twin_path_rotate(path, rotation);
+    twin_path_translate(path, -param.cx, -param.cy);
+    twin_path_arc(path, param.cx, param.cy, radius_x, radius_y, param.start,
+                  param.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,83 @@ 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 TWIN_ANGLE_0;
+    if (x == 0)
+        return TWIN_ANGLE_90;
+    if (y == 0)
+        return TWIN_ANGLE_0;
+    twin_angle_t angle = 0;
+    /* CORDIC iteration */
+    for (int i = 0; i < 12; i++) {
+        twin_fixed_t temp_x = x;
+        if (y > 0) {
+            x += (y >> i);
+            y -= (temp_x >> i);
+            angle += atan_table[i];
+        } else {
+            x -= (y >> i);
+            y += (temp_x >> i);
+            angle -= atan_table[i];
+        }
+    }
+    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;
+    twin_fixed_t x_sign_mask = x >> 31;
+    twin_fixed_t abs_x = (x ^ x_sign_mask) - x_sign_mask;
+    twin_fixed_t y_sign_mask = y >> 31;
+    twin_fixed_t abs_y = (y ^ y_sign_mask) - y_sign_mask;
+    twin_fixed_t m = ((~x_sign_mask & ~y_sign_mask) * 0) +
+                     ((x_sign_mask & ~y_sign_mask) * 1) +
+                     ((x_sign_mask & y_sign_mask) * 1) +
+                     ((~x_sign_mask & y_sign_mask) * 2);
+    twin_fixed_t sign = 1 - 2 * (x_sign_mask ^ y_sign_mask);
+    twin_angle_t angle = twin_atan2_first_quadrant(abs_y, abs_x);
+    /* First quadrant  : angle
+     * Second quadrant : 180 - angle
+     * Third quadrant  : 180 + angle
+     * Fourth quadrant : 360 - angle
+     */
+    return TWIN_ANGLE_180 * m + sign * angle;
+}
+
+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));
+    if (x >= 0) {
+        return twin_atan2_first_quadrant(y, x);
+    } else {
+        return TWIN_ANGLE_180 - twin_atan2_first_quadrant(y, -x);
+    }
+}