Pillow

1
/*
2
 * The Python Imaging Library.
3
 * $Id$
4
 *
5
 * a simple drawing package for the Imaging library
6
 *
7
 * history:
8
 * 1996-04-13 fl  Created.
9
 * 1996-04-30 fl  Added transforms and polygon support.
10
 * 1996-08-12 fl  Added filled polygons.
11
 * 1996-11-05 fl  Fixed float/int confusion in polygon filler
12
 * 1997-07-04 fl  Support 32-bit images (C++ would have been nice)
13
 * 1998-09-09 fl  Eliminated qsort casts; improved rectangle clipping
14
 * 1998-09-10 fl  Fixed fill rectangle to include lower edge (!)
15
 * 1998-12-29 fl  Added arc, chord, and pieslice primitives
16
 * 1999-01-10 fl  Added some level 2 ("arrow") stuff (experimental)
17
 * 1999-02-06 fl  Added bitmap primitive
18
 * 1999-07-26 fl  Eliminated a compiler warning
19
 * 1999-07-31 fl  Pass ink as void* instead of int
20
 * 2002-12-10 fl  Added experimental RGBA-on-RGB drawing
21
 * 2004-09-04 fl  Support simple wide lines (no joins)
22
 * 2005-05-25 fl  Fixed line width calculation
23
 *
24
 * Copyright (c) 1996-2006 by Fredrik Lundh
25
 * Copyright (c) 1997-2006 by Secret Labs AB.
26
 *
27
 * See the README file for information on usage and redistribution.
28
 */
29

30
/* FIXME: support fill/outline attribute for all filled shapes */
31
/* FIXME: support zero-winding fill */
32
/* FIXME: add drawing context, support affine transforms */
33
/* FIXME: support clip window (and mask?) */
34

35
#include "Imaging.h"
36

37
#include <math.h>
38
#include <stdint.h>
39

40
#define CEIL(v) (int)ceil(v)
41
#define FLOOR(v) ((v) >= 0.0 ? (int)(v) : (int)floor(v))
42

43
#define INK8(ink) (*(UINT8 *)ink)
44
#define INK16(ink) (*(UINT16 *)ink)
45

46
/*
47
 * Rounds around zero (up=away from zero, down=towards zero)
48
 * This guarantees that ROUND_UP|DOWN(f) == -ROUND_UP|DOWN(-f)
49
 */
50
#define ROUND_UP(f) ((int)((f) >= 0.0 ? floor((f) + 0.5F) : -floor(fabs(f) + 0.5F)))
51
#define ROUND_DOWN(f) ((int)((f) >= 0.0 ? ceil((f) - 0.5F) : -ceil(fabs(f) - 0.5F)))
52

53
/* -------------------------------------------------------------------- */
54
/* Primitives                                                           */
55
/* -------------------------------------------------------------------- */
56

57
typedef struct {
58
    /* edge descriptor for polygon engine */
59
    int d;
60
    int x0, y0;
61
    int xmin, ymin, xmax, ymax;
62
    float dx;
63
} Edge;
64

65
/* Type used in "polygon*" functions */
66
typedef void (*hline_handler)(Imaging, int, int, int, int);
67

68
static inline void
69
point8(Imaging im, int x, int y, int ink) {
70
    if (x >= 0 && x < im->xsize && y >= 0 && y < im->ysize) {
71
        if (strncmp(im->mode, "I;16", 4) == 0) {
72
#ifdef WORDS_BIGENDIAN
73
            im->image8[y][x * 2] = (UINT8)(ink >> 8);
74
            im->image8[y][x * 2 + 1] = (UINT8)ink;
75
#else
76
            im->image8[y][x * 2] = (UINT8)ink;
77
            im->image8[y][x * 2 + 1] = (UINT8)(ink >> 8);
78
#endif
79
        } else {
80
            im->image8[y][x] = (UINT8)ink;
81
        }
82
    }
83
}
84

85
static inline void
86
point32(Imaging im, int x, int y, int ink) {
87
    if (x >= 0 && x < im->xsize && y >= 0 && y < im->ysize) {
88
        im->image32[y][x] = ink;
89
    }
90
}
91

92
static inline void
93
point32rgba(Imaging im, int x, int y, int ink) {
94
    unsigned int tmp;
95

96
    if (x >= 0 && x < im->xsize && y >= 0 && y < im->ysize) {
97
        UINT8 *out = (UINT8 *)im->image[y] + x * 4;
98
        UINT8 *in = (UINT8 *)&ink;
99
        out[0] = BLEND(in[3], out[0], in[0], tmp);
100
        out[1] = BLEND(in[3], out[1], in[1], tmp);
101
        out[2] = BLEND(in[3], out[2], in[2], tmp);
102
    }
103
}
104

105
static inline void
106
hline8(Imaging im, int x0, int y0, int x1, int ink) {
107
    int pixelwidth;
108

109
    if (y0 >= 0 && y0 < im->ysize) {
110
        if (x0 < 0) {
111
            x0 = 0;
112
        } else if (x0 >= im->xsize) {
113
            return;
114
        }
115
        if (x1 < 0) {
116
            return;
117
        } else if (x1 >= im->xsize) {
118
            x1 = im->xsize - 1;
119
        }
120
        if (x0 <= x1) {
121
            pixelwidth = strncmp(im->mode, "I;16", 4) == 0 ? 2 : 1;
122
            memset(
123
                im->image8[y0] + x0 * pixelwidth, (UINT8)ink, (x1 - x0 + 1) * pixelwidth
124
            );
125
        }
126
    }
127
}
128

129
static inline void
130
hline32(Imaging im, int x0, int y0, int x1, int ink) {
131
    INT32 *p;
132

133
    if (y0 >= 0 && y0 < im->ysize) {
134
        if (x0 < 0) {
135
            x0 = 0;
136
        } else if (x0 >= im->xsize) {
137
            return;
138
        }
139
        if (x1 < 0) {
140
            return;
141
        } else if (x1 >= im->xsize) {
142
            x1 = im->xsize - 1;
143
        }
144
        p = im->image32[y0];
145
        while (x0 <= x1) {
146
            p[x0++] = ink;
147
        }
148
    }
149
}
150

151
static inline void
152
hline32rgba(Imaging im, int x0, int y0, int x1, int ink) {
153
    unsigned int tmp;
154

155
    if (y0 >= 0 && y0 < im->ysize) {
156
        if (x0 < 0) {
157
            x0 = 0;
158
        } else if (x0 >= im->xsize) {
159
            return;
160
        }
161
        if (x1 < 0) {
162
            return;
163
        } else if (x1 >= im->xsize) {
164
            x1 = im->xsize - 1;
165
        }
166
        if (x0 <= x1) {
167
            UINT8 *out = (UINT8 *)im->image[y0] + x0 * 4;
168
            UINT8 *in = (UINT8 *)&ink;
169
            while (x0 <= x1) {
170
                out[0] = BLEND(in[3], out[0], in[0], tmp);
171
                out[1] = BLEND(in[3], out[1], in[1], tmp);
172
                out[2] = BLEND(in[3], out[2], in[2], tmp);
173
                x0++;
174
                out += 4;
175
            }
176
        }
177
    }
178
}
179

180
static inline void
181
line8(Imaging im, int x0, int y0, int x1, int y1, int ink) {
182
    int i, n, e;
183
    int dx, dy;
184
    int xs, ys;
185

186
    /* normalize coordinates */
187
    dx = x1 - x0;
188
    if (dx < 0) {
189
        dx = -dx, xs = -1;
190
    } else {
191
        xs = 1;
192
    }
193
    dy = y1 - y0;
194
    if (dy < 0) {
195
        dy = -dy, ys = -1;
196
    } else {
197
        ys = 1;
198
    }
199

200
    n = (dx > dy) ? dx : dy;
201

202
    if (dx == 0) {
203
        /* vertical */
204
        for (i = 0; i < dy; i++) {
205
            point8(im, x0, y0, ink);
206
            y0 += ys;
207
        }
208

209
    } else if (dy == 0) {
210
        /* horizontal */
211
        for (i = 0; i < dx; i++) {
212
            point8(im, x0, y0, ink);
213
            x0 += xs;
214
        }
215

216
    } else if (dx > dy) {
217
        /* bresenham, horizontal slope */
218
        n = dx;
219
        dy += dy;
220
        e = dy - dx;
221
        dx += dx;
222

223
        for (i = 0; i < n; i++) {
224
            point8(im, x0, y0, ink);
225
            if (e >= 0) {
226
                y0 += ys;
227
                e -= dx;
228
            }
229
            e += dy;
230
            x0 += xs;
231
        }
232

233
    } else {
234
        /* bresenham, vertical slope */
235
        n = dy;
236
        dx += dx;
237
        e = dx - dy;
238
        dy += dy;
239

240
        for (i = 0; i < n; i++) {
241
            point8(im, x0, y0, ink);
242
            if (e >= 0) {
243
                x0 += xs;
244
                e -= dy;
245
            }
246
            e += dx;
247
            y0 += ys;
248
        }
249
    }
250
}
251

