Pillow

1
#include "Imaging.h"
2

3
/* For large images rotation is an inefficient operation in terms of CPU cache.
4
   One row in the source image affects each column in destination.
5
   Rotating in chunks that fit in the cache can speed up rotation
6
   8x on a modern CPU. A chunk size of 128 requires only 65k and is large enough
7
   that the overhead from the extra loops are not apparent. */
8
#define ROTATE_CHUNK 512
9
#define ROTATE_SMALL_CHUNK 8
10

11
#define COORD(v) ((v) < 0.0 ? -1 : ((int)(v)))
12
#define FLOOR(v) ((v) < 0.0 ? ((int)floor(v)) : ((int)(v)))
13

14
/* -------------------------------------------------------------------- */
15
/* Transpose operations                                                 */
16

17
Imaging
18
ImagingFlipLeftRight(Imaging imOut, Imaging imIn) {
19
    ImagingSectionCookie cookie;
20
    int x, y, xr;
21

22
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
23
        return (Imaging)ImagingError_ModeError();
24
    }
25
    if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize) {
26
        return (Imaging)ImagingError_Mismatch();
27
    }
28

29
    ImagingCopyPalette(imOut, imIn);
30

31
#define FLIP_LEFT_RIGHT(INT, image)               \
32
    for (y = 0; y < imIn->ysize; y++) {           \
33
        INT *in = (INT *)imIn->image[y];          \
34
        INT *out = (INT *)imOut->image[y];        \
35
        xr = imIn->xsize - 1;                     \
36
        for (x = 0; x < imIn->xsize; x++, xr--) { \
37
            out[xr] = in[x];                      \
38
        }                                         \
39
    }
40

41
    ImagingSectionEnter(&cookie);
42

43
    if (imIn->image8) {
44
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
45
            FLIP_LEFT_RIGHT(UINT16, image8)
46
        } else {
47
            FLIP_LEFT_RIGHT(UINT8, image8)
48
        }
49
    } else {
50
        FLIP_LEFT_RIGHT(INT32, image32)
51
    }
52

53
    ImagingSectionLeave(&cookie);
54

55
#undef FLIP_LEFT_RIGHT
56

57
    return imOut;
58
}
59

60
Imaging
61
ImagingFlipTopBottom(Imaging imOut, Imaging imIn) {
62
    ImagingSectionCookie cookie;
63
    int y, yr;
64

65
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
66
        return (Imaging)ImagingError_ModeError();
67
    }
68
    if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize) {
69
        return (Imaging)ImagingError_Mismatch();
70
    }
71

72
    ImagingCopyPalette(imOut, imIn);
73

74
    ImagingSectionEnter(&cookie);
75

76
    yr = imIn->ysize - 1;
77
    for (y = 0; y < imIn->ysize; y++, yr--) {
78
        memcpy(imOut->image[yr], imIn->image[y], imIn->linesize);
79
    }
80

81
    ImagingSectionLeave(&cookie);
82

83
    return imOut;
84
}
85

86
Imaging
87
ImagingRotate90(Imaging imOut, Imaging imIn) {
88
    ImagingSectionCookie cookie;
89
    int x, y, xx, yy, xr, xxsize, yysize;
90
    int xxx, yyy, xxxsize, yyysize;
91

92
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
93
        return (Imaging)ImagingError_ModeError();
94
    }
95
    if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize) {
96
        return (Imaging)ImagingError_Mismatch();
97
    }
98

99
    ImagingCopyPalette(imOut, imIn);
100

101
#define ROTATE_90(INT, image)                                                         \
102
    for (y = 0; y < imIn->ysize; y += ROTATE_CHUNK) {                                 \
103
        for (x = 0; x < imIn->xsize; x += ROTATE_CHUNK) {                             \
104
            yysize = y + ROTATE_CHUNK < imIn->ysize ? y + ROTATE_CHUNK : imIn->ysize; \
105
            xxsize = x + ROTATE_CHUNK < imIn->xsize ? x + ROTATE_CHUNK : imIn->xsize; \
106
            for (yy = y; yy < yysize; yy += ROTATE_SMALL_CHUNK) {                     \
107
                for (xx = x; xx < xxsize; xx += ROTATE_SMALL_CHUNK) {                 \
108
                    yyysize = yy + ROTATE_SMALL_CHUNK < imIn->ysize                   \
109
                                  ? yy + ROTATE_SMALL_CHUNK                           \
110
                                  : imIn->ysize;                                      \
111
                    xxxsize = xx + ROTATE_SMALL_CHUNK < imIn->xsize                   \
112
                                  ? xx + ROTATE_SMALL_CHUNK                           \
113
                                  : imIn->xsize;                                      \
114
                    for (yyy = yy; yyy < yyysize; yyy++) {                            \
115
                        INT *in = (INT *)imIn->image[yyy];                            \
116
                        xr = imIn->xsize - 1 - xx;                                    \
117
                        for (xxx = xx; xxx < xxxsize; xxx++, xr--) {                  \
118
                            INT *out = (INT *)imOut->image[xr];                       \
119
                            out[yyy] = in[xxx];                                       \
120
                        }                                                             \
121
                    }                                                                 \
122
                }                                                                     \
123
            }                                                                         \
124
        }                                                                             \
125
    }
126

127
    ImagingSectionEnter(&cookie);
128

129
    if (imIn->image8) {
130
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
131
            ROTATE_90(UINT16, image8);
132
        } else {
133
            ROTATE_90(UINT8, image8);
134
        }
135
    } else {
136
        ROTATE_90(INT32, image32);
137
    }
