Merge pull request #4523 from xtsm/ellipse

Implemented another ellipse drawing algorithm

Tests/test_imagechops.py

@@ -311,8 +311,8 @@ def test_subtract():
 
     # Assert
     assert new.getbbox() == (25, 50, 76, 76)
-    assert new.getpixel((50, 50)) == GREEN
-    assert new.getpixel((50, 51)) == BLACK
+    assert new.getpixel((50, 51)) == GREEN
+    assert new.getpixel((50, 52)) == BLACK
 
 
 def test_subtract_scale_offset():
@@ -350,8 +350,8 @@ def test_subtract_modulo():
 
     # Assert
     assert new.getbbox() == (25, 50, 76, 76)
-    assert new.getpixel((50, 50)) == GREEN
-    assert new.getpixel((50, 51)) == BLACK
+    assert new.getpixel((50, 51)) == GREEN
+    assert new.getpixel((50, 52)) == BLACK
 
 
 def test_subtract_modulo_no_clip():

Tests/test_imagedraw.py

@@ -164,6 +164,19 @@ def test_arc_width_non_whole_angle():
     assert_image_similar_tofile(im, expected, 1)
 
 
+def test_arc_high():
+    # Arrange
+    im = Image.new("RGB", (200, 200))
+    draw = ImageDraw.Draw(im)
+
+    # Act
+    draw.arc([10, 10, 89, 189], 20, 330, width=20, fill="white")
+    draw.arc([110, 10, 189, 189], 20, 150, width=20, fill="white")
+
+    # Assert
+    assert_image_equal(im, Image.open("Tests/images/imagedraw_arc_high.png"))
+
+
 def test_bitmap():
     # Arrange
     im = Image.new("RGB", (W, H))
@@ -194,13 +207,11 @@ def helper_chord(mode, bbox, start, end):
 def test_chord1():
     for mode in ["RGB", "L"]:
         helper_chord(mode, BBOX1, 0, 180)
-        helper_chord(mode, BBOX1, 0.5, 180.4)
 
 
 def test_chord2():
     for mode in ["RGB", "L"]:
         helper_chord(mode, BBOX2, 0, 180)
-        helper_chord(mode, BBOX2, 0.5, 180.4)
 
 
 def test_chord_width():
@@ -240,6 +251,18 @@ def test_chord_zero_width():
         assert_image_equal(im, expected)
 
 
+def test_chord_too_fat():
+    # Arrange
+    im = Image.new("RGB", (100, 100))
+    draw = ImageDraw.Draw(im)
+
+    # Act
+    draw.chord([-150, -150, 99, 99], 15, 60, width=10, fill="white", outline="red")
+
+    # Assert
+    assert_image_equal(im, Image.open("Tests/images/imagedraw_chord_too_fat.png"))
+
+
 def helper_ellipse(mode, bbox):
     # Arrange
     im = Image.new(mode, (W, H))
@@ -282,15 +305,18 @@ def test_ellipse_edge():
     draw = ImageDraw.Draw(im)
 
     # Act
-    draw.ellipse(((0, 0), (W - 1, H)), fill="white")
+    draw.ellipse(((0, 0), (W - 1, H - 1)), fill="white")
 
     # Assert
     assert_image_similar_tofile(im, "Tests/images/imagedraw_ellipse_edge.png", 1)
 
 
 def test_ellipse_symmetric():
-    for bbox in [(25, 25, 76, 76), (25, 25, 75, 75)]:
-        im = Image.new("RGB", (101, 101))
+    for width, bbox in (
+        (100, (24, 24, 75, 75)),
+        (101, (25, 25, 75, 75)),
+    ):
+        im = Image.new("RGB", (width, 100))
         draw = ImageDraw.Draw(im)
         draw.ellipse(bbox, fill="green", outline="blue")
         assert_image_equal(im, im.transpose(Image.FLIP_LEFT_RIGHT))
@@ -345,6 +371,43 @@ def test_ellipse_zero_width():
         assert_image_equal(im, expected)
 
 
+def ellipse_various_sizes_helper(filled):
+    ellipse_sizes = range(32)
+    image_size = sum(ellipse_sizes) + len(ellipse_sizes) + 1
+    im = Image.new("RGB", (image_size, image_size))
+    draw = ImageDraw.Draw(im)
+
+    x = 1
+    for w in ellipse_sizes:
+        y = 1
+        for h in ellipse_sizes:
+            border = [x, y, x + w - 1, y + h - 1]
+            if filled:
+                draw.ellipse(border, fill="white")
+            else:
+                draw.ellipse(border, outline="white")
+            y += h + 1
+        x += w + 1
+
+    return im
+
+
+def test_ellipse_various_sizes():
+    im = ellipse_various_sizes_helper(False)
+
+    with Image.open("Tests/images/imagedraw_ellipse_various_sizes.png") as expected:
+        assert_image_equal(im, expected)
+
+
+def test_ellipse_various_sizes_filled():
+    im = ellipse_various_sizes_helper(True)
+
+    with Image.open(
+        "Tests/images/imagedraw_ellipse_various_sizes_filled.png"
+    ) as expected:
+        assert_image_equal(im, expected)
+
+
 def helper_line(points):
     # Arrange
     im = Image.new("RGB", (W, H))