252
static inline void
253
line32(Imaging im, int x0, int y0, int x1, int y1, int ink) {
254
    int i, n, e;
255
    int dx, dy;
256
    int xs, ys;
257

258
    /* normalize coordinates */
259
    dx = x1 - x0;
260
    if (dx < 0) {
261
        dx = -dx, xs = -1;
262
    } else {
263
        xs = 1;
264
    }
265
    dy = y1 - y0;
266
    if (dy < 0) {
267
        dy = -dy, ys = -1;
268
    } else {
269
        ys = 1;
270
    }
271

272
    n = (dx > dy) ? dx : dy;
273

274
    if (dx == 0) {
275
        /* vertical */
276
        for (i = 0; i < dy; i++) {
277
            point32(im, x0, y0, ink);
278
            y0 += ys;
279
        }
280

281
    } else if (dy == 0) {
282
        /* horizontal */
283
        for (i = 0; i < dx; i++) {
284
            point32(im, x0, y0, ink);
285
            x0 += xs;
286
        }
287

288
    } else if (dx > dy) {
289
        /* bresenham, horizontal slope */
290
        n = dx;
291
        dy += dy;
292
        e = dy - dx;
293
        dx += dx;
294

295
        for (i = 0; i < n; i++) {
296
            point32(im, x0, y0, ink);
297
            if (e >= 0) {
298
                y0 += ys;
299
                e -= dx;
300
            }
301
            e += dy;
302
            x0 += xs;
303
        }
304

305
    } else {
306
        /* bresenham, vertical slope */
307
        n = dy;
308
        dx += dx;
309
        e = dx - dy;
310
        dy += dy;
311

312
        for (i = 0; i < n; i++) {
313
            point32(im, x0, y0, ink);
314
            if (e >= 0) {
315
                x0 += xs;
316
                e -= dy;
317
            }
318
            e += dx;
319
            y0 += ys;
320
        }
321
    }
322
}
323

324
static inline void
325
line32rgba(Imaging im, int x0, int y0, int x1, int y1, int ink) {
326
    int i, n, e;
327
    int dx, dy;
328
    int xs, ys;
329

330
    /* normalize coordinates */
331
    dx = x1 - x0;
332
    if (dx < 0) {
333
        dx = -dx, xs = -1;
334
    } else {
335
        xs = 1;
336
    }
337
    dy = y1 - y0;
338
    if (dy < 0) {
339
        dy = -dy, ys = -1;
340
    } else {
341
        ys = 1;
342
    }
343

344
    n = (dx > dy) ? dx : dy;
345

346
    if (dx == 0) {
347
        /* vertical */
348
        for (i = 0; i < dy; i++) {
349
            point32rgba(im, x0, y0, ink);
350
            y0 += ys;
351
        }
352

353
    } else if (dy == 0) {
354
        /* horizontal */
355
        for (i = 0; i < dx; i++) {
356
            point32rgba(im, x0, y0, ink);
357
            x0 += xs;
358
        }
359

360
    } else if (dx > dy) {
361
        /* bresenham, horizontal slope */
362
        n = dx;
363
        dy += dy;
364
        e = dy - dx;
365
        dx += dx;
366

367
        for (i = 0; i < n; i++) {
368
            point32rgba(im, x0, y0, ink);
369
            if (e >= 0) {
370
                y0 += ys;
371
                e -= dx;
372
            }
373
            e += dy;
374
            x0 += xs;
375
        }
376

377
    } else {
378
        /* bresenham, vertical slope */
379
        n = dy;
380
        dx += dx;
381
        e = dx - dy;
382
        dy += dy;
383

384
        for (i = 0; i < n; i++) {
385
            point32rgba(im, x0, y0, ink);
386
            if (e >= 0) {
387
                x0 += xs;
388
                e -= dy;
389
            }
390
            e += dx;
391
            y0 += ys;
392
        }
393
    }
394
}
395

396
static int
397
x_cmp(const void *x0, const void *x1) {
398
    float diff = *((float *)x0) - *((float *)x1);
399
    if (diff < 0) {
400
        return -1;
401
    } else if (diff > 0) {
402
        return 1;
403
    } else {
404
        return 0;
405
    }
406
}
407

408
static void
409
draw_horizontal_lines(
410
    Imaging im, int n, Edge *e, int ink, int *x_pos, int y, hline_handler hline
411
) {
412
    int i;
413
    for (i = 0; i < n; i++) {
414
        if (e[i].ymin == y && e[i].ymin == e[i].ymax) {
415
            int xmax;
416
            int xmin = e[i].xmin;
417
            if (*x_pos != -1 && *x_pos < xmin) {
418
                // Line would be after the current position
419
                continue;
420
            }
421

422
            xmax = e[i].xmax;
423
            if (*x_pos > xmin) {
424
                // Line would be partway through x_pos, so increase the starting point
425
                xmin = *x_pos;
426
                if (xmax < xmin) {
427
                    // Line would now end before it started
428
                    continue;
429
                }
430
            }
431

432
            (*hline)(im, xmin, e[i].ymin, xmax, ink);
433
            *x_pos = xmax + 1;
434
        }
435
    }
436
}
437

438
/*
439
 * Filled polygon draw function using scan line algorithm.
440
 */
441
static inline int
442
polygon_generic(
443
    Imaging im, int n, Edge *e, int ink, int eofill, hline_handler hline, int hasAlpha
444
) {
445
    Edge **edge_table;
446
    float *xx;
447
    int edge_count = 0;
448
    int ymin = im->ysize - 1;
449
    int ymax = 0;
450
    int i, j, k;
451
    float adjacent_line_x, adjacent_line_x_other_edge;
452

453
    if (n <= 0) {
454
        return 0;
455
    }
456

457
    /* Initialize the edge table and find polygon boundaries */
458
    /* malloc check ok, using calloc */
459
    edge_table = calloc(n, sizeof(Edge *));
460
    if (!edge_table) {
461
        return -1;
462
    }
463

464
    for (i = 0; i < n; i++) {
465
        if (ymin > e[i].ymin) {
466
            ymin = e[i].ymin;
467
        }
468
        if (ymax < e[i].ymax) {
469
            ymax = e[i].ymax;
470
        }
471
        if (e[i].ymin == e[i].ymax) {
472
            if (hasAlpha != 1) {
473
                (*hline)(im, e[i].xmin, e[i].ymin, e[i].xmax, ink);
474
            }
475
            continue;
476
        }
477
        edge_table[edge_count++] = (e + i);
478
    }
479
    if (ymin < 0) {
480
        ymin = 0;
481
    }
482
    if (ymax > im->ysize) {
483
        ymax = im->ysize;
484
    }
485

486
    /* Process the edge table with a scan line searching for intersections */
487
    /* malloc check ok, using calloc */
488
    xx = calloc(edge_count * 2, sizeof(float));
489
    if (!xx) {
490
        free(edge_table);
491
        return -1;
492
    }
493
    for (; ymin <= ymax; ymin++) {
494
        j = 0;
495
        for (i = 0; i < edge_count; i++) {
496
            Edge *current = edge_table[i];
497
            if (ymin >= current->ymin && ymin <= current->ymax) {
498
                xx[j++] = (ymin - current->y0) * current->dx + current->x0;
499

500
                if (ymin == current->ymax && ymin < ymax) {
501
                    // Needed to draw consistent polygons
502
                    xx[j] = xx[j - 1];
503
                    j++;
504
                } else if (current->dx != 0 && roundf(xx[j - 1]) == xx[j - 1]) {
505
                    // Connect discontiguous corners
506
                    for (k = 0; k < i; k++) {
507
                        Edge *other_edge = edge_table[k];
508
                        if ((current->dx > 0 && other_edge->dx <= 0) ||
509
                            (current->dx < 0 && other_edge->dx >= 0)) {
510
                            continue;
511
                        }
512
                        // Check if the two edges join to make a corner
513
                        if (((ymin == current->ymin && ymin == other_edge->ymin) ||
514
                             (ymin == current->ymax && ymin == other_edge->ymax)) &&
515
                            xx[j - 1] == (ymin - other_edge->y0) * other_edge->dx +
516
                                             other_edge->x0) {
517
                            // Determine points from the edges on the next row
518
                            // Or if this is the last row, check the previous row
519
                            int offset = ymin == ymax ? -1 : 1;
520
                            adjacent_line_x =
521
                                (ymin + offset - current->y0) * current->dx +
522
                                current->x0;
523
                            adjacent_line_x_other_edge =
524
                                (ymin + offset - other_edge->y0) * other_edge->dx +
525
                                other_edge->x0;
526
                            if (ymin == current->ymax) {
527
                                if (current->dx > 0) {
528
                                    xx[k] =
529
                                        fmax(
530
                                            adjacent_line_x, adjacent_line_x_other_edge
531
                                        ) +
532
                                        1;
533
                                } else {
534
                                    xx[k] =
535
                                        fmin(
536
                                            adjacent_line_x, adjacent_line_x_other_edge
537
                                        ) -
538
                                        1;
539
                                }
540
                            } else {
541
                                if (current->dx > 0) {
542
                                    xx[k] = fmin(
543
                                        adjacent_line_x, adjacent_line_x_other_edge
544
                                    );
545
                                } else {
546
                                    xx[k] =
547
                                        fmax(
548
                                            adjacent_line_x, adjacent_line_x_other_edge
549
                                        ) +
550
                                        1;
551
                                }
552
                            }
553
                            break;
554
                        }
555
                    }
556
                }
557
            }
558
        }
559
        qsort(xx, j, sizeof(float), x_cmp);
560
        if (hasAlpha == 1) {
561
            int x_pos = j == 0 ? -1 : 0;
562
            for (i = 1; i < j; i += 2) {
563
                int x_end = ROUND_DOWN(xx[i]);
564
                if (x_end < x_pos) {
565
                    // Line would be before the current position
566
                    continue;
567
                }
568
                draw_horizontal_lines(im, n, e, ink, &x_pos, ymin, hline);
569
                if (x_end < x_pos) {
570
                    // Line would be before the current position
571
                    continue;
572
                }
573

574
                int x_start = ROUND_UP(xx[i - 1]);
575
                if (x_pos > x_start) {
576
                    // Line would be partway through x_pos, so increase the starting
577
                    // point
578
                    x_start = x_pos;
579
                    if (x_end < x_start) {
580
                        // Line would now end before it started
581
                        continue;
582
                    }
583
                }
584
                (*hline)(im, x_start, ymin, x_end, ink);
585
                x_pos = x_end + 1;
586
            }
587
            draw_horizontal_lines(im, n, e, ink, &x_pos, ymin, hline);
588
        } else {
589
            for (i = 1; i < j; i += 2) {
590
                (*hline)(im, ROUND_UP(xx[i - 1]), ymin, ROUND_DOWN(xx[i]), ink);
591
            }
592
        }
593
    }
594

595
    free(xx);
596
    free(edge_table);
597
    return 0;
598
}
599