138

139
    ImagingSectionLeave(&cookie);
140

141
#undef ROTATE_90
142

143
    return imOut;
144
}
145

146
Imaging
147
ImagingTranspose(Imaging imOut, Imaging imIn) {
148
    ImagingSectionCookie cookie;
149
    int x, y, xx, yy, xxsize, yysize;
150
    int xxx, yyy, xxxsize, yyysize;
151

152
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
153
        return (Imaging)ImagingError_ModeError();
154
    }
155
    if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize) {
156
        return (Imaging)ImagingError_Mismatch();
157
    }
158

159
    ImagingCopyPalette(imOut, imIn);
160

161
#define TRANSPOSE(INT, image)                                                         \
162
    for (y = 0; y < imIn->ysize; y += ROTATE_CHUNK) {                                 \
163
        for (x = 0; x < imIn->xsize; x += ROTATE_CHUNK) {                             \
164
            yysize = y + ROTATE_CHUNK < imIn->ysize ? y + ROTATE_CHUNK : imIn->ysize; \
165
            xxsize = x + ROTATE_CHUNK < imIn->xsize ? x + ROTATE_CHUNK : imIn->xsize; \
166
            for (yy = y; yy < yysize; yy += ROTATE_SMALL_CHUNK) {                     \
167
                for (xx = x; xx < xxsize; xx += ROTATE_SMALL_CHUNK) {                 \
168
                    yyysize = yy + ROTATE_SMALL_CHUNK < imIn->ysize                   \
169
                                  ? yy + ROTATE_SMALL_CHUNK                           \
170
                                  : imIn->ysize;                                      \
171
                    xxxsize = xx + ROTATE_SMALL_CHUNK < imIn->xsize                   \
172
                                  ? xx + ROTATE_SMALL_CHUNK                           \
173
                                  : imIn->xsize;                                      \
174
                    for (yyy = yy; yyy < yyysize; yyy++) {                            \
175
                        INT *in = (INT *)imIn->image[yyy];                            \
176
                        for (xxx = xx; xxx < xxxsize; xxx++) {                        \
177
                            INT *out = (INT *)imOut->image[xxx];                      \
178
                            out[yyy] = in[xxx];                                       \
179
                        }                                                             \
180
                    }                                                                 \
181
                }                                                                     \
182
            }                                                                         \
183
        }                                                                             \
184
    }
185

186
    ImagingSectionEnter(&cookie);
187

188
    if (imIn->image8) {
189
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
190
            TRANSPOSE(UINT16, image8);
191
        } else {
192
            TRANSPOSE(UINT8, image8);
193
        }
194
    } else {
195
        TRANSPOSE(INT32, image32);
196
    }
197

198
    ImagingSectionLeave(&cookie);
199

200
#undef TRANSPOSE
201

202
    return imOut;
203
}
204

205
Imaging
206
ImagingTransverse(Imaging imOut, Imaging imIn) {
207
    ImagingSectionCookie cookie;
208
    int x, y, xr, yr, xx, yy, xxsize, yysize;
209
    int xxx, yyy, xxxsize, yyysize;
210

211
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
212
        return (Imaging)ImagingError_ModeError();
213
    }
214
    if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize) {
215
        return (Imaging)ImagingError_Mismatch();
216
    }
217

218
    ImagingCopyPalette(imOut, imIn);
219

220
#define TRANSVERSE(INT, image)                                                        \
221
    for (y = 0; y < imIn->ysize; y += ROTATE_CHUNK) {                                 \
222
        for (x = 0; x < imIn->xsize; x += ROTATE_CHUNK) {                             \
223
            yysize = y + ROTATE_CHUNK < imIn->ysize ? y + ROTATE_CHUNK : imIn->ysize; \
224
            xxsize = x + ROTATE_CHUNK < imIn->xsize ? x + ROTATE_CHUNK : imIn->xsize; \
225
            for (yy = y; yy < yysize; yy += ROTATE_SMALL_CHUNK) {                     \
226
                for (xx = x; xx < xxsize; xx += ROTATE_SMALL_CHUNK) {                 \
227
                    yyysize = yy + ROTATE_SMALL_CHUNK < imIn->ysize                   \
228
                                  ? yy + ROTATE_SMALL_CHUNK                           \
229
                                  : imIn->ysize;                                      \
230
                    xxxsize = xx + ROTATE_SMALL_CHUNK < imIn->xsize                   \
231
                                  ? xx + ROTATE_SMALL_CHUNK                           \
232
                                  : imIn->xsize;                                      \
233
                    yr = imIn->ysize - 1 - yy;                                        \
234
                    for (yyy = yy; yyy < yyysize; yyy++, yr--) {                      \
235
                        INT *in = (INT *)imIn->image[yyy];                            \
236
                        xr = imIn->xsize - 1 - xx;                                    \
237
                        for (xxx = xx; xxx < xxxsize; xxx++, xr--) {                  \
238
                            INT *out = (INT *)imOut->image[xr];                       \
239
                            out[yr] = in[xxx];                                        \
240
                        }                                                             \
241
                    }                                                                 \
242
                }                                                                     \
243
            }                                                                         \
244
        }                                                                             \
245
    }
246

247
    ImagingSectionEnter(&cookie);
248

249
    if (imIn->image8) {
250
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
251
            TRANSVERSE(UINT16, image8);
252
        } else {
253
            TRANSVERSE(UINT8, image8);
254
        }
255
    } else {
256
        TRANSVERSE(INT32, image32);
257
    }
258