@@ -420,13 +483,13 @@ def helper_pieslice(bbox, start, end):
 
 
 def test_pieslice1():
-    helper_pieslice(BBOX1, -90, 45)
-    helper_pieslice(BBOX1, -90.5, 45.4)
+    helper_pieslice(BBOX1, -92, 46)
+    helper_pieslice(BBOX1, -92.2, 46.2)
 
 
 def test_pieslice2():
-    helper_pieslice(BBOX2, -90, 45)
-    helper_pieslice(BBOX2, -90.5, 45.4)
+    helper_pieslice(BBOX2, -92, 46)
+    helper_pieslice(BBOX2, -92.2, 46.2)
 
 
 def test_pieslice_width():
@@ -467,6 +530,18 @@ def test_pieslice_zero_width():
         assert_image_equal(im, expected)
 
 
+def test_pieslice_wide():
+    # Arrange
+    im = Image.new("RGB", (200, 100))
+    draw = ImageDraw.Draw(im)
+
+    # Act
+    draw.pieslice([0, 0, 199, 99], 190, 170, width=10, fill="white", outline="red")
+
+    # Assert
+    assert_image_equal(im, Image.open("Tests/images/imagedraw_pieslice_wide.png"))
+
+
 def helper_point(points):
     # Arrange
     im = Image.new("RGB", (W, H))

Tests/test_imagedraw2.py

@@ -79,7 +79,7 @@ def test_ellipse_edge():
     brush = ImageDraw2.Brush("white")
 
     # Act
-    draw.ellipse(((0, 0), (W - 1, H)), brush)
+    draw.ellipse(((0, 0), (W - 1, H - 1)), brush)
 
     # Assert
     assert_image_similar(im, Image.open("Tests/images/imagedraw_ellipse_edge.png"), 1)

src/PIL/ImageDraw.py

@@ -118,7 +118,7 @@ class ImageDraw:
                 fill = self.draw.draw_ink(fill)
         return ink, fill
 
-    def arc(self, xy, start, end, fill=None, width=0):
+    def arc(self, xy, start, end, fill=None, width=1):
         """Draw an arc."""
         ink, fill = self._getink(fill)
         if ink is not None:

src/_imaging.c

@@ -2836,8 +2836,7 @@ _draw_arc(ImagingDrawObject* self, PyObject* args)
     int ink;
     int width = 0;
     float start, end;