600
static inline int
601
polygon8(Imaging im, int n, Edge *e, int ink, int eofill) {
602
    return polygon_generic(im, n, e, ink, eofill, hline8, 0);
603
}
604

605
static inline int
606
polygon32(Imaging im, int n, Edge *e, int ink, int eofill) {
607
    return polygon_generic(im, n, e, ink, eofill, hline32, 0);
608
}
609

610
static inline int
611
polygon32rgba(Imaging im, int n, Edge *e, int ink, int eofill) {
612
    return polygon_generic(im, n, e, ink, eofill, hline32rgba, 1);
613
}
614

615
static inline void
616
add_edge(Edge *e, int x0, int y0, int x1, int y1) {
617
    /* printf("edge %d %d %d %d\n", x0, y0, x1, y1); */
618

619
    if (x0 <= x1) {
620
        e->xmin = x0, e->xmax = x1;
621
    } else {
622
        e->xmin = x1, e->xmax = x0;
623
    }
624

625
    if (y0 <= y1) {
626
        e->ymin = y0, e->ymax = y1;
627
    } else {
628
        e->ymin = y1, e->ymax = y0;
629
    }
630

631
    if (y0 == y1) {
632
        e->d = 0;
633
        e->dx = 0.0;
634
    } else {
635
        e->dx = ((float)(x1 - x0)) / (y1 - y0);
636
        if (y0 == e->ymin) {
637
            e->d = 1;
638
        } else {
639
            e->d = -1;
640
        }
641
    }
642

643
    e->x0 = x0;
644
    e->y0 = y0;
645
}
646

647
typedef struct {
648
    void (*point)(Imaging im, int x, int y, int ink);
649
    void (*hline)(Imaging im, int x0, int y0, int x1, int ink);
650
    void (*line)(Imaging im, int x0, int y0, int x1, int y1, int ink);
651
    int (*polygon)(Imaging im, int n, Edge *e, int ink, int eofill);
652
} DRAW;
653

654
DRAW draw8 = {point8, hline8, line8, polygon8};
655
DRAW draw32 = {point32, hline32, line32, polygon32};
656
DRAW draw32rgba = {point32rgba, hline32rgba, line32rgba, polygon32rgba};
657

658
/* -------------------------------------------------------------------- */
659
/* Interface                                                            */
660
/* -------------------------------------------------------------------- */
661

662
#define DRAWINIT()                               \
663
    if (im->image8) {                            \
664
        draw = &draw8;                           \
665
        if (strncmp(im->mode, "I;16", 4) == 0) { \
666
            ink = INK16(ink_);                   \
667
        } else {                                 \
668
            ink = INK8(ink_);                    \
669
        }                                        \
670
    } else {                                     \
671
        draw = (op) ? &draw32rgba : &draw32;     \
672
        memcpy(&ink, ink_, sizeof(ink));         \
673
    }
674

675
int
676
ImagingDrawPoint(Imaging im, int x0, int y0, const void *ink_, int op) {
677
    DRAW *draw;
678
    INT32 ink;
679

680
    DRAWINIT();
681

682
    draw->point(im, x0, y0, ink);
683

684
    return 0;
685
}
686

687
int
688
ImagingDrawLine(Imaging im, int x0, int y0, int x1, int y1, const void *ink_, int op) {
689
    DRAW *draw;
690
    INT32 ink;
691

692
    DRAWINIT();
693

694
    draw->line(im, x0, y0, x1, y1, ink);
695

696
    return 0;
697
}
698

699
int
700
ImagingDrawWideLine(
701
    Imaging im, int x0, int y0, int x1, int y1, const void *ink_, int width, int op
702
) {
703
    DRAW *draw;
704
    INT32 ink;
705
    int dx, dy;
706
    double big_hypotenuse, small_hypotenuse, ratio_max, ratio_min;
707
    int dxmin, dxmax, dymin, dymax;
708
    Edge e[4];
709

710
    DRAWINIT();
711

712
    dx = x1 - x0;
713
    dy = y1 - y0;
714
    if (dx == 0 && dy == 0) {
715
        draw->point(im, x0, y0, ink);
716
        return 0;
717
    }
718

719
    big_hypotenuse = hypot(dx, dy);
720
    small_hypotenuse = (width - 1) / 2.0;
721
    ratio_max = ROUND_UP(small_hypotenuse) / big_hypotenuse;
722
    ratio_min = ROUND_DOWN(small_hypotenuse) / big_hypotenuse;
723

724
    dxmin = ROUND_DOWN(ratio_min * dy);
725
    dxmax = ROUND_DOWN(ratio_max * dy);
726
    dymin = ROUND_DOWN(ratio_min * dx);
727
    dymax = ROUND_DOWN(ratio_max * dx);
728
    {
729
        int vertices[4][2] = {
730
            {x0 - dxmin, y0 + dymax},
731
            {x1 - dxmin, y1 + dymax},
732
            {x1 + dxmax, y1 - dymin},
733
            {x0 + dxmax, y0 - dymin}
734
        };
735

736
        add_edge(e + 0, vertices[0][0], vertices[0][1], vertices[1][0], vertices[1][1]);
737
        add_edge(e + 1, vertices[1][0], vertices[1][1], vertices[2][0], vertices[2][1]);
738
        add_edge(e + 2, vertices[2][0], vertices[2][1], vertices[3][0], vertices[3][1]);
739
        add_edge(e + 3, vertices[3][0], vertices[3][1], vertices[0][0], vertices[0][1]);
740

741
        draw->polygon(im, 4, e, ink, 0);
742
    }
743
    return 0;
744
}
745

746
int
747
ImagingDrawRectangle(
748
    Imaging im,
749
    int x0,
750
    int y0,
751
    int x1,
752
    int y1,
753
    const void *ink_,
754
    int fill,
755
    int width,
756
    int op
757
) {
758
    int i;
759
    int y;
760
    int tmp;
761
    DRAW *draw;
762
    INT32 ink;
763

764
    DRAWINIT();
765

766
    if (y0 > y1) {
767
        tmp = y0, y0 = y1, y1 = tmp;
768
    }
769

770
    if (fill) {
771
        if (y0 < 0) {
772
            y0 = 0;
773
        } else if (y0 >= im->ysize) {
774
            return 0;
775
        }
776

777
        if (y1 < 0) {
778
            return 0;
779
        } else if (y1 > im->ysize) {
780
            y1 = im->ysize;
781
        }
782

783
        for (y = y0; y <= y1; y++) {
784
            draw->hline(im, x0, y, x1, ink);
785
        }
786

787
    } else {
788
        /* outline */
789
        if (width == 0) {
790
            width = 1;
791
        }
792
        for (i = 0; i < width; i++) {
793
            draw->hline(im, x0, y0 + i, x1, ink);
794
            draw->hline(im, x0, y1 - i, x1, ink);
795
            draw->line(im, x1 - i, y0 + width, x1 - i, y1 - width + 1, ink);
796
            draw->line(im, x0 + i, y0 + width, x0 + i, y1 - width + 1, ink);
797
        }
798
    }
799

800
    return 0;
801
}
802