259
    ImagingSectionLeave(&cookie);
260

261
#undef TRANSVERSE
262

263
    return imOut;
264
}
265

266
Imaging
267
ImagingRotate180(Imaging imOut, Imaging imIn) {
268
    ImagingSectionCookie cookie;
269
    int x, y, xr, yr;
270

271
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
272
        return (Imaging)ImagingError_ModeError();
273
    }
274
    if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize) {
275
        return (Imaging)ImagingError_Mismatch();
276
    }
277

278
    ImagingCopyPalette(imOut, imIn);
279

280
#define ROTATE_180(INT, image)                    \
281
    for (y = 0; y < imIn->ysize; y++, yr--) {     \
282
        INT *in = (INT *)imIn->image[y];          \
283
        INT *out = (INT *)imOut->image[yr];       \
284
        xr = imIn->xsize - 1;                     \
285
        for (x = 0; x < imIn->xsize; x++, xr--) { \
286
            out[xr] = in[x];                      \
287
        }                                         \
288
    }
289

290
    ImagingSectionEnter(&cookie);
291

292
    yr = imIn->ysize - 1;
293
    if (imIn->image8) {
294
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
295
            ROTATE_180(UINT16, image8)
296
        } else {
297
            ROTATE_180(UINT8, image8)
298
        }
299
    } else {
300
        ROTATE_180(INT32, image32)
301
    }
302

303
    ImagingSectionLeave(&cookie);
304

305
#undef ROTATE_180
306

307
    return imOut;
308
}
309

310
Imaging
311
ImagingRotate270(Imaging imOut, Imaging imIn) {
312
    ImagingSectionCookie cookie;
313
    int x, y, xx, yy, yr, xxsize, yysize;
314
    int xxx, yyy, xxxsize, yyysize;
315

316
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
317
        return (Imaging)ImagingError_ModeError();
318
    }
319
    if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize) {
320
        return (Imaging)ImagingError_Mismatch();
321
    }
322

323
    ImagingCopyPalette(imOut, imIn);
324

325
#define ROTATE_270(INT, image)                                                        \
326
    for (y = 0; y < imIn->ysize; y += ROTATE_CHUNK) {                                 \
327
        for (x = 0; x < imIn->xsize; x += ROTATE_CHUNK) {                             \
328
            yysize = y + ROTATE_CHUNK < imIn->ysize ? y + ROTATE_CHUNK : imIn->ysize; \
329
            xxsize = x + ROTATE_CHUNK < imIn->xsize ? x + ROTATE_CHUNK : imIn->xsize; \
330
            for (yy = y; yy < yysize; yy += ROTATE_SMALL_CHUNK) {                     \
331
                for (xx = x; xx < xxsize; xx += ROTATE_SMALL_CHUNK) {                 \
332
                    yyysize = yy + ROTATE_SMALL_CHUNK < imIn->ysize                   \
333
                                  ? yy + ROTATE_SMALL_CHUNK                           \
334
                                  : imIn->ysize;                                      \
335
                    xxxsize = xx + ROTATE_SMALL_CHUNK < imIn->xsize                   \
336
                                  ? xx + ROTATE_SMALL_CHUNK                           \
337
                                  : imIn->xsize;                                      \
338
                    yr = imIn->ysize - 1 - yy;                                        \
339
                    for (yyy = yy; yyy < yyysize; yyy++, yr--) {                      \
340
                        INT *in = (INT *)imIn->image[yyy];                            \
341
                        for (xxx = xx; xxx < xxxsize; xxx++) {                        \
342
                            INT *out = (INT *)imOut->image[xxx];                      \
343
                            out[yr] = in[xxx];                                        \
344
                        }                                                             \
345
                    }                                                                 \
346
                }                                                                     \
347
            }                                                                         \
348
        }                                                                             \
349
    }
350

351
    ImagingSectionEnter(&cookie);
352

353
    if (imIn->image8) {
354
        if (strncmp(imIn->mode, "I;16", 4) == 0) {
355
            ROTATE_270(UINT16, image8);
356
        } else {
357
            ROTATE_270(UINT8, image8);
358
        }
359
    } else {
360
        ROTATE_270(INT32, image32);
361
    }
362

363
    ImagingSectionLeave(&cookie);
364

365
#undef ROTATE_270
366

367
    return imOut;
368
}
369

370
/* -------------------------------------------------------------------- */
371
/* Transforms                                                           */
372

373
/* transform primitives (ImagingTransformMap) */
374

375
static int
376
affine_transform(double *xout, double *yout, int x, int y, void *data) {
377
    /* full moon tonight.  your compiler will generate bogus code
378
       for simple expressions, unless you reorganize the code, or
379
       install Service Pack 3 */
380

381
    double *a = (double *)data;
382
    double a0 = a[0];
383
    double a1 = a[1];
384
    double a2 = a[2];
385
    double a3 = a[3];
386
    double a4 = a[4];
387
    double a5 = a[5];
388

389
    double xin = x + 0.5;
390
    double yin = y + 0.5;
391

392
    xout[0] = a0 * xin + a1 * yin + a2;
393
    yout[0] = a3 * xin + a4 * yin + a5;
394

395
    return 1;
396
}
397