-    int op = 0;
-    if (!PyArg_ParseTuple(args, "Offi|ii", &data, &start, &end, &ink, &width)) {
+    if (!PyArg_ParseTuple(args, "Offi|i", &data, &start, &end, &ink, &width)) {
         return NULL;
     }
 
@@ -2854,7 +2853,7 @@ _draw_arc(ImagingDrawObject* self, PyObject* args)
     n = ImagingDrawArc(self->image->image,
                        (int) xy[0], (int) xy[1],
                        (int) xy[2], (int) xy[3],
-                       start, end, &ink, width, op
+                       start, end, &ink, width, self->blend
                        );
 
     free(xy);

src/libImaging/Draw.c

@@ -35,6 +35,7 @@
 #include "Imaging.h"
 
 #include <math.h>
+#include <stdint.h>
 
 #define CEIL(v)  (int) ceil(v)
 #define FLOOR(v) ((v) >= 0.0 ? (int) (v) : (int) floor(v))
@@ -818,221 +819,750 @@ ImagingDrawBitmap(Imaging im, int x0, int y0, Imaging bitmap, const void* ink,
 /* -------------------------------------------------------------------- */
 /* standard shapes */
 
-#define ARC 0
-#define CHORD 1
-#define PIESLICE 2
+// Imagine 2D plane and ellipse with center in (0, 0) and semi-major axes a and b.
+// Then quarter_* stuff approximates its top right quarter (x, y >= 0) with integer
+// points from set {(2x+x0, 2y+y0) | x,y in Z} where x0, y0 are from {0, 1} and
+// are such that point (a, b) is in the set.
 
-static void
-ellipsePoint(int cx, int cy, int w, int h,
-             float i, int *x, int *y)
-{
-    float i_cos, i_sin;
-    float x_f, y_f;
-    double modf_int;
-    i_cos = cos(i*M_PI/180);
-    i_sin = sin(i*M_PI/180);
-    x_f = (i_cos * w/2) + cx;
-    y_f = (i_sin * h/2) + cy;
-    if (modf(x_f, &modf_int) == 0.5) {
-        *x = i_cos > 0 ? FLOOR(x_f) : CEIL(x_f);
+typedef struct {
+    int32_t a, b, cx, cy, ex, ey;
+    int64_t a2, b2, a2b2;
+    int8_t finished;
+} quarter_state;
+
+void quarter_init(quarter_state* s, int32_t a, int32_t b) {
+    if (a < 0 || b < 0) {
+        s->finished = 1;
     } else {
-        *x = FLOOR(x_f + 0.5);
+        s->a = a;
+        s->b = b;
+        s->cx = a;
+        s->cy = b % 2;
+        s->ex = a % 2;
+        s->ey = b;
+        s->a2 = a * a;
+        s->b2 = b * b;
+        s->a2b2 = s->a2 * s->b2;
+        s->finished = 0;
     }
-    if (modf(y_f, &modf_int) == 0.5) {
-        *y = i_sin > 0 ? FLOOR(y_f) : CEIL(y_f);
+}
+
+// deviation of the point from ellipse curve, basically a substitution
+// of the point into the ellipse equation
+int64_t quarter_delta(quarter_state* s, int64_t x, int64_t y) {
+    return llabs(s->a2 * y * y + s->b2 * x * x - s->a2b2);
+}
+
+int8_t quarter_next(quarter_state* s, int32_t* ret_x, int32_t* ret_y) {
+    if (s->finished) {
+        return -1;
+    }
+    *ret_x = s->cx;
+    *ret_y = s->cy;
+    if (s->cx == s->ex && s->cy == s->ey) {
+        s->finished = 1;
     } else {
-        *y = FLOOR(y_f + 0.5);
+        // Bresenham's algorithm, possible optimization: only consider 2 of 3
+        // next points depending on current slope
+        int32_t nx = s->cx;
+        int32_t ny = s->cy + 2;
+        int64_t ndelta = quarter_delta(s, nx, ny);
+        if (nx > 1) {
+            int64_t newdelta = quarter_delta(s, s->cx - 2, s->cy + 2);
+            if (ndelta > newdelta) {
+                nx = s->cx - 2;
+                ny = s->cy + 2;
+                ndelta = newdelta;
+            }
+            newdelta = quarter_delta(s, s->cx - 2, s->cy);
+            if (ndelta > newdelta) {
+                nx = s->cx - 2;
+                ny = s->cy;
             }
         }
-
-static int
-ellipse(Imaging im, int x0, int y0, int x1, int y1,
-        float start, float end, const void* ink_, int fill,
-        int width, int mode, int op)
-{
-    float i;
-    int inner;
-    int n;
-    int maxEdgeCount;
-    int w, h;
-    int x, y;
-    int cx, cy;
-    int lx = 0, ly = 0;
-    int sx = 0, sy = 0;
-    int lx_inner = 0, ly_inner = 0;
-    int sx_inner = 0, sy_inner = 0;
-    DRAW* draw;
-    INT32 ink;
-    Edge* e;
-
-    DRAWINIT();
-
-    while (end < start) {
-        end += 360;
+        s->cx = nx;
+        s->cy = ny;
     }
-
-    if (end - start > 360) {
-        // no need to go in loops
-        end = start + 361;
-    }
-
-    w = x1 - x0;
-    h = y1 - y0;
-    if (w <= 0 || h <= 0) {
     return 0;
 }
 
-    cx = (x0 + x1) / 2;
-    cy = (y0 + y1) / 2;
+// quarter_* stuff can "draw" a quarter of an ellipse with thickness 1, great.
+// Now we use ellipse_* stuff to join all four quarters of two different sized
+// ellipses and receive horizontal segments of a complete ellipse with
+// specified thickness.
+//
+// Still using integer grid with step 2 at this point (like in quarter_*)
+// to ease angle clipping in future.
 
-    if (!fill && width <= 1) {
-        for (i = start; i < end+1; i++) {
-            if (i > end) {
-                i = end;
-            }
-            ellipsePoint(cx, cy, w, h, i, &x, &y);
-            if (i != start) {
-                draw->line(im, lx, ly, x, y, ink);
+typedef struct {
+    quarter_state st_o, st_i;
+    int32_t py, pl, pr;
+    int32_t cy[4], cl[4], cr[4];
+    int8_t bufcnt;
+    int8_t finished;
+    int8_t leftmost;
+} ellipse_state;
+
+void ellipse_init(ellipse_state* s, int32_t a, int32_t b, int32_t w) {
+    s->bufcnt = 0;
+    s->leftmost = a % 2;
+    quarter_init(&s->st_o, a, b);
+    if (w < 1 || quarter_next(&s->st_o, &s->pr, &s->py) == -1) {
+        s->finished = 1;
     } else {
-                sx = x, sy = y;
+        s->finished = 0;
+        quarter_init(&s->st_i, a - 2 * (w - 1), b - 2 * (w - 1));
+        s->pl = s->leftmost;
     }
-            lx = x, ly = y;
 }
 
-        if (i != start) {
-            if (mode == PIESLICE) {
-                if (x != cx || y != cy) {
-                    draw->line(im, x, y, cx, cy, ink);
-                    draw->line(im, cx, cy, sx, sy, ink);
-                }
-            } else if (mode == CHORD) {
-                if (x != sx || y != sy) {
-                    draw->line(im, x, y, sx, sy, ink);
-                }
+int8_t ellipse_next(ellipse_state* s, int32_t* ret_x0, int32_t* ret_y, int32_t* ret_x1) {
+    if (s->bufcnt == 0) {
+        if (s->finished) {
+            return -1;
         }
+        int32_t y = s->py;
+        int32_t l = s->pl;
+        int32_t r = s->pr;
+        int32_t cx = 0, cy = 0;
+        int8_t next_ret;
+        while ((next_ret = quarter_next(&s->st_o, &cx, &cy)) != -1 && cy <= y) {
         }
+        if (next_ret == -1) {
+            s->finished = 1;
         } else {
-        inner = (mode == ARC || !fill) ? 1 : 0;
-
-        // Build edge list
-        // malloc check UNDONE, FLOAT?
-        maxEdgeCount = ceil(end - start);
-        if (inner) {
-            maxEdgeCount *= 2;
+            s->pr = cx;
+            s->py = cy;
         }
-        maxEdgeCount += 3;
-        e = calloc(maxEdgeCount, sizeof(Edge));
-        if (!e) {
+        while ((next_ret = quarter_next(&s->st_i, &cx, &cy)) != -1 && cy <= y) {
+            l = cx;
+        }
+        s->pl = next_ret == -1 ? s->leftmost : cx;
+
+        if ((l > 0 || l < r) && y > 0) {
+            s->cl[s->bufcnt] = l == 0 ? 2 : l;
+            s->cy[s->bufcnt] = y;
+            s->cr[s->bufcnt] = r;
+            ++s->bufcnt;
+        }
+        if (y > 0) {
+            s->cl[s->bufcnt] = -r;
+            s->cy[s->bufcnt] = y;
+            s->cr[s->bufcnt] = -l;
+            ++s->bufcnt;
+        }
+        if (l > 0 || l < r) {
+            s->cl[s->bufcnt] = l == 0 ? 2 : l;
+            s->cy[s->bufcnt] = -y;
+            s->cr[s->bufcnt] = r;
+            ++s->bufcnt;
+        }
+        s->cl[s->bufcnt] = -r;
+        s->cy[s->bufcnt] = -y;
+        s->cr[s->bufcnt] = -l;
+        ++s->bufcnt;
+    }
+    --s->bufcnt;
+    *ret_x0 = s->cl[s->bufcnt];
+    *ret_y = s->cy[s->bufcnt];
+    *ret_x1 = s->cr[s->bufcnt];
+    return 0;
+}
+
+// Clipping tree consists of half-plane clipping nodes and combining nodes.
+// We can throw a horizontal segment in such a tree and collect an ordered set
+// of resulting disjoint clipped segments organized into a sorted linked list
+// of their end points.
+typedef enum {
+  CT_AND, // intersection
+  CT_OR, // union
+  CT_CLIP // half-plane clipping
+} clip_type;
+
+typedef struct clip_node {
+  clip_type type;
+  double a, b, c; // half-plane coeffs, only used in clipping nodes
+  struct clip_node* l; // child pointers, are only non-NULL in combining nodes
+  struct clip_node* r;
+} clip_node;
+
+// Linked list for the ends of the clipped horizontal segments.
+// Since the segment is always horizontal, we don't need to store Y coordinate.
+typedef struct event_list {
+  int32_t x;
+  int8_t type; // used internally, 1 for the left end (smaller X), -1 for the
+               // right end; pointless in output since the output segments
+               // are disjoint, therefore the types would always come in pairs
+               // and interchange (1 -1 1 -1 ...)
+  struct event_list* next;
+} event_list;
+
+// Mirrors all the clipping nodes of the tree relative to the y = x line.
+void clip_tree_transpose(clip_node* root) {
+  if (root != NULL) {
+    if (root->type == CT_CLIP) {
+      double t = root->a;
+      root->a = root->b;
+      root->b = t;
+    }
+    clip_tree_transpose(root->l);
+    clip_tree_transpose(root->r);
+  }
+}
+
+// Outputs a sequence of open-close events (types -1 and 1) for
+// non-intersecting segments sorted by X coordinate.
+// Combining nodes (AND, OR) may also accept sequences for intersecting
+// segments, i.e. something like correct bracket sequences.
+int clip_tree_do_clip(clip_node* root, int32_t x0, int32_t y, int32_t x1, event_list** ret) {
+  if (root == NULL) {
+    event_list* start = malloc(sizeof(event_list));
+    if (!start) {
       ImagingError_MemoryError();
       return -1;
     }
-
-        // Outer circle
-        n = 0;
-        for (i = start; i < end+1; i++) {
-            if (i > end) {
-                i = end;
+    event_list* end = malloc(sizeof(event_list));
+    if (!end) {
+      free(start);
+      ImagingError_MemoryError();
+      return -1;
+    }
+    start->x = x0;
+    start->type = 1;
+    start->next = end;
+    end->x = x1;
+    end->type = -1;
+    end->next = NULL;
+    *ret = start;
+    return 0;
+  }
+  if (root->type == CT_CLIP) {
+    double eps = 1e-9;
+    double A = root->a;
+    double B = root->b;
+    double C = root->c;
+    if (fabs(A) < eps) {
+      if (B * y + C < -eps) {
+        x0 = 1;
+        x1 = 0;
       }
-            ellipsePoint(cx, cy, w, h, i, &x, &y);
-            if (i == start) {
-                sx = x, sy = y;
     } else {
-                add_edge(&e[n++], lx, ly, x, y);
+      // X of intersection
+      double ix = - (B * y + C) / A;
+      if (A * x0 + B * y + C < eps) {
+        x0 = lround(fmax(x0, ix));
       }
-            lx = x, ly = y;
+      if (A * x1 + B * y + C < eps) {
+        x1 = lround(fmin(x1, ix));
       }
-        if (n == 0) {
+    }
+    if (x0 <= x1) {
+      event_list* start = malloc(sizeof(event_list));
+      if (!start) {
+        ImagingError_MemoryError();
+        return -1;
+      }
+      event_list* end = malloc(sizeof(event_list));
+      if (!end) {
+        free(start);
+        ImagingError_MemoryError();
+        return -1;
+      }
+      start->x = x0;
+      start->type = 1;
+      start->next = end;
+      end->x = x1;
+      end->type = -1;
+      end->next = NULL;
+      *ret = start;
+    } else {
+      *ret = NULL;
+    }
+    return 0;
+  }
+  if (root->type == CT_OR || root->type == CT_AND) {
+    event_list* l1;
+    event_list* l2;
+    if (clip_tree_do_clip(root->l, x0, y, x1, &l1) < 0) {
+      return -1;
+    }
+    if (clip_tree_do_clip(root->r, x0, y, x1, &l2) < 0) {
+      while (l1) {
+        l2 = l1->next;
+        free(l1);
+        l1 = l2;
+      }
+      return -1;
+    }
+    *ret = NULL;
+    event_list* tail = NULL;
+    int32_t k1 = 0;
+    int32_t k2 = 0;
+    while (l1 != NULL || l2 != NULL) {
+      event_list* t;
+      if (l2 == NULL || (l1 != NULL && (l1->x < l2->x || (l1->x == l2->x && l1->type > l2->type)))) {
+        t = l1;
+        k1 += t->type;
+        assert(k1 >= 0);
+        l1 = l1->next;
+      } else {
+        t = l2;
+        k2 += t->type;
+        assert(k2 >= 0);
+        l2 = l2->next;
+      }
+      t->next = NULL;
+      if ((root->type == CT_OR && (
+              (t->type == 1 && (tail == NULL || tail->type == -1)) ||
+              (t->type == -1 && k1 == 0 && k2 == 0)
+          )) || 
+          (root->type == CT_AND && (
+              (t->type == 1 && (tail == NULL || tail->type == -1) && k1 > 0 && k2 > 0) ||
+              (t->type == -1 && tail != NULL && tail->type == 1 && (k1 == 0 || k2 == 0))
+          ))) {
+        if (tail == NULL) {
+          *ret = t;
+        } else {
+          tail->next = t;
+        }
+        tail = t;
+      } else {
+        free(t);
+      }
+    }
+    return 0;
+  }
+  *ret = NULL;
   return 0;
 }
 
-        if (inner) {
-            // Inner circle
-            x0 += width - 1;
-            y0 += width - 1;
-            x1 -= width - 1;
-            y1 -= width - 1;
+// One more layer of processing on top of the regular ellipse.
+// Uses the clipping tree.
+// Used for producing ellipse derivatives such as arc, chord, pie, etc.
+typedef struct {
+  ellipse_state st;
+  clip_node* root;
+  clip_node nodes[7];
+  int32_t node_count;
+  event_list* head;
+  int32_t y;
+} clip_ellipse_state;
 
-            w = x1 - x0;
-            h = y1 - y0;
-            if (w <= 0 || h <= 0) {
-                // ARC with no gap in the middle is a PIESLICE
-                mode = PIESLICE;
-                inner = 0;
+typedef void (*clip_ellipse_init)(clip_ellipse_state*, int32_t, int32_t, int32_t, float, float);
+
+void debug_clip_tree(clip_node* root, int space) {
+  if (root == NULL) {
+    return;
+  }
+  if (root->type == CT_CLIP) {
+    int t = space;
+    while (t--) {
+      fputc(' ', stderr);
+    }
+    fprintf(stderr, "clip %+fx%+fy%+f > 0\n", root->a, root->b, root->c);
   } else {
-                for (i = start; i < end+1; i++) {
-                    if (i > end) {
-                        i = end;
+    debug_clip_tree(root->l, space + 2);
+    int t = space;
+    while (t--) {
+      fputc(' ', stderr);
     }
-                    ellipsePoint(cx, cy, w, h, i, &x, &y);
-                    if (i == start) {
-                        sx_inner = x, sy_inner = y;
+    fprintf(stderr, "%s\n", root->type == CT_AND ? "and" : "or");
+    debug_clip_tree(root->r, space + 2);
+  }
+  if (space == 0) {
+    fputc('\n', stderr);
+  }
+}
+
+// Resulting angles will satisfy 0 <= al < 360, al <= ar <= al + 360
+void normalize_angles(float* al, float* ar) {
+  if (*ar - *al >= 360) {
+    *al = 0;
+    *ar = 360;
   } else {
-                        add_edge(&e[n++], lx_inner, ly_inner, x, y);
+    *al = fmod(*al < 0 ? 360 - (fmod(-*al, 360)) : *al, 360);
+    *ar = *al + fmod(*ar < *al ? 360 - fmod(*al - *ar, 360) : *ar - *al, 360);
+  }
+}
+
+// An arc with caps orthogonal to the ellipse curve.
+void arc_init(clip_ellipse_state* s, int32_t a, int32_t b, int32_t w, float al, float ar) {
+  if (a < b) {
+    // transpose the coordinate system
+    arc_init(s, b, a, w, 90 - ar, 90 - al);
+    ellipse_init(&s->st, a, b, w);
+    clip_tree_transpose(s->root);
+  } else {
+    // a >= b, based on "wide" ellipse
+    ellipse_init(&s->st, a, b, w);
+
+    s->head = NULL;
+    s->node_count = 0;
+    normalize_angles(&al, &ar);
+
+    // building clipping tree, a lot of different cases
+    if (ar == al + 360) {
+      s->root = NULL;
+    } else {
+      clip_node* lc = s->nodes + s->node_count++;
+      clip_node* rc = s->nodes + s->node_count++;
+      lc->l = lc->r = rc->l = rc->r = NULL;
+      lc->type = rc->type = CT_CLIP;
+      lc->a = -a * sin(al * M_PI / 180.0);
+      lc->b = b * cos(al * M_PI / 180.0);
+      lc->c = (a * a - b * b) * sin(al * M_PI / 90.0) / 2.0;
+      rc->a = a * sin(ar * M_PI / 180.0);
+      rc->b = -b * cos(ar * M_PI / 180.0);
+      rc->c = (b * b - a * a) * sin(ar * M_PI / 90.0) / 2.0;
+      if (fmod(al, 180) == 0 || fmod(ar, 180) == 0) {
+        s->root = s->nodes + s->node_count++;
+        s->root->l = lc;
+        s->root->r = rc;
+        s->root->type = ar - al < 180 ? CT_AND : CT_OR;
+      } else if (((int)(al / 180) + (int)(ar / 180)) % 2 == 1) {
+        s->root = s->nodes + s->node_count++;
+        s->root->l = s->nodes + s->node_count++;
+        s->root->l->l = s->nodes + s->node_count++;
+        s->root->l->r = lc;
+        s->root->r = s->nodes + s->node_count++;
+        s->root->r->l = s->nodes + s->node_count++;
+        s->root->r->r = rc;
+        s->root->type = CT_OR;
+        s->root->l->type = CT_AND;
+        s->root->r->type = CT_AND;
+        s->root->l->l->type = CT_CLIP;
+        s->root->r->l->type = CT_CLIP;
+        s->root->l->l->l = s->root->l->l->r = NULL;
+        s->root->r->l->l = s->root->r->l->r = NULL;
+        s->root->l->l->a = s->root->l->l->c = 0;
+        s->root->r->l->a = s->root->r->l->c = 0;
+        s->root->l->l->b = (int)(al / 180) % 2 == 0 ? 1 : -1;
+        s->root->r->l->b = (int)(ar / 180) % 2 == 0 ? 1 : -1;
+      } else {
+        s->root = s->nodes + s->node_count++;
+        s->root->l = s->nodes + s->node_count++;
+        s->root->r = s->nodes + s->node_count++;
+        s->root->type = s->root->l->type = ar - al < 180 ? CT_AND : CT_OR;
+        s->root->l->l = lc;
+        s->root->l->r = rc;
+        s->root->r->type = CT_CLIP;
+        s->root->r->l = s->root->r->r = NULL;
+        s->root->r->a = s->root->r->c = 0;
+        s->root->r->b = ar < 180 || ar > 540 ? 1 : -1;
       }
-                    lx_inner = x, ly_inner = y;
     }
   }
 }
 
-        if (end - start < 360) {
-            // Close polygon
-            if (mode == PIESLICE) {
-                if (x != cx || y != cy) {
-                    add_edge(&e[n++], sx, sy, cx, cy);
-                    add_edge(&e[n++], cx, cy, lx, ly);
-                    if (inner) {
-                        ImagingDrawWideLine(im, sx, sy, cx, cy, &ink, width, op);
-                        ImagingDrawWideLine(im, cx, cy, lx, ly, &ink, width, op);
+// A chord line.
+void chord_line_init(clip_ellipse_state* s, int32_t a, int32_t b, int32_t w, float al, float ar) {
+  ellipse_init(&s->st, a, b, a + b + 1);
+
+  s->head = NULL;
+  s->node_count = 0;
+
+  // line equation for chord
+  double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
+  double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
+  s->root = s->nodes + s->node_count++;
+  s->root->l = s->nodes + s->node_count++;
+  s->root->r = s->nodes + s->node_count++;
+  s->root->type = CT_AND;
+  s->root->l->type = s->root->r->type = CT_CLIP;
+  s->root->l->l = s->root->l->r = s->root->r->l = s->root->r->r = NULL;
+  s->root->l->a = yr - yl;
+  s->root->l->b = xl - xr;
+  s->root->l->c = -(s->root->l->a * xl + s->root->l->b * yl);
+  s->root->r->a = -s->root->l->a;
+  s->root->r->b = -s->root->l->b;
+  s->root->r->c = 2 * w * sqrt(pow(s->root->l->a, 2.0) + pow(s->root->l->b, 2.0)) - s->root->l->c;
 }
+
+// Pie side.
+void pie_side_init(clip_ellipse_state* s, int32_t a, int32_t b, int32_t w, float al, float _) {
+  ellipse_init(&s->st, a, b, a + b + 1);
+
+  s->head = NULL;
+  s->node_count = 0;
+
+  double xl = a * cos(al * M_PI / 180.0);
+  double yl = b * sin(al * M_PI / 180.0);
+  double a1 = -yl;
+  double b1 = xl;
+  double c1 = w * sqrt(a1 * a1 + b1 * b1);
+
+  s->root = s->nodes + s->node_count++;
+  s->root->type = CT_AND;
+  s->root->l = s->nodes + s->node_count++;
+  s->root->l->type = CT_AND;
+
+  clip_node* cnode;
+  cnode = s->nodes + s->node_count++;
+  cnode->l = cnode->r = NULL;
+  cnode->type = CT_CLIP;
+  cnode->a = a1;
+  cnode->b = b1;
+  cnode->c = c1;
+  s->root->l->l = cnode;
+  cnode = s->nodes + s->node_count++;
+  cnode->l = cnode->r = NULL;
+  cnode->type = CT_CLIP;
+  cnode->a = -a1;
+  cnode->b = -b1;
+  cnode->c = c1;
+  s->root->l->r = cnode;
+  cnode = s->nodes + s->node_count++;
+  cnode->l = cnode->r = NULL;
+  cnode->type = CT_CLIP;
+  cnode->a = b1;
+  cnode->b = -a1;
+  cnode->c = 0;
+  s->root->r = cnode;
 }
-            } else if (mode == CHORD) {
-                add_edge(&e[n++], sx, sy, lx, ly);
-                if (inner) {
-                    add_edge(&e[n++], sx_inner, sy_inner, lx_inner, ly_inner);
+
+// A chord.
+void chord_init(clip_ellipse_state* s, int32_t a, int32_t b, int32_t w, float al, float ar) {
+  ellipse_init(&s->st, a, b, w);
+
+  s->head = NULL;
+  s->node_count = 0;
+
+  // line equation for chord
+  double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
+  double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
+  s->root = s->nodes + s->node_count++;
+  s->root->l = s->root->r = NULL;
+  s->root->type = CT_CLIP;
+  s->root->a = yr - yl;
+  s->root->b = xl - xr;
+  s->root->c = -(s->root->a * xl + s->root->b * yl);
 }
-            } else if (mode == ARC) {
-                add_edge(&e[n++], sx, sy, sx_inner, sy_inner);
-                add_edge(&e[n++], lx, ly, lx_inner, ly_inner);
+
+// A pie. Can also be used to draw an arc with ugly sharp caps.
+void pie_init(clip_ellipse_state* s, int32_t a, int32_t b, int32_t w, float al, float ar) {
+  ellipse_init(&s->st, a, b, w);
+
+  s->head = NULL;
+  s->node_count = 0;
+
+  // line equations for pie sides
+  double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
+  double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
+
+  clip_node* lc = s->nodes + s->node_count++;
+  clip_node* rc = s->nodes + s->node_count++;
+  lc->l = lc->r = rc->l = rc->r = NULL;
+  lc->type = rc->type = CT_CLIP;
+  lc->a = -yl;
+  lc->b = xl;
+  lc->c = 0;
+  rc->a = yr;
+  rc->b = -xr;
+  rc->c = 0;
+        
+  s->root = s->nodes + s->node_count++;
+  s->root->l = lc;
+  s->root->r = rc;
+  s->root->type = ar - al < 180 ? CT_AND : CT_OR;
+}
+
+void clip_ellipse_free(clip_ellipse_state* s) {
+  while (s->head != NULL) {
+    event_list* t = s->head;
+    s->head = s->head->next;
+    free(t);
   }
 }
 
-        draw->polygon(im, n, e, ink, 0);
-
-        free(e);
+int8_t clip_ellipse_next(clip_ellipse_state* s, int32_t* ret_x0, int32_t* ret_y, int32_t* ret_x1) {
+  int32_t x0, y, x1;
+  while (s->head == NULL && ellipse_next(&s->st, &x0, &y, &x1) >= 0) {
+    if (clip_tree_do_clip(s->root, x0, y, x1, &s->head) < 0) {
+      return -2;
     }
-
+    s->y = y;
+  }
+  if (s->head != NULL) {
+    *ret_y = s->y;
+    event_list* t = s->head;
+    s->head = s->head->next;
+    *ret_x0 = t->x;
+    free(t);
+    t = s->head;
+    assert(t != NULL);
+    s->head = s->head->next;
+    *ret_x1 = t->x;
+    free(t);
     return 0;
   }
-
-int
-ImagingDrawArc(Imaging im, int x0, int y0, int x1, int y1,
-               float start, float end, const void* ink, int width, int op)
-{
-    return ellipse(im, x0, y0, x1, y1, start, end, ink, 0, width, ARC, op);
+  return -1;
 }
 
-int
-ImagingDrawChord(Imaging im, int x0, int y0, int x1, int y1,
-                 float start, float end, const void* ink, int fill,
+static int
+ellipseNew(Imaging im, int x0, int y0, int x1, int y1,
+        const void* ink_, int fill,
         int width, int op)
 {
-    return ellipse(im, x0, y0, x1, y1, start, end, ink, fill, width, CHORD, op);
+    DRAW* draw;
+    INT32 ink;
+    DRAWINIT();
+
+    int a = x1 - x0;
+    int b = y1 - y0;
+    if (a < 0 || b < 0) {
+      return 0;
+    }
+    if (fill) {
+      width = a + b;
+    }
+
+    ellipse_state st;
+    ellipse_init(&st, a, b, width);
+    int32_t X0, Y, X1;
+    while (ellipse_next(&st, &X0, &Y, &X1) != -1) {
+      draw->hline(im, x0 + (X0 + a) / 2, y0 + (Y + b) / 2, x0 + (X1 + a) / 2, ink);
+    }
+    return 0;
+}
+
+static int
+clipEllipseNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start, float end,
+        const void* ink_, int width, int op, clip_ellipse_init init)
+{
+    DRAW* draw;
+    INT32 ink;
+    DRAWINIT();
+
+    int a = x1 - x0;
+    int b = y1 - y0;
+    if (a < 0 || b < 0) {
+      return 0;
+    }
+
+    clip_ellipse_state st;
+    init(&st, a, b, width, start, end);
+    // debug_clip_tree(st.root, 0);
+    int32_t X0, Y, X1;
+    int next_code;
+    while ((next_code = clip_ellipse_next(&st, &X0, &Y, &X1)) >= 0) {
+      draw->hline(im, x0 + (X0 + a) / 2, y0 + (Y + b) / 2, x0 + (X1 + a) / 2, ink);
+    }
+    clip_ellipse_free(&st);
+    return next_code == -1 ? 0 : -1;
+}
+static int
+arcNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start, float end,
+        const void* ink_, int width, int op)
+{
+  return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, arc_init);
+}
+
+static int
+chordNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start, float end,
+        const void* ink_, int width, int op)
+{
+  return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, chord_init);
+}
+
+static int
+chordLineNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start, float end,
+        const void* ink_, int width, int op)
+{
+  return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, chord_line_init);
+}
+
+static int
+pieNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start, float end,
+        const void* ink_, int width, int op)
+{
+  return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, pie_init);
+}
+
+static int
+pieSideNew(Imaging im, int x0, int y0, int x1, int y1,
+        float start,
+        const void* ink_, int width, int op)
+{
+  return clipEllipseNew(im, x0, y0, x1, y1, start, 0, ink_, width, op, pie_side_init);
 }
 
 int
 ImagingDrawEllipse(Imaging im, int x0, int y0, int x1, int y1,
                    const void* ink, int fill, int width, int op)
 {
-    return ellipse(im, x0, y0, x1, y1, 0, 360, ink, fill, width, CHORD, op);
+    return ellipseNew(im, x0, y0, x1, y1, ink, fill, width, op);
 }
 
+int
+ImagingDrawArc(Imaging im, int x0, int y0, int x1, int y1,
+               float start, float end, const void* ink, int width, int op)
+{
+    normalize_angles(&start, &end);
+    if (start + 360 == end) {
+      return ImagingDrawEllipse(im, x0, y0, x1, y1, ink, 0, width, op);
+    }
+    if (start == end) {
+      return 0;
+    }
+    return arcNew(im, x0, y0, x1, y1, start, end, ink, width, op);
+}
+
+
+int
+ImagingDrawChord(Imaging im, int x0, int y0, int x1, int y1,
+                 float start, float end, const void* ink, int fill,
+                 int width, int op)
+{
+    normalize_angles(&start, &end);
+    if (start + 360 == end) {
+      return ImagingDrawEllipse(im, x0, y0, x1, y1, ink, fill, width, op);
+    }
+    if (start == end) {
+      return 0;
+    }
+    if (fill) {
+      return chordNew(im, x0, y0, x1, y1, start, end, ink, x1 - x0 + y1 - y0 + 1, op);
+    } else {
+      if (chordLineNew(im, x0, y0, x1, y1, start, end, ink, width, op)) {
+        return -1;
+      }
+      return chordNew(im, x0, y0, x1, y1, start, end, ink, width, op);
+    }
+}
+
+
 int
 ImagingDrawPieslice(Imaging im, int x0, int y0, int x1, int y1,
                     float start, float end, const void* ink, int fill,
                     int width, int op)
 {
-    return ellipse(im, x0, y0, x1, y1, start, end, ink, fill, width, PIESLICE, op);
+    normalize_angles(&start, &end);
+    if (start + 360 == end) {
+      return ellipseNew(im, x0, y0, x1, y1, ink, fill, width, op);
     }
+    if (start == end) {
+      return 0;
+    }
+    if (fill) {
+      return pieNew(im, x0, y0, x1, y1, start, end, ink, x1 + y1 - x0 - y0, op);
+    } else {
+      if (pieSideNew(im, x0, y0, x1, y1, start, ink, width, op)) {
+        return -1;
+      }
+      if (pieSideNew(im, x0, y0, x1, y1, end, ink, width, op)) {
+        return -1;
+      }
+      int xc = lround((x0 + x1 - width) / 2.0), yc = lround((y0 + y1 - width) / 2.0);
+      ellipseNew(im, xc, yc, xc + width - 1, yc + width - 1, ink, 1, 0, op);
+      return pieNew(im, x0, y0, x1, y1, start, end, ink, width, op);
+    }
+}
+
 
 /* -------------------------------------------------------------------- */