803
int
804
ImagingDrawPolygon(
805
    Imaging im, int count, int *xy, const void *ink_, int fill, int width, int op
806
) {
807
    int i, n, x0, y0, x1, y1;
808
    DRAW *draw;
809
    INT32 ink;
810

811
    if (count <= 0) {
812
        return 0;
813
    }
814

815
    DRAWINIT();
816

817
    if (fill) {
818
        /* Build edge list */
819
        /* malloc check ok, using calloc */
820
        Edge *e = calloc(count, sizeof(Edge));
821
        if (!e) {
822
            (void)ImagingError_MemoryError();
823
            return -1;
824
        }
825
        for (i = n = 0; i < count - 1; i++) {
826
            x0 = xy[i * 2];
827
            y0 = xy[i * 2 + 1];
828
            x1 = xy[i * 2 + 2];
829
            y1 = xy[i * 2 + 3];
830
            if (y0 == y1 && i != 0 && y0 == xy[i * 2 - 1]) {
831
                // This is a horizontal line,
832
                // that immediately follows another horizontal line
833
                Edge *last_e = &e[n - 1];
834
                if (x1 > x0 && x0 > xy[i * 2 - 2]) {
835
                    // They are both increasing in x
836
                    last_e->xmax = x1;
837
                    continue;
838
                } else if (x1 < x0 && x0 < xy[i * 2 - 2]) {
839
                    // They are both decreasing in x
840
                    last_e->xmin = x1;
841
                    continue;
842
                }
843
            }
844
            add_edge(&e[n++], x0, y0, x1, y1);
845
        }
846
        if (xy[i * 2] != xy[0] || xy[i * 2 + 1] != xy[1]) {
847
            add_edge(&e[n++], xy[i * 2], xy[i * 2 + 1], xy[0], xy[1]);
848
        }
849
        draw->polygon(im, n, e, ink, 0);
850
        free(e);
851

852
    } else {
853
        /* Outline */
854
        if (width == 1) {
855
            for (i = 0; i < count - 1; i++) {
856
                draw->line(
857
                    im, xy[i * 2], xy[i * 2 + 1], xy[i * 2 + 2], xy[i * 2 + 3], ink
858
                );
859
            }
860
            draw->line(im, xy[i * 2], xy[i * 2 + 1], xy[0], xy[1], ink);
861
        } else {
862
            for (i = 0; i < count - 1; i++) {
863
                ImagingDrawWideLine(
864
                    im,
865
                    xy[i * 2],
866
                    xy[i * 2 + 1],
867
                    xy[i * 2 + 2],
868
                    xy[i * 2 + 3],
869
                    ink_,
870
                    width,
871
                    op
872
                );
873
            }
874
            ImagingDrawWideLine(
875
                im, xy[i * 2], xy[i * 2 + 1], xy[0], xy[1], ink_, width, op
876
            );
877
        }
878
    }
879

880
    return 0;
881
}
882

883
int
884
ImagingDrawBitmap(Imaging im, int x0, int y0, Imaging bitmap, const void *ink, int op) {
885
    return ImagingFill2(
886
        im, ink, bitmap, x0, y0, x0 + bitmap->xsize, y0 + bitmap->ysize
887
    );
888
}
889

890
/* -------------------------------------------------------------------- */
891
/* standard shapes */
892

893
// Imagine 2D plane and ellipse with center in (0, 0) and semi-major axes a and b.
894
// Then quarter_* stuff approximates its top right quarter (x, y >= 0) with integer
895
// points from set {(2x+x0, 2y+y0) | x,y in Z} where x0, y0 are from {0, 1} and
896
// are such that point (a, b) is in the set.
897

898
typedef struct {
899
    int32_t a, b, cx, cy, ex, ey;
900
    int64_t a2, b2, a2b2;
901
    int8_t finished;
902
} quarter_state;
903

904
void
905
quarter_init(quarter_state *s, int32_t a, int32_t b) {
906
    if (a < 0 || b < 0) {
907
        s->finished = 1;
908
    } else {
909
        s->a = a;
910
        s->b = b;
911
        s->cx = a;
912
        s->cy = b % 2;
913
        s->ex = a % 2;
914
        s->ey = b;
915
        s->a2 = a * a;
916
        s->b2 = b * b;
917
        s->a2b2 = s->a2 * s->b2;
918
        s->finished = 0;
919
    }
920
}
921

922
// deviation of the point from ellipse curve, basically a substitution
923
// of the point into the ellipse equation
924
int64_t
925
quarter_delta(quarter_state *s, int64_t x, int64_t y) {
926
    return llabs(s->a2 * y * y + s->b2 * x * x - s->a2b2);
927
}
928

929
int8_t
930
quarter_next(quarter_state *s, int32_t *ret_x, int32_t *ret_y) {
931
    if (s->finished) {
932
        return -1;
933
    }
934
    *ret_x = s->cx;
935
    *ret_y = s->cy;
936
    if (s->cx == s->ex && s->cy == s->ey) {
937
        s->finished = 1;
938
    } else {
939
        // Bresenham's algorithm, possible optimization: only consider 2 of 3
940
        // next points depending on current slope
941
        int32_t nx = s->cx;
942
        int32_t ny = s->cy + 2;
943
        int64_t ndelta = quarter_delta(s, nx, ny);
944
        if (nx > 1) {
945
            int64_t newdelta = quarter_delta(s, s->cx - 2, s->cy + 2);
946
            if (ndelta > newdelta) {
947
                nx = s->cx - 2;
948
                ny = s->cy + 2;
949
                ndelta = newdelta;
950
            }
951
            newdelta = quarter_delta(s, s->cx - 2, s->cy);
952
            if (ndelta > newdelta) {
953
                nx = s->cx - 2;
954
                ny = s->cy;
955
            }
956
        }
957
        s->cx = nx;
958
        s->cy = ny;
959
    }
960
    return 0;
961
}
962

963
// quarter_* stuff can "draw" a quarter of an ellipse with thickness 1, great.
964
// Now we use ellipse_* stuff to join all four quarters of two different sized
965
// ellipses and receive horizontal segments of a complete ellipse with
966
// specified thickness.
967
//
968
// Still using integer grid with step 2 at this point (like in quarter_*)
969
// to ease angle clipping in future.
970

971
typedef struct {
972
    quarter_state st_o, st_i;
973
    int32_t py, pl, pr;
974
    int32_t cy[4], cl[4], cr[4];
975
    int8_t bufcnt;
976
    int8_t finished;
977
    int8_t leftmost;
978
} ellipse_state;
979

980
void
981
ellipse_init(ellipse_state *s, int32_t a, int32_t b, int32_t w) {
982
    s->bufcnt = 0;
983
    s->leftmost = a % 2;
984
    quarter_init(&s->st_o, a, b);
985
    if (w < 1 || quarter_next(&s->st_o, &s->pr, &s->py) == -1) {
986
        s->finished = 1;
987
    } else {
988
        s->finished = 0;
989
        quarter_init(&s->st_i, a - 2 * (w - 1), b - 2 * (w - 1));
990
        s->pl = s->leftmost;
991
    }
992
}
993