398
static int
399
perspective_transform(double *xout, double *yout, int x, int y, void *data) {
400
    double *a = (double *)data;
401
    double a0 = a[0];
402
    double a1 = a[1];
403
    double a2 = a[2];
404
    double a3 = a[3];
405
    double a4 = a[4];
406
    double a5 = a[5];
407
    double a6 = a[6];
408
    double a7 = a[7];
409

410
    double xin = x + 0.5;
411
    double yin = y + 0.5;
412

413
    xout[0] = (a0 * xin + a1 * yin + a2) / (a6 * xin + a7 * yin + 1);
414
    yout[0] = (a3 * xin + a4 * yin + a5) / (a6 * xin + a7 * yin + 1);
415

416
    return 1;
417
}
418

419
static int
420
quad_transform(double *xout, double *yout, int x, int y, void *data) {
421
    /* quad warp: map quadrilateral to rectangle */
422

423
    double *a = (double *)data;
424
    double a0 = a[0];
425
    double a1 = a[1];
426
    double a2 = a[2];
427
    double a3 = a[3];
428
    double a4 = a[4];
429
    double a5 = a[5];
430
    double a6 = a[6];
431
    double a7 = a[7];
432

433
    double xin = x + 0.5;
434
    double yin = y + 0.5;
435

436
    xout[0] = a0 + a1 * xin + a2 * yin + a3 * xin * yin;
437
    yout[0] = a4 + a5 * xin + a6 * yin + a7 * xin * yin;
438

439
    return 1;
440
}
441

442
/* transform filters (ImagingTransformFilter) */
443

444
static int
445
nearest_filter8(void *out, Imaging im, double xin, double yin) {
446
    int x = COORD(xin);
447
    int y = COORD(yin);
448
    if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize) {
449
        return 0;
450
    }
451
    ((UINT8 *)out)[0] = im->image8[y][x];
452
    return 1;
453
}
454

455
static int
456
nearest_filter16(void *out, Imaging im, double xin, double yin) {
457
    int x = COORD(xin);
458
    int y = COORD(yin);
459
    if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize) {
460
        return 0;
461
    }
462
    memcpy(out, im->image8[y] + x * sizeof(INT16), sizeof(INT16));
463
    return 1;
464
}
465

466
static int
467
nearest_filter32(void *out, Imaging im, double xin, double yin) {
468
    int x = COORD(xin);
469
    int y = COORD(yin);
470
    if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize) {
471
        return 0;
472
    }
473
    memcpy(out, &im->image32[y][x], sizeof(INT32));
474
    return 1;
475
}
476

477
#define XCLIP(im, x) (((x) < 0) ? 0 : ((x) < im->xsize) ? (x) : im->xsize - 1)
478
#define YCLIP(im, y) (((y) < 0) ? 0 : ((y) < im->ysize) ? (y) : im->ysize - 1)
479

480
#define BILINEAR(v, a, b, d) (v = (a) + ((b) - (a)) * (d))
481

482
#define BILINEAR_HEAD(type)                                               \
483
    int x, y;                                                             \
484
    int x0, x1;                                                           \
485
    double v1, v2;                                                        \
486
    double dx, dy;                                                        \
487
    type *in;                                                             \
488
    if (xin < 0.0 || xin >= im->xsize || yin < 0.0 || yin >= im->ysize) { \
489
        return 0;                                                         \
490
    }                                                                     \
491
    xin -= 0.5;                                                           \
492
    yin -= 0.5;                                                           \
493
    x = FLOOR(xin);                                                       \
494
    y = FLOOR(yin);                                                       \
495
    dx = xin - x;                                                         \
496
    dy = yin - y;
497

498
#define BILINEAR_BODY(type, image, step, offset)       \
499
    {                                                  \
500
        in = (type *)((image)[YCLIP(im, y)] + offset); \
501
        x0 = XCLIP(im, x + 0) * step;                  \
502
        x1 = XCLIP(im, x + 1) * step;                  \
503
        BILINEAR(v1, in[x0], in[x1], dx);              \
504
        if (y + 1 >= 0 && y + 1 < im->ysize) {         \
505
            in = (type *)((image)[y + 1] + offset);    \
506
            BILINEAR(v2, in[x0], in[x1], dx);          \
507
        } else {                                       \
508
            v2 = v1;                                   \
509
        }                                              \
510
        BILINEAR(v1, v1, v2, dy);                      \
511
    }
512

513
static int
514
bilinear_filter8(void *out, Imaging im, double xin, double yin) {
515
    BILINEAR_HEAD(UINT8);
516
    BILINEAR_BODY(UINT8, im->image8, 1, 0);
517
    ((UINT8 *)out)[0] = (UINT8)v1;
518
    return 1;
519
}
520

521
static int
522
bilinear_filter32I(void *out, Imaging im, double xin, double yin) {
523
    INT32 k;
524
    BILINEAR_HEAD(INT32);
525
    BILINEAR_BODY(INT32, im->image32, 1, 0);
526
    k = v1;
527
    memcpy(out, &k, sizeof(k));
528
    return 1;
529
}
530

531
static int
532
bilinear_filter32F(void *out, Imaging im, double xin, double yin) {
533
    FLOAT32 k;
534
    BILINEAR_HEAD(FLOAT32);
535
    BILINEAR_BODY(FLOAT32, im->image32, 1, 0);
536
    k = v1;
537
    memcpy(out, &k, sizeof(k));
538
    return 1;
539
}
540

541
static int
542
bilinear_filter32LA(void *out, Imaging im, double xin, double yin) {
543
    BILINEAR_HEAD(UINT8);
544
    BILINEAR_BODY(UINT8, im->image, 4, 0);
545
    ((UINT8 *)out)[0] = (UINT8)v1;
546
    ((UINT8 *)out)[1] = (UINT8)v1;
547
    ((UINT8 *)out)[2] = (UINT8)v1;
548
    BILINEAR_BODY(UINT8, im->image, 4, 3);
549
    ((UINT8 *)out)[3] = (UINT8)v1;
550
    return 1;
551
}
552