994
int8_t
995
ellipse_next(ellipse_state *s, int32_t *ret_x0, int32_t *ret_y, int32_t *ret_x1) {
996
    if (s->bufcnt == 0) {
997
        if (s->finished) {
998
            return -1;
999
        }
1000
        int32_t y = s->py;
1001
        int32_t l = s->pl;
1002
        int32_t r = s->pr;
1003
        int32_t cx = 0, cy = 0;
1004
        int8_t next_ret;
1005
        while ((next_ret = quarter_next(&s->st_o, &cx, &cy)) != -1 && cy <= y) {
1006
        }
1007
        if (next_ret == -1) {
1008
            s->finished = 1;
1009
        } else {
1010
            s->pr = cx;
1011
            s->py = cy;
1012
        }
1013
        while ((next_ret = quarter_next(&s->st_i, &cx, &cy)) != -1 && cy <= y) {
1014
            l = cx;
1015
        }
1016
        s->pl = next_ret == -1 ? s->leftmost : cx;
1017

1018
        if ((l > 0 || l < r) && y > 0) {
1019
            s->cl[s->bufcnt] = l == 0 ? 2 : l;
1020
            s->cy[s->bufcnt] = y;
1021
            s->cr[s->bufcnt] = r;
1022
            ++s->bufcnt;
1023
        }
1024
        if (y > 0) {
1025
            s->cl[s->bufcnt] = -r;
1026
            s->cy[s->bufcnt] = y;
1027
            s->cr[s->bufcnt] = -l;
1028
            ++s->bufcnt;
1029
        }
1030
        if (l > 0 || l < r) {
1031
            s->cl[s->bufcnt] = l == 0 ? 2 : l;
1032
            s->cy[s->bufcnt] = -y;
1033
            s->cr[s->bufcnt] = r;
1034
            ++s->bufcnt;
1035
        }
1036
        s->cl[s->bufcnt] = -r;
1037
        s->cy[s->bufcnt] = -y;
1038
        s->cr[s->bufcnt] = -l;
1039
        ++s->bufcnt;
1040
    }
1041
    --s->bufcnt;
1042
    *ret_x0 = s->cl[s->bufcnt];
1043
    *ret_y = s->cy[s->bufcnt];
1044
    *ret_x1 = s->cr[s->bufcnt];
1045
    return 0;
1046
}
1047

1048
// Clipping tree consists of half-plane clipping nodes and combining nodes.
1049
// We can throw a horizontal segment in such a tree and collect an ordered set
1050
// of resulting disjoint clipped segments organized into a sorted linked list
1051
// of their end points.
1052
typedef enum {
1053
    CT_AND,  // intersection
1054
    CT_OR,   // union
1055
    CT_CLIP  // half-plane clipping
1056
} clip_type;
1057

1058
typedef struct clip_node {
1059
    clip_type type;
1060
    double a, b, c;       // half-plane coeffs, only used in clipping nodes
1061
    struct clip_node *l;  // child pointers, are only non-NULL in combining nodes
1062
    struct clip_node *r;
1063
} clip_node;
1064

1065
// Linked list for the ends of the clipped horizontal segments.
1066
// Since the segment is always horizontal, we don't need to store Y coordinate.
1067
typedef struct event_list {
1068
    int32_t x;
1069
    int8_t type;  // used internally, 1 for the left end (smaller X), -1 for the
1070
                  // right end; pointless in output since the output segments
1071
                  // are disjoint, therefore the types would always come in pairs
1072
                  // and interchange (1 -1 1 -1 ...)
1073
    struct event_list *next;
1074
} event_list;
1075

1076
// Mirrors all the clipping nodes of the tree relative to the y = x line.
1077
void
1078
clip_tree_transpose(clip_node *root) {
1079
    if (root != NULL) {
1080
        if (root->type == CT_CLIP) {
1081
            double t = root->a;
1082
            root->a = root->b;
1083
            root->b = t;
1084
        }
1085
        clip_tree_transpose(root->l);
1086
        clip_tree_transpose(root->r);
1087
    }
1088
}
1089

1090
// Outputs a sequence of open-close events (types -1 and 1) for
1091
// non-intersecting segments sorted by X coordinate.
1092
// Combining nodes (AND, OR) may also accept sequences for intersecting
1093
// segments, i.e. something like correct bracket sequences.
1094
int
1095
clip_tree_do_clip(
1096
    clip_node *root, int32_t x0, int32_t y, int32_t x1, event_list **ret
1097
) {
1098
    if (root == NULL) {
1099
        event_list *start = malloc(sizeof(event_list));
1100
        if (!start) {
1101
            ImagingError_MemoryError();
1102
            return -1;
1103
        }
1104
        event_list *end = malloc(sizeof(event_list));
1105
        if (!end) {
1106
            free(start);
1107
            ImagingError_MemoryError();
1108
            return -1;
1109
        }
1110
        start->x = x0;
1111
        start->type = 1;
1112
        start->next = end;
1113
        end->x = x1;
1114
        end->type = -1;
1115
        end->next = NULL;
1116
        *ret = start;
1117
        return 0;
1118
    }
1119
    if (root->type == CT_CLIP) {
1120
        double eps = 1e-9;
1121
        double A = root->a;
1122
        double B = root->b;
1123
        double C = root->c;
1124
        if (fabs(A) < eps) {
1125
            if (B * y + C < -eps) {
1126
                x0 = 1;
1127
                x1 = 0;
1128
            }
1129
        } else {
1130
            // X of intersection
1131
            double ix = -(B * y + C) / A;
1132
            if (A * x0 + B * y + C < eps) {
1133
                x0 = lround(fmax(x0, ix));
1134
            }
1135
            if (A * x1 + B * y + C < eps) {
1136
                x1 = lround(fmin(x1, ix));
1137
            }
1138
        }
1139
        if (x0 <= x1) {
1140
            event_list *start = malloc(sizeof(event_list));
1141
            if (!start) {
1142
                ImagingError_MemoryError();
1143
                return -1;
1144
            }
1145
            event_list *end = malloc(sizeof(event_list));
1146
            if (!end) {
1147
                free(start);
1148
                ImagingError_MemoryError();
1149
                return -1;
1150
            }
1151
            start->x = x0;
1152
            start->type = 1;
1153
            start->next = end;
1154
            end->x = x1;
1155
            end->type = -1;
1156
            end->next = NULL;
1157
            *ret = start;
1158
        } else {
1159
            *ret = NULL;
1160
        }
1161
        return 0;
1162
    }
1163
    if (root->type == CT_OR || root->type == CT_AND) {
1164
        event_list *l1;
1165
        event_list *l2;
1166
        if (clip_tree_do_clip(root->l, x0, y, x1, &l1) < 0) {
1167
            return -1;
1168
        }
1169
        if (clip_tree_do_clip(root->r, x0, y, x1, &l2) < 0) {
1170
            while (l1) {
1171
                l2 = l1->next;
1172
                free(l1);
1173
                l1 = l2;
1174
            }
1175
            return -1;
1176
        }
1177
        *ret = NULL;
1178
        event_list *tail = NULL;
1179
        int32_t k1 = 0;
1180
        int32_t k2 = 0;
1181
        while (l1 != NULL || l2 != NULL) {
1182
            event_list *t;
1183
            if (l2 == NULL ||
1184
                (l1 != NULL &&
1185
                 (l1->x < l2->x || (l1->x == l2->x && l1->type > l2->type)))) {
1186
                t = l1;
1187
                k1 += t->type;
1188
                assert(k1 >= 0);
1189
                l1 = l1->next;
1190
            } else {
1191
                t = l2;
1192
                k2 += t->type;
1193
                assert(k2 >= 0);
1194
                l2 = l2->next;
1195
            }
1196
            t->next = NULL;
1197
            if ((root->type == CT_OR &&
1198
                 ((t->type == 1 && (tail == NULL || tail->type == -1)) ||
1199
                  (t->type == -1 && k1 == 0 && k2 == 0))) ||
1200
                (root->type == CT_AND &&
1201
                 ((t->type == 1 && (tail == NULL || tail->type == -1) && k1 > 0 &&
1202
                   k2 > 0) ||
1203
                  (t->type == -1 && tail != NULL && tail->type == 1 &&
1204
                   (k1 == 0 || k2 == 0))))) {
1205
                if (tail == NULL) {
1206
                    *ret = t;
1207
                } else {
1208
                    tail->next = t;
1209
                }
1210
                tail = t;
1211
            } else {
1212
                free(t);
1213
            }
1214
        }
1215
        return 0;
1216
    }
1217
    *ret = NULL;
1218
    return 0;
1219
}
1220

1221
// One more layer of processing on top of the regular ellipse.
1222
// Uses the clipping tree.
1223
// Used for producing ellipse derivatives such as arc, chord, pie, etc.
1224
typedef struct {
1225
    ellipse_state st;
1226
    clip_node *root;
1227
    clip_node nodes[7];
1228
    int32_t node_count;
1229
    event_list *head;
1230
    int32_t y;
1231
} clip_ellipse_state;
1232

1233
typedef void (*clip_ellipse_init)(
1234
    clip_ellipse_state *, int32_t, int32_t, int32_t, float, float
1235
);
1236

1237
void
1238
debug_clip_tree(clip_node *root, int space) {
1239
    if (root == NULL) {
1240
        return;
1241
    }
1242
    if (root->type == CT_CLIP) {
1243
        int t = space;
1244
        while (t--) {
1245
            fputc(' ', stderr);
1246
        }
1247
        fprintf(stderr, "clip %+fx%+fy%+f > 0\n", root->a, root->b, root->c);
1248
    } else {
1249
        debug_clip_tree(root->l, space + 2);
1250
        int t = space;
1251
        while (t--) {
1252
            fputc(' ', stderr);
1253
        }
1254
        fprintf(stderr, "%s\n", root->type == CT_AND ? "and" : "or");
1255
        debug_clip_tree(root->r, space + 2);
1256
    }
1257
    if (space == 0) {
1258
        fputc('\n', stderr);
1259
    }
1260
}
1261

1262
// Resulting angles will satisfy 0 <= al < 360, al <= ar <= al + 360
1263
void
1264
normalize_angles(float *al, float *ar) {
1265
    if (*ar - *al >= 360) {
1266
        *al = 0;
1267
        *ar = 360;
1268
    } else {
1269
        *al = fmod(*al < 0 ? 360 - (fmod(-*al, 360)) : *al, 360);
1270
        *ar = *al + fmod(*ar < *al ? 360 - fmod(*al - *ar, 360) : *ar - *al, 360);
1271
    }
1272
}
1273

1274
// An arc with caps orthogonal to the ellipse curve.
1275
void
1276
arc_init(clip_ellipse_state *s, int32_t a, int32_t b, int32_t w, float al, float ar) {
1277
    if (a < b) {
1278
        // transpose the coordinate system
1279
        arc_init(s, b, a, w, 90 - ar, 90 - al);
1280
        ellipse_init(&s->st, a, b, w);
1281
        clip_tree_transpose(s->root);
1282
    } else {
1283
        // a >= b, based on "wide" ellipse
1284
        ellipse_init(&s->st, a, b, w);
1285

1286
        s->head = NULL;
1287
        s->node_count = 0;
1288
        normalize_angles(&al, &ar);
1289

1290
        // building clipping tree, a lot of different cases
1291
        if (ar == al + 360) {
1292
            s->root = NULL;
1293
        } else {
1294
            clip_node *lc = s->nodes + s->node_count++;
1295
            clip_node *rc = s->nodes + s->node_count++;
1296
            lc->l = lc->r = rc->l = rc->r = NULL;
1297
            lc->type = rc->type = CT_CLIP;
1298
            lc->a = -a * sin(al * M_PI / 180.0);
1299
            lc->b = b * cos(al * M_PI / 180.0);
1300
            lc->c = (a * a - b * b) * sin(al * M_PI / 90.0) / 2.0;
1301
            rc->a = a * sin(ar * M_PI / 180.0);
1302
            rc->b = -b * cos(ar * M_PI / 180.0);
1303
            rc->c = (b * b - a * a) * sin(ar * M_PI / 90.0) / 2.0;
1304
            if (fmod(al, 180) == 0 || fmod(ar, 180) == 0) {
1305
                s->root = s->nodes + s->node_count++;
1306
                s->root->l = lc;
1307
                s->root->r = rc;
1308
                s->root->type = ar - al < 180 ? CT_AND : CT_OR;
1309
            } else if (((int)(al / 180) + (int)(ar / 180)) % 2 == 1) {
1310
                s->root = s->nodes + s->node_count++;
1311
                s->root->l = s->nodes + s->node_count++;
1312
                s->root->l->l = s->nodes + s->node_count++;
1313
                s->root->l->r = lc;
1314
                s->root->r = s->nodes + s->node_count++;
1315
                s->root->r->l = s->nodes + s->node_count++;
1316
                s->root->r->r = rc;
1317
                s->root->type = CT_OR;
1318
                s->root->l->type = CT_AND;
1319
                s->root->r->type = CT_AND;
1320
                s->root->l->l->type = CT_CLIP;
1321
                s->root->r->l->type = CT_CLIP;
1322
                s->root->l->l->l = s->root->l->l->r = NULL;
1323
                s->root->r->l->l = s->root->r->l->r = NULL;
1324
                s->root->l->l->a = s->root->l->l->c = 0;
1325
                s->root->r->l->a = s->root->r->l->c = 0;
1326
                s->root->l->l->b = (int)(al / 180) % 2 == 0 ? 1 : -1;
1327
                s->root->r->l->b = (int)(ar / 180) % 2 == 0 ? 1 : -1;
1328
            } else {
1329
                s->root = s->nodes + s->node_count++;
1330
                s->root->l = s->nodes + s->node_count++;
1331
                s->root->r = s->nodes + s->node_count++;
1332
                s->root->type = s->root->l->type = ar - al < 180 ? CT_AND : CT_OR;
1333
                s->root->l->l = lc;
1334
                s->root->l->r = rc;
1335
                s->root->r->type = CT_CLIP;
1336
                s->root->r->l = s->root->r->r = NULL;
1337
                s->root->r->a = s->root->r->c = 0;
1338
                s->root->r->b = ar < 180 || ar > 540 ? 1 : -1;
1339
            }
1340
        }
1341
    }
1342
}
1343

1344
// A chord line.
1345
void
1346
chord_line_init(
1347
    clip_ellipse_state *s, int32_t a, int32_t b, int32_t w, float al, float ar
1348
) {
1349
    ellipse_init(&s->st, a, b, a + b + 1);
1350

1351
    s->head = NULL;
1352
    s->node_count = 0;
1353

1354
    // line equation for chord
1355
    double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
1356
    double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
1357
    s->root = s->nodes + s->node_count++;
1358
    s->root->l = s->nodes + s->node_count++;
1359
    s->root->r = s->nodes + s->node_count++;
1360
    s->root->type = CT_AND;
1361
    s->root->l->type = s->root->r->type = CT_CLIP;
1362
    s->root->l->l = s->root->l->r = s->root->r->l = s->root->r->r = NULL;
1363
    s->root->l->a = yr - yl;
1364
    s->root->l->b = xl - xr;
1365
    s->root->l->c = -(s->root->l->a * xl + s->root->l->b * yl);
1366
    s->root->r->a = -s->root->l->a;
1367
    s->root->r->b = -s->root->l->b;
1368
    s->root->r->c =
1369
        2 * w * sqrt(pow(s->root->l->a, 2.0) + pow(s->root->l->b, 2.0)) - s->root->l->c;
1370
}
1371

1372
// Pie side.
1373
void
1374
pie_side_init(
1375
    clip_ellipse_state *s, int32_t a, int32_t b, int32_t w, float al, float _
1376
) {
1377
    ellipse_init(&s->st, a, b, a + b + 1);
1378

1379
    s->head = NULL;
1380
    s->node_count = 0;
1381

1382
    double xl = a * cos(al * M_PI / 180.0);
1383
    double yl = b * sin(al * M_PI / 180.0);
1384
    double a1 = -yl;
1385
    double b1 = xl;
1386
    double c1 = w * sqrt(a1 * a1 + b1 * b1);
1387

1388
    s->root = s->nodes + s->node_count++;
1389
    s->root->type = CT_AND;
1390
    s->root->l = s->nodes + s->node_count++;
1391
    s->root->l->type = CT_AND;
1392

1393
    clip_node *cnode;
1394
    cnode = s->nodes + s->node_count++;
1395
    cnode->l = cnode->r = NULL;
1396
    cnode->type = CT_CLIP;
1397
    cnode->a = a1;
1398
    cnode->b = b1;
1399
    cnode->c = c1;
1400
    s->root->l->l = cnode;
1401
    cnode = s->nodes + s->node_count++;
1402
    cnode->l = cnode->r = NULL;
1403
    cnode->type = CT_CLIP;
1404
    cnode->a = -a1;
1405
    cnode->b = -b1;
1406
    cnode->c = c1;
1407
    s->root->l->r = cnode;
1408
    cnode = s->nodes + s->node_count++;
1409
    cnode->l = cnode->r = NULL;
1410
    cnode->type = CT_CLIP;
1411
    cnode->a = b1;
1412
    cnode->b = -a1;
1413
    cnode->c = 0;
1414
    s->root->r = cnode;
1415
}
1416

1417
// A chord.
1418
void
1419
chord_init(clip_ellipse_state *s, int32_t a, int32_t b, int32_t w, float al, float ar) {
1420
    ellipse_init(&s->st, a, b, w);
1421

1422
    s->head = NULL;
1423
    s->node_count = 0;
1424

1425
    // line equation for chord
1426
    double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
1427
    double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
1428
    s->root = s->nodes + s->node_count++;
1429
    s->root->l = s->root->r = NULL;
1430
    s->root->type = CT_CLIP;
1431
    s->root->a = yr - yl;
1432
    s->root->b = xl - xr;
1433
    s->root->c = -(s->root->a * xl + s->root->b * yl);
1434
}
1435