553
static int
554
bilinear_filter32RGB(void *out, Imaging im, double xin, double yin) {
555
    int b;
556
    BILINEAR_HEAD(UINT8);
557
    for (b = 0; b < im->bands; b++) {
558
        BILINEAR_BODY(UINT8, im->image, 4, b);
559
        ((UINT8 *)out)[b] = (UINT8)v1;
560
    }
561
    return 1;
562
}
563

564
#undef BILINEAR
565
#undef BILINEAR_HEAD
566
#undef BILINEAR_BODY
567

568
#define BICUBIC(v, v1, v2, v3, v4, d)                \
569
    {                                                \
570
        double p1 = v2;                              \
571
        double p2 = -v1 + v3;                        \
572
        double p3 = 2 * (v1 - v2) + v3 - v4;         \
573
        double p4 = -v1 + v2 - v3 + v4;              \
574
        v = p1 + (d) * (p2 + (d) * (p3 + (d) * p4)); \
575
    }
576

577
#define BICUBIC_HEAD(type)                                                \
578
    int x = FLOOR(xin);                                                   \
579
    int y = FLOOR(yin);                                                   \
580
    int x0, x1, x2, x3;                                                   \
581
    double v1, v2, v3, v4;                                                \
582
    double dx, dy;                                                        \
583
    type *in;                                                             \
584
    if (xin < 0.0 || xin >= im->xsize || yin < 0.0 || yin >= im->ysize) { \
585
        return 0;                                                         \
586
    }                                                                     \
587
    xin -= 0.5;                                                           \
588
    yin -= 0.5;                                                           \
589
    x = FLOOR(xin);                                                       \
590
    y = FLOOR(yin);                                                       \
591
    dx = xin - x;                                                         \
592
    dy = yin - y;                                                         \
593
    x--;                                                                  \
594
    y--;
595

596
#define BICUBIC_BODY(type, image, step, offset)              \
597
    {                                                        \
598
        in = (type *)((image)[YCLIP(im, y)] + offset);       \
599
        x0 = XCLIP(im, x + 0) * step;                        \
600
        x1 = XCLIP(im, x + 1) * step;                        \
601
        x2 = XCLIP(im, x + 2) * step;                        \
602
        x3 = XCLIP(im, x + 3) * step;                        \
603
        BICUBIC(v1, in[x0], in[x1], in[x2], in[x3], dx);     \
604
        if (y + 1 >= 0 && y + 1 < im->ysize) {               \
605
            in = (type *)((image)[y + 1] + offset);          \
606
            BICUBIC(v2, in[x0], in[x1], in[x2], in[x3], dx); \
607
        } else {                                             \
608
            v2 = v1;                                         \
609
        }                                                    \
610
        if (y + 2 >= 0 && y + 2 < im->ysize) {               \
611
            in = (type *)((image)[y + 2] + offset);          \
612
            BICUBIC(v3, in[x0], in[x1], in[x2], in[x3], dx); \
613
        } else {                                             \
614
            v3 = v2;                                         \
615
        }                                                    \
616
        if (y + 3 >= 0 && y + 3 < im->ysize) {               \
617
            in = (type *)((image)[y + 3] + offset);          \
618
            BICUBIC(v4, in[x0], in[x1], in[x2], in[x3], dx); \
619
        } else {                                             \
620
            v4 = v3;                                         \
621
        }                                                    \
622
        BICUBIC(v1, v1, v2, v3, v4, dy);                     \
623
    }
624

625
static int
626
bicubic_filter8(void *out, Imaging im, double xin, double yin) {
627
    BICUBIC_HEAD(UINT8);
628
    BICUBIC_BODY(UINT8, im->image8, 1, 0);
629
    if (v1 <= 0.0) {
630
        ((UINT8 *)out)[0] = 0;
631
    } else if (v1 >= 255.0) {
632
        ((UINT8 *)out)[0] = 255;
633
    } else {
634
        ((UINT8 *)out)[0] = (UINT8)v1;
635
    }
636
    return 1;
637
}
638

639
static int
640
bicubic_filter32I(void *out, Imaging im, double xin, double yin) {
641
    INT32 k;
642
    BICUBIC_HEAD(INT32);
643
    BICUBIC_BODY(INT32, im->image32, 1, 0);
644
    k = v1;
645
    memcpy(out, &k, sizeof(k));
646
    return 1;
647
}
648

649
static int
650
bicubic_filter32F(void *out, Imaging im, double xin, double yin) {
651
    FLOAT32 k;
652
    BICUBIC_HEAD(FLOAT32);
653
    BICUBIC_BODY(FLOAT32, im->image32, 1, 0);
654
    k = v1;
655
    memcpy(out, &k, sizeof(k));
656
    return 1;
657
}
658