1436
// A pie. Can also be used to draw an arc with ugly sharp caps.
1437
void
1438
pie_init(clip_ellipse_state *s, int32_t a, int32_t b, int32_t w, float al, float ar) {
1439
    ellipse_init(&s->st, a, b, w);
1440

1441
    s->head = NULL;
1442
    s->node_count = 0;
1443

1444
    // line equations for pie sides
1445
    double xl = a * cos(al * M_PI / 180.0), xr = a * cos(ar * M_PI / 180.0);
1446
    double yl = b * sin(al * M_PI / 180.0), yr = b * sin(ar * M_PI / 180.0);
1447

1448
    clip_node *lc = s->nodes + s->node_count++;
1449
    clip_node *rc = s->nodes + s->node_count++;
1450
    lc->l = lc->r = rc->l = rc->r = NULL;
1451
    lc->type = rc->type = CT_CLIP;
1452
    lc->a = -yl;
1453
    lc->b = xl;
1454
    lc->c = 0;
1455
    rc->a = yr;
1456
    rc->b = -xr;
1457
    rc->c = 0;
1458

1459
    s->root = s->nodes + s->node_count++;
1460
    s->root->l = lc;
1461
    s->root->r = rc;
1462
    s->root->type = ar - al < 180 ? CT_AND : CT_OR;
1463

1464
    // add one more semiplane to avoid spikes
1465
    if (ar - al < 90) {
1466
        clip_node *old_root = s->root;
1467
        clip_node *spike_clipper = s->nodes + s->node_count++;
1468
        s->root = s->nodes + s->node_count++;
1469
        s->root->l = old_root;
1470
        s->root->r = spike_clipper;
1471
        s->root->type = CT_AND;
1472

1473
        spike_clipper->l = spike_clipper->r = NULL;
1474
        spike_clipper->type = CT_CLIP;
1475
        spike_clipper->a = (xl + xr) / 2.0;
1476
        spike_clipper->b = (yl + yr) / 2.0;
1477
        spike_clipper->c = 0;
1478
    }
1479
}
1480

1481
void
1482
clip_ellipse_free(clip_ellipse_state *s) {
1483
    while (s->head != NULL) {
1484
        event_list *t = s->head;
1485
        s->head = s->head->next;
1486
        free(t);
1487
    }
1488
}
1489

1490
int8_t
1491
clip_ellipse_next(
1492
    clip_ellipse_state *s, int32_t *ret_x0, int32_t *ret_y, int32_t *ret_x1
1493
) {
1494
    int32_t x0, y, x1;
1495
    while (s->head == NULL && ellipse_next(&s->st, &x0, &y, &x1) >= 0) {
1496
        if (clip_tree_do_clip(s->root, x0, y, x1, &s->head) < 0) {
1497
            return -2;
1498
        }
1499
        s->y = y;
1500
    }
1501
    if (s->head != NULL) {
1502
        *ret_y = s->y;
1503
        event_list *t = s->head;
1504
        s->head = s->head->next;
1505
        *ret_x0 = t->x;
1506
        free(t);
1507
        t = s->head;
1508
        assert(t != NULL);
1509
        s->head = s->head->next;
1510
        *ret_x1 = t->x;
1511
        free(t);
1512
        return 0;
1513
    }
1514
    return -1;
1515
}
1516

1517
static int
1518
ellipseNew(
1519
    Imaging im,
1520
    int x0,
1521
    int y0,
1522
    int x1,
1523
    int y1,
1524
    const void *ink_,
1525
    int fill,
1526
    int width,
1527
    int op
1528
) {
1529
    DRAW *draw;
1530
    INT32 ink;
1531
    DRAWINIT();
1532

1533
    int a = x1 - x0;
1534
    int b = y1 - y0;
1535
    if (a < 0 || b < 0) {
1536
        return 0;
1537
    }
1538
    if (fill) {
1539
        width = a + b;
1540
    }
1541

1542
    ellipse_state st;
1543
    ellipse_init(&st, a, b, width);
1544
    int32_t X0, Y, X1;
1545
    while (ellipse_next(&st, &X0, &Y, &X1) != -1) {
1546
        draw->hline(im, x0 + (X0 + a) / 2, y0 + (Y + b) / 2, x0 + (X1 + a) / 2, ink);
1547
    }
1548
    return 0;
1549
}
1550

1551
static int
1552
clipEllipseNew(
1553
    Imaging im,
1554
    int x0,
1555
    int y0,
1556
    int x1,
1557
    int y1,
1558
    float start,
1559
    float end,
1560
    const void *ink_,
1561
    int width,
1562
    int op,
1563
    clip_ellipse_init init
1564
) {
1565
    DRAW *draw;
1566
    INT32 ink;
1567
    DRAWINIT();
1568

1569
    int a = x1 - x0;
1570
    int b = y1 - y0;
1571
    if (a < 0 || b < 0) {
1572
        return 0;
1573
    }
1574

1575
    clip_ellipse_state st;
1576
    init(&st, a, b, width, start, end);
1577
    // debug_clip_tree(st.root, 0);
1578
    int32_t X0, Y, X1;
1579
    int next_code;
1580
    while ((next_code = clip_ellipse_next(&st, &X0, &Y, &X1)) >= 0) {
1581
        draw->hline(im, x0 + (X0 + a) / 2, y0 + (Y + b) / 2, x0 + (X1 + a) / 2, ink);
1582
    }
1583
    clip_ellipse_free(&st);
1584
    return next_code == -1 ? 0 : -1;
1585
}
1586
static int
1587
arcNew(
1588
    Imaging im,
1589
    int x0,
1590
    int y0,
1591
    int x1,
1592
    int y1,
1593
    float start,
1594
    float end,
1595
    const void *ink_,
1596
    int width,
1597
    int op
1598
) {
1599
    return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, arc_init);
1600
}
1601

1602
static int
1603
chordNew(
1604
    Imaging im,
1605
    int x0,
1606
    int y0,
1607
    int x1,
1608
    int y1,
1609
    float start,
1610
    float end,
1611
    const void *ink_,
1612
    int width,
1613
    int op
1614
) {
1615
    return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, chord_init);
1616
}
1617

1618
static int
1619
chordLineNew(
1620
    Imaging im,
1621
    int x0,
1622
    int y0,
1623
    int x1,
1624
    int y1,
1625
    float start,
1626
    float end,
1627
    const void *ink_,
1628
    int width,
1629
    int op
1630
) {
1631
    return clipEllipseNew(
1632
        im, x0, y0, x1, y1, start, end, ink_, width, op, chord_line_init
1633
    );
1634
}
1635

1636
static int
1637
pieNew(
1638
    Imaging im,
1639
    int x0,
1640
    int y0,
1641
    int x1,
1642
    int y1,
1643
    float start,
1644
    float end,
1645
    const void *ink_,
1646
    int width,
1647
    int op
1648
) {
1649
    return clipEllipseNew(im, x0, y0, x1, y1, start, end, ink_, width, op, pie_init);
1650
}
1651

1652
static int
1653
pieSideNew(
1654
    Imaging im,
1655
    int x0,
1656
    int y0,
1657
    int x1,
1658
    int y1,
1659
    float start,
1660
    const void *ink_,
1661
    int width,
1662
    int op
1663
) {
1664
    return clipEllipseNew(im, x0, y0, x1, y1, start, 0, ink_, width, op, pie_side_init);
1665
}
1666

1667
int
1668
ImagingDrawEllipse(
1669
    Imaging im,
1670
    int x0,
1671
    int y0,
1672
    int x1,
1673
    int y1,
1674
    const void *ink,
1675
    int fill,
1676
    int width,
1677
    int op
1678
) {
1679
    return ellipseNew(im, x0, y0, x1, y1, ink, fill, width, op);
1680
}
1681

1682
int
1683
ImagingDrawArc(
1684
    Imaging im,
1685
    int x0,
1686
    int y0,
1687
    int x1,
1688
    int y1,
1689
    float start,
1690
    float end,
1691
    const void *ink,
1692
    int width,
1693
    int op
1694
) {
1695
    normalize_angles(&start, &end);
1696
    if (start + 360 == end) {
1697
        return ImagingDrawEllipse(im, x0, y0, x1, y1, ink, 0, width, op);
1698
    }
1699
    if (start == end) {
1700
        return 0;
1701
    }
1702
    return arcNew(im, x0, y0, x1, y1, start, end, ink, width, op);
1703
}
1704

1705
int
1706
ImagingDrawChord(
1707
    Imaging im,
1708
    int x0,
1709
    int y0,
1710
    int x1,
1711
    int y1,
1712
    float start,
1713
    float end,
1714
    const void *ink,
1715
    int fill,
1716
    int width,
1717
    int op
1718
) {
1719
    normalize_angles(&start, &end);
1720
    if (start + 360 == end) {
1721
        return ImagingDrawEllipse(im, x0, y0, x1, y1, ink, fill, width, op);
1722
    }
1723
    if (start == end) {
1724
        return 0;
1725
    }
1726
    if (fill) {
1727
        return chordNew(im, x0, y0, x1, y1, start, end, ink, x1 - x0 + y1 - y0 + 1, op);
1728
    } else {
1729
        if (chordLineNew(im, x0, y0, x1, y1, start, end, ink, width, op)) {
1730
            return -1;
1731
        }
1732
        return chordNew(im, x0, y0, x1, y1, start, end, ink, width, op);
1733
    }
1734
}
1735