659
static int
660
bicubic_filter32LA(void *out, Imaging im, double xin, double yin) {
661
    BICUBIC_HEAD(UINT8);
662
    BICUBIC_BODY(UINT8, im->image, 4, 0);
663
    if (v1 <= 0.0) {
664
        ((UINT8 *)out)[0] = 0;
665
        ((UINT8 *)out)[1] = 0;
666
        ((UINT8 *)out)[2] = 0;
667
    } else if (v1 >= 255.0) {
668
        ((UINT8 *)out)[0] = 255;
669
        ((UINT8 *)out)[1] = 255;
670
        ((UINT8 *)out)[2] = 255;
671
    } else {
672
        ((UINT8 *)out)[0] = (UINT8)v1;
673
        ((UINT8 *)out)[1] = (UINT8)v1;
674
        ((UINT8 *)out)[2] = (UINT8)v1;
675
    }
676
    BICUBIC_BODY(UINT8, im->image, 4, 3);
677
    if (v1 <= 0.0) {
678
        ((UINT8 *)out)[3] = 0;
679
    } else if (v1 >= 255.0) {
680
        ((UINT8 *)out)[3] = 255;
681
    } else {
682
        ((UINT8 *)out)[3] = (UINT8)v1;
683
    }
684
    return 1;
685
}
686

687
static int
688
bicubic_filter32RGB(void *out, Imaging im, double xin, double yin) {
689
    int b;
690
    BICUBIC_HEAD(UINT8);
691
    for (b = 0; b < im->bands; b++) {
692
        BICUBIC_BODY(UINT8, im->image, 4, b);
693
        if (v1 <= 0.0) {
694
            ((UINT8 *)out)[b] = 0;
695
        } else if (v1 >= 255.0) {
696
            ((UINT8 *)out)[b] = 255;
697
        } else {
698
            ((UINT8 *)out)[b] = (UINT8)v1;
699
        }
700
    }
701
    return 1;
702
}
703

704
#undef BICUBIC
705
#undef BICUBIC_HEAD
706
#undef BICUBIC_BODY
707

708
static ImagingTransformFilter
709
getfilter(Imaging im, int filterid) {
710
    switch (filterid) {
711
        case IMAGING_TRANSFORM_NEAREST:
712
            if (im->image8) {
713
                switch (im->type) {
714
                    case IMAGING_TYPE_UINT8:
715
                        return nearest_filter8;
716
                    case IMAGING_TYPE_SPECIAL:
717
                        switch (im->pixelsize) {
718
                            case 1:
719
                                return nearest_filter8;
720
                            case 2:
721
                                return nearest_filter16;
722
                            case 4:
723
                                return nearest_filter32;
724
                        }
725
                }
726
            } else {
727
                return nearest_filter32;
728
            }
729
            break;
730
        case IMAGING_TRANSFORM_BILINEAR:
731
            if (im->image8) {
732
                return bilinear_filter8;
733
            } else if (im->image32) {
734
                switch (im->type) {
735
                    case IMAGING_TYPE_UINT8:
736
                        if (im->bands == 2) {
737
                            return bilinear_filter32LA;
738
                        } else {
739
                            return bilinear_filter32RGB;
740
                        }
741
                    case IMAGING_TYPE_INT32:
742
                        return bilinear_filter32I;
743
                    case IMAGING_TYPE_FLOAT32:
744
                        return bilinear_filter32F;
745
                }
746
            }
747
            break;
748
        case IMAGING_TRANSFORM_BICUBIC:
749
            if (im->image8) {
750
                return bicubic_filter8;
751
            } else if (im->image32) {
752
                switch (im->type) {
753
                    case IMAGING_TYPE_UINT8:
754
                        if (im->bands == 2) {
755
                            return bicubic_filter32LA;
756
                        } else {
757
                            return bicubic_filter32RGB;
758
                        }
759
                    case IMAGING_TYPE_INT32:
760
                        return bicubic_filter32I;
761
                    case IMAGING_TYPE_FLOAT32:
762
                        return bicubic_filter32F;
763
                }
764
            }
765
            break;
766
    }
767
    /* no such filter */
768
    return NULL;
769
}
770

771
/* transformation engines */
772

773
Imaging
774
ImagingGenericTransform(
775
    Imaging imOut,
776
    Imaging imIn,
777
    int x0,
778
    int y0,
779
    int x1,
780
    int y1,
781
    ImagingTransformMap transform,
782
    void *transform_data,
783
    int filterid,
784
    int fill
785
) {
786
    /* slow generic transformation.  use ImagingTransformAffine or
787
       ImagingScaleAffine where possible. */
788

789
    ImagingSectionCookie cookie;
790
    int x, y;
791
    char *out;
792
    double xx, yy;
793

794
    ImagingTransformFilter filter = getfilter(imIn, filterid);
795
    if (!filter) {
796
        return (Imaging)ImagingError_ValueError("bad filter number");
797
    }
798

799
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
800
        return (Imaging)ImagingError_ModeError();
801
    }
802

803
    ImagingCopyPalette(imOut, imIn);
804

805
    ImagingSectionEnter(&cookie);
806

807
    if (x0 < 0) {
808
        x0 = 0;
809
    }
810
    if (y0 < 0) {
811
        y0 = 0;
812
    }
813
    if (x1 > imOut->xsize) {
814
        x1 = imOut->xsize;
815
    }
816
    if (y1 > imOut->ysize) {
817
        y1 = imOut->ysize;
818
    }
819

820
    for (y = y0; y < y1; y++) {
821
        out = imOut->image[y] + x0 * imOut->pixelsize;
822
        for (x = x0; x < x1; x++) {
823
            if (!transform(&xx, &yy, x - x0, y - y0, transform_data) ||
824
                !filter(out, imIn, xx, yy)) {
825
                if (fill) {
826
                    memset(out, 0, imOut->pixelsize);
827
                }
828
            }
829
            out += imOut->pixelsize;
830
        }
831
    }
832

833
    ImagingSectionLeave(&cookie);
834

835
    return imOut;
836
}
837

838
static Imaging
839
ImagingScaleAffine(
840
    Imaging imOut, Imaging imIn, int x0, int y0, int x1, int y1, double a[6], int fill
841
) {
842
    /* scale, nearest neighbour resampling */
843

844
    ImagingSectionCookie cookie;
845
    int x, y;
846
    int xin;
847
    double xo, yo;
848
    int xmin, xmax;
849
    int *xintab;
850

851
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
852
        return (Imaging)ImagingError_ModeError();
853
    }
854

855
    ImagingCopyPalette(imOut, imIn);
856

857
    if (x0 < 0) {
858
        x0 = 0;
859
    }
860
    if (y0 < 0) {
861
        y0 = 0;
862
    }
863
    if (x1 > imOut->xsize) {
864
        x1 = imOut->xsize;
865
    }
866
    if (y1 > imOut->ysize) {
867
        y1 = imOut->ysize;
868
    }
869

870
    /* malloc check ok, uses calloc for overflow */
871
    xintab = (int *)calloc(imOut->xsize, sizeof(int));
872
    if (!xintab) {
873
        ImagingDelete(imOut);
874
        return (Imaging)ImagingError_MemoryError();
875
    }
876

877
    xo = a[2] + a[0] * 0.5;
878
    yo = a[5] + a[4] * 0.5;
879

880
    xmin = x1;
881
    xmax = x0;
882

883
    /* Pretabulate horizontal pixel positions */
884
    for (x = x0; x < x1; x++) {
885
        xin = COORD(xo);
886
        if (xin >= 0 && xin < (int)imIn->xsize) {
887
            xmax = x + 1;
888
            if (x < xmin) {
889
                xmin = x;
890
            }
891
            xintab[x] = xin;
892
        }
893
        xo += a[0];
894
    }
895

896
#define AFFINE_SCALE(pixel, image)                          \
897
    for (y = y0; y < y1; y++) {                             \
898
        int yi = COORD(yo);                                 \
899
        pixel *in, *out;                                    \
900
        out = imOut->image[y];                              \
901
        if (fill && x1 > x0) {                              \
902
            memset(out + x0, 0, (x1 - x0) * sizeof(pixel)); \
903
        }                                                   \
904
        if (yi >= 0 && yi < imIn->ysize) {                  \
905
            in = imIn->image[yi];                           \
906
            for (x = xmin; x < xmax; x++) {                 \
907
                out[x] = in[xintab[x]];                     \
908
            }                                               \
909
        }                                                   \
910
        yo += a[4];                                         \
911
    }
912

913
    ImagingSectionEnter(&cookie);
914

915
    if (imIn->image8) {
916
        AFFINE_SCALE(UINT8, image8);
917
    } else {
918
        AFFINE_SCALE(INT32, image32);
919
    }
920

921
    ImagingSectionLeave(&cookie);
922

923
#undef AFFINE_SCALE
924

925
    free(xintab);
926

927
    return imOut;
928
}
929

930
static inline int
931
check_fixed(double a[6], int x, int y) {
932
    return (
933
        fabs(x * a[0] + y * a[1] + a[2]) < 32768.0 &&
934
        fabs(x * a[3] + y * a[4] + a[5]) < 32768.0
935
    );
936
}
937

938
static inline Imaging
939
affine_fixed(
940
    Imaging imOut,
941
    Imaging imIn,
942
    int x0,
943
    int y0,
944
    int x1,
945
    int y1,
946
    double a[6],
947
    int filterid,
948
    int fill
949
) {
950
    /* affine transform, nearest neighbour resampling, fixed point
951
       arithmetics */
952

953
    ImagingSectionCookie cookie;
954
    int x, y;
955
    int xin, yin;
956
    int xsize, ysize;
957
    int xx, yy;
958
    int a0, a1, a2, a3, a4, a5;
959

960
    ImagingCopyPalette(imOut, imIn);
961

962
    xsize = (int)imIn->xsize;
963
    ysize = (int)imIn->ysize;
964

965
/* use 16.16 fixed point arithmetics */
966
#define FIX(v) FLOOR((v) * 65536.0 + 0.5)
967

968
    a0 = FIX(a[0]);
969
    a1 = FIX(a[1]);
970
    a3 = FIX(a[3]);
971
    a4 = FIX(a[4]);
972
    a2 = FIX(a[2] + a[0] * 0.5 + a[1] * 0.5);
973
    a5 = FIX(a[5] + a[3] * 0.5 + a[4] * 0.5);
974

975
#undef FIX
976

977
#define AFFINE_TRANSFORM_FIXED(pixel, image)                \
978
    for (y = y0; y < y1; y++) {                             \
979
        pixel *out;                                         \
980
        xx = a2;                                            \
981
        yy = a5;                                            \
982
        out = imOut->image[y];                              \
983
        if (fill && x1 > x0) {                              \
984
            memset(out + x0, 0, (x1 - x0) * sizeof(pixel)); \
985
        }                                                   \
986
        for (x = x0; x < x1; x++, out++) {                  \
987
            xin = xx >> 16;                                 \
988
            if (xin >= 0 && xin < xsize) {                  \
989
                yin = yy >> 16;                             \
990
                if (yin >= 0 && yin < ysize) {              \
991
                    *out = imIn->image[yin][xin];           \
992
                }                                           \
993
            }                                               \
994
            xx += a0;                                       \
995
            yy += a3;                                       \
996
        }                                                   \
997
        a2 += a1;                                           \
998
        a5 += a4;                                           \
999
    }
1000

1001
    ImagingSectionEnter(&cookie);
1002

1003
    if (imIn->image8) {
1004
        AFFINE_TRANSFORM_FIXED(UINT8, image8)
1005
    } else {
1006
        AFFINE_TRANSFORM_FIXED(INT32, image32)
1007
    }
1008

1009
    ImagingSectionLeave(&cookie);
1010

1011
#undef AFFINE_TRANSFORM_FIXED
1012

1013
    return imOut;
1014
}
1015