1736
int
1737
ImagingDrawPieslice(
1738
    Imaging im,
1739
    int x0,
1740
    int y0,
1741
    int x1,
1742
    int y1,
1743
    float start,
1744
    float end,
1745
    const void *ink,
1746
    int fill,
1747
    int width,
1748
    int op
1749
) {
1750
    normalize_angles(&start, &end);
1751
    if (start + 360 == end) {
1752
        return ellipseNew(im, x0, y0, x1, y1, ink, fill, width, op);
1753
    }
1754
    if (start == end) {
1755
        return 0;
1756
    }
1757
    if (fill) {
1758
        return pieNew(im, x0, y0, x1, y1, start, end, ink, x1 + y1 - x0 - y0, op);
1759
    } else {
1760
        if (pieSideNew(im, x0, y0, x1, y1, start, ink, width, op)) {
1761
            return -1;
1762
        }
1763
        if (pieSideNew(im, x0, y0, x1, y1, end, ink, width, op)) {
1764
            return -1;
1765
        }
1766
        int xc = lround((x0 + x1 - width) / 2.0), yc = lround((y0 + y1 - width) / 2.0);
1767
        ellipseNew(im, xc, yc, xc + width - 1, yc + width - 1, ink, 1, 0, op);
1768
        return pieNew(im, x0, y0, x1, y1, start, end, ink, width, op);
1769
    }
1770
}
1771

1772
/* -------------------------------------------------------------------- */
1773

1774
/* experimental level 2 ("arrow") graphics stuff.  this implements
1775
   portions of the arrow api on top of the Edge structure.  the
1776
   semantics are ok, except that "curve" flattens the bezier curves by
1777
   itself */
1778

1779
struct ImagingOutlineInstance {
1780
    float x0, y0;
1781

1782
    float x, y;
1783

1784
    int count;
1785
    Edge *edges;
1786

1787
    int size;
1788
};
1789

1790
ImagingOutline
1791
ImagingOutlineNew(void) {
1792
    ImagingOutline outline;
1793

1794
    outline = calloc(1, sizeof(struct ImagingOutlineInstance));
1795
    if (!outline) {
1796
        return (ImagingOutline)ImagingError_MemoryError();
1797
    }
1798

1799
    outline->edges = NULL;
1800
    outline->count = outline->size = 0;
1801

1802
    ImagingOutlineMove(outline, 0, 0);
1803

1804
    return outline;
1805
}
1806

1807
void
1808
ImagingOutlineDelete(ImagingOutline outline) {
1809
    if (!outline) {
1810
        return;
1811
    }
1812

1813
    if (outline->edges) {
1814
        free(outline->edges);
1815
    }
1816

1817
    free(outline);
1818
}
1819

1820
static Edge *
1821
allocate(ImagingOutline outline, int extra) {
1822
    Edge *e;
1823

1824
    if (outline->count + extra > outline->size) {
1825
        /* expand outline buffer */
1826
        outline->size += extra + 25;
1827
        if (!outline->edges) {
1828
            /* malloc check ok, uses calloc for overflow */
1829
            e = calloc(outline->size, sizeof(Edge));
1830
        } else {
1831
            if (outline->size > INT_MAX / (int)sizeof(Edge)) {
1832
                return NULL;
1833
            }
1834
            /* malloc check ok, overflow checked above */
1835
            e = realloc(outline->edges, outline->size * sizeof(Edge));
1836
        }
1837
        if (!e) {
1838
            return NULL;
1839
        }
1840
        outline->edges = e;
1841
    }
1842

1843
    e = outline->edges + outline->count;
1844

1845
    outline->count += extra;
1846

1847
    return e;
1848
}
1849

1850
int
1851
ImagingOutlineMove(ImagingOutline outline, float x0, float y0) {
1852
    outline->x = outline->x0 = x0;
1853
    outline->y = outline->y0 = y0;
1854

1855
    return 0;
1856
}
1857

1858
int
1859
ImagingOutlineLine(ImagingOutline outline, float x1, float y1) {
1860
    Edge *e;
1861

1862
    e = allocate(outline, 1);
1863
    if (!e) {
1864
        return -1; /* out of memory */
1865
    }
1866

1867
    add_edge(e, (int)outline->x, (int)outline->y, (int)x1, (int)y1);
1868

1869
    outline->x = x1;
1870
    outline->y = y1;
1871

1872
    return 0;
1873
}
1874

1875
int
1876
ImagingOutlineCurve(
1877
    ImagingOutline outline, float x1, float y1, float x2, float y2, float x3, float y3
1878
) {
1879
    Edge *e;
1880
    int i;
1881
    float xo, yo;
1882

1883
#define STEPS 32
1884

1885
    e = allocate(outline, STEPS);
1886
    if (!e) {
1887
        return -1; /* out of memory */
1888
    }
1889

1890
    xo = outline->x;
1891
    yo = outline->y;
1892

1893
    /* flatten the bezier segment */
1894

1895
    for (i = 1; i <= STEPS; i++) {
1896
        float t = ((float)i) / STEPS;
1897
        float t2 = t * t;
1898
        float t3 = t2 * t;
1899

1900
        float u = 1.0F - t;
1901
        float u2 = u * u;
1902
        float u3 = u2 * u;
1903

1904
        float x = outline->x * u3 + 3 * (x1 * t * u2 + x2 * t2 * u) + x3 * t3 + 0.5;
1905
        float y = outline->y * u3 + 3 * (y1 * t * u2 + y2 * t2 * u) + y3 * t3 + 0.5;
1906

1907
        add_edge(e++, xo, yo, (int)x, (int)y);
1908

1909
        xo = x, yo = y;
1910
    }
1911

1912
    outline->x = xo;
1913
    outline->y = yo;
1914

1915
    return 0;
1916
}
1917

1918
int
1919
ImagingOutlineClose(ImagingOutline outline) {
1920
    if (outline->x == outline->x0 && outline->y == outline->y0) {
1921
        return 0;
1922
    }
1923
    return ImagingOutlineLine(outline, outline->x0, outline->y0);
1924
}
1925

1926
int
1927
ImagingOutlineTransform(ImagingOutline outline, double a[6]) {
1928
    Edge *eIn;
1929
    Edge *eOut;
1930
    int i, n;
1931
    int x0, y0, x1, y1;
1932
    int X0, Y0, X1, Y1;
1933

1934
    double a0 = a[0];
1935
    double a1 = a[1];
1936
    double a2 = a[2];
1937
    double a3 = a[3];
1938
    double a4 = a[4];
1939
    double a5 = a[5];
1940

1941
    eIn = outline->edges;
1942
    n = outline->count;
1943

1944
    eOut = allocate(outline, n);
1945
    if (!eOut) {
1946
        ImagingError_MemoryError();
1947
        return -1;
1948
    }
1949

1950
    for (i = 0; i < n; i++) {
1951
        x0 = eIn->x0;
1952
        y0 = eIn->y0;
1953

1954
        /* FIXME: ouch! */
1955
        if (eIn->x0 == eIn->xmin) {
1956
            x1 = eIn->xmax;
1957
        } else {
1958
            x1 = eIn->xmin;
1959
        }
1960
        if (eIn->y0 == eIn->ymin) {
1961
            y1 = eIn->ymax;
1962
        } else {
1963
            y1 = eIn->ymin;
1964
        }
1965

1966
        /* full moon tonight!  if this doesn't work, you may need to
1967
           upgrade your compiler (make sure you have the right service
1968
           pack) */
1969

1970
        X0 = (int)(a0 * x0 + a1 * y0 + a2);
1971
        Y0 = (int)(a3 * x0 + a4 * y0 + a5);
1972
        X1 = (int)(a0 * x1 + a1 * y1 + a2);
1973
        Y1 = (int)(a3 * x1 + a4 * y1 + a5);
1974

1975
        add_edge(eOut, X0, Y0, X1, Y1);
1976

1977
        eIn++;
1978
        eOut++;
1979
    }
1980

1981
    free(outline->edges);
1982

1983
    /* FIXME: ugly! */
1984
    outline->edges = NULL;
1985
    outline->count = outline->size = 0;
1986

1987
    return 0;
1988
}
1989

1990
int
1991
ImagingDrawOutline(
1992
    Imaging im, ImagingOutline outline, const void *ink_, int fill, int op
1993
) {
1994
    DRAW *draw;
1995
    INT32 ink;
1996

1997
    DRAWINIT();
1998

1999
    draw->polygon(im, outline->count, outline->edges, ink, 0);
2000

2001
    return 0;
2002
}
2003