1016
Imaging
1017
ImagingTransformAffine(
1018
    Imaging imOut,
1019
    Imaging imIn,
1020
    int x0,
1021
    int y0,
1022
    int x1,
1023
    int y1,
1024
    double a[6],
1025
    int filterid,
1026
    int fill
1027
) {
1028
    /* affine transform, nearest neighbour resampling, floating point
1029
       arithmetics*/
1030

1031
    ImagingSectionCookie cookie;
1032
    int x, y;
1033
    int xin, yin;
1034
    int xsize, ysize;
1035
    double xx, yy;
1036
    double xo, yo;
1037

1038
    if (filterid || imIn->type == IMAGING_TYPE_SPECIAL) {
1039
        return ImagingGenericTransform(
1040
            imOut, imIn, x0, y0, x1, y1, affine_transform, a, filterid, fill
1041
        );
1042
    }
1043

1044
    if (a[1] == 0 && a[3] == 0) {
1045
        /* Scaling */
1046
        return ImagingScaleAffine(imOut, imIn, x0, y0, x1, y1, a, fill);
1047
    }
1048

1049
    if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0) {
1050
        return (Imaging)ImagingError_ModeError();
1051
    }
1052

1053
    if (x0 < 0) {
1054
        x0 = 0;
1055
    }
1056
    if (y0 < 0) {
1057
        y0 = 0;
1058
    }
1059
    if (x1 > imOut->xsize) {
1060
        x1 = imOut->xsize;
1061
    }
1062
    if (y1 > imOut->ysize) {
1063
        y1 = imOut->ysize;
1064
    }
1065

1066
    /* translate all four corners to check if they are within the
1067
       range that can be represented by the fixed point arithmetics */
1068

1069
    if (check_fixed(a, 0, 0) && check_fixed(a, x1 - x0, y1 - y0) &&
1070
        check_fixed(a, 0, y1 - y0) && check_fixed(a, x1 - x0, 0)) {
1071
        return affine_fixed(imOut, imIn, x0, y0, x1, y1, a, filterid, fill);
1072
    }
1073

1074
    /* FIXME: cannot really think of any reasonable case when the
1075
       following code is used.  maybe we should fall back on the slow
1076
       generic transform engine in this case? */
1077

1078
    ImagingCopyPalette(imOut, imIn);
1079

1080
    xsize = (int)imIn->xsize;
1081
    ysize = (int)imIn->ysize;
1082

1083
    xo = a[2] + a[1] * 0.5 + a[0] * 0.5;
1084
    yo = a[5] + a[4] * 0.5 + a[3] * 0.5;
1085

1086
#define AFFINE_TRANSFORM(pixel, image)                      \
1087
    for (y = y0; y < y1; y++) {                             \
1088
        pixel *out;                                         \
1089
        xx = xo;                                            \
1090
        yy = yo;                                            \
1091
        out = imOut->image[y];                              \
1092
        if (fill && x1 > x0) {                              \
1093
            memset(out + x0, 0, (x1 - x0) * sizeof(pixel)); \
1094
        }                                                   \
1095
        for (x = x0; x < x1; x++, out++) {                  \
1096
            xin = COORD(xx);                                \
1097
            if (xin >= 0 && xin < xsize) {                  \
1098
                yin = COORD(yy);                            \
1099
                if (yin >= 0 && yin < ysize) {              \
1100
                    *out = imIn->image[yin][xin];           \
1101
                }                                           \
1102
            }                                               \
1103
            xx += a[0];                                     \
1104
            yy += a[3];                                     \
1105
        }                                                   \
1106
        xo += a[1];                                         \
1107
        yo += a[4];                                         \
1108
    }
1109

1110
    ImagingSectionEnter(&cookie);
1111

1112
    if (imIn->image8) {
1113
        AFFINE_TRANSFORM(UINT8, image8)
1114
    } else {
1115
        AFFINE_TRANSFORM(INT32, image32)
1116
    }
1117

1118
    ImagingSectionLeave(&cookie);
1119

1120
#undef AFFINE_TRANSFORM
1121

1122
    return imOut;
1123
}
1124

1125
Imaging
1126
ImagingTransform(
1127
    Imaging imOut,
1128
    Imaging imIn,
1129
    int method,
1130
    int x0,
1131
    int y0,
1132
    int x1,
1133
    int y1,
1134
    double a[8],
1135
    int filterid,
1136
    int fill
1137
) {
1138
    ImagingTransformMap transform;
1139

1140
    switch (method) {
1141
        case IMAGING_TRANSFORM_AFFINE:
1142
            return ImagingTransformAffine(
1143
                imOut, imIn, x0, y0, x1, y1, a, filterid, fill
1144
            );
1145
            break;
1146
        case IMAGING_TRANSFORM_PERSPECTIVE:
1147
            transform = perspective_transform;
1148
            break;
1149
        case IMAGING_TRANSFORM_QUAD:
1150
            transform = quad_transform;
1151
            break;
1152
        default:
1153
            return (Imaging)ImagingError_ValueError("bad transform method");
1154
    }
1155

1156
    return ImagingGenericTransform(
1157
        imOut, imIn, x0, y0, x1, y1, transform, a, filterid, fill
1158
    );
1159
}
1160