Web   ·   Wiki   ·   Activities   ·   Blog   ·   Lists   ·   Chat   ·   Meeting   ·   Bugs   ·   Git   ·   Translate   ·   Archive   ·   People   ·   Donate
summaryrefslogtreecommitdiffstats
path: root/Imaging/libImaging/Geometry.c
diff options
context:
space:
mode:
Diffstat (limited to 'Imaging/libImaging/Geometry.c')
-rw-r--r--Imaging/libImaging/Geometry.c957
1 files changed, 957 insertions, 0 deletions
diff --git a/Imaging/libImaging/Geometry.c b/Imaging/libImaging/Geometry.c
new file mode 100644
index 0000000..6b623df
--- /dev/null
+++ b/Imaging/libImaging/Geometry.c
@@ -0,0 +1,957 @@
+/*
+ * The Python Imaging Library
+ * $Id: Geometry.c 2308 2005-03-02 12:00:55Z fredrik $
+ *
+ * the imaging geometry methods
+ *
+ * history:
+ * 1995-06-15 fl Created
+ * 1996-04-15 fl Changed origin
+ * 1996-05-18 fl Fixed rotate90/270 for rectangular images
+ * 1996-05-27 fl Added general purpose transform
+ * 1996-11-22 fl Don't crash when resizing from outside source image
+ * 1997-08-09 fl Fixed rounding error in resize
+ * 1998-09-21 fl Incorporated transformation patches (from Zircon #2)
+ * 1998-09-22 fl Added bounding box to transform engines
+ * 1999-02-03 fl Fixed bicubic filtering for RGB images
+ * 1999-02-16 fl Added fixed-point version of affine transform
+ * 2001-03-28 fl Fixed transform(EXTENT) for xoffset < 0
+ * 2003-03-10 fl Compiler tweaks
+ * 2004-09-19 fl Fixed bilinear/bicubic filtering of LA images
+ *
+ * Copyright (c) 1997-2003 by Secret Labs AB
+ * Copyright (c) 1995-1997 by Fredrik Lundh
+ *
+ * See the README file for information on usage and redistribution.
+ */
+
+#include "Imaging.h"
+
+/* Undef if you don't need resampling filters */
+#define WITH_FILTERS
+
+#define COORD(v) ((v) < 0.0 ? -1 : ((int)(v)))
+#define FLOOR(v) ((v) < 0.0 ? ((int)floor(v)) : ((int)(v)))
+
+/* -------------------------------------------------------------------- */
+/* Transpose operations */
+
+Imaging
+ImagingFlipLeftRight(Imaging imOut, Imaging imIn)
+{
+ ImagingSectionCookie cookie;
+ int x, y, xr;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+ if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize)
+ return (Imaging) ImagingError_Mismatch();
+
+ ImagingCopyInfo(imOut, imIn);
+
+#define FLIP_HORIZ(image)\
+ for (y = 0; y < imIn->ysize; y++) {\
+ xr = imIn->xsize-1;\
+ for (x = 0; x < imIn->xsize; x++, xr--)\
+ imOut->image[y][x] = imIn->image[y][xr];\
+ }
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8)
+ FLIP_HORIZ(image8)
+ else
+ FLIP_HORIZ(image32)
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+
+Imaging
+ImagingFlipTopBottom(Imaging imOut, Imaging imIn)
+{
+ ImagingSectionCookie cookie;
+ int y, yr;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+ if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize)
+ return (Imaging) ImagingError_Mismatch();
+
+ ImagingCopyInfo(imOut, imIn);
+
+ ImagingSectionEnter(&cookie);
+
+ yr = imIn->ysize-1;
+ for (y = 0; y < imIn->ysize; y++, yr--)
+ memcpy(imOut->image[yr], imIn->image[y], imIn->linesize);
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+
+Imaging
+ImagingRotate90(Imaging imOut, Imaging imIn)
+{
+ ImagingSectionCookie cookie;
+ int x, y, xr;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+ if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize)
+ return (Imaging) ImagingError_Mismatch();
+
+ ImagingCopyInfo(imOut, imIn);
+
+#define ROTATE_90(image)\
+ for (y = 0; y < imIn->ysize; y++) {\
+ xr = imIn->xsize-1;\
+ for (x = 0; x < imIn->xsize; x++, xr--)\
+ imOut->image[xr][y] = imIn->image[y][x];\
+ }
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8)
+ ROTATE_90(image8)
+ else
+ ROTATE_90(image32)
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+
+Imaging
+ImagingRotate180(Imaging imOut, Imaging imIn)
+{
+ ImagingSectionCookie cookie;
+ int x, y, xr, yr;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+ if (imIn->xsize != imOut->xsize || imIn->ysize != imOut->ysize)
+ return (Imaging) ImagingError_Mismatch();
+
+ ImagingCopyInfo(imOut, imIn);
+
+ yr = imIn->ysize-1;
+
+#define ROTATE_180(image)\
+ for (y = 0; y < imIn->ysize; y++, yr--) {\
+ xr = imIn->xsize-1;\
+ for (x = 0; x < imIn->xsize; x++, xr--)\
+ imOut->image[y][x] = imIn->image[yr][xr];\
+ }
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8)
+ ROTATE_180(image8)
+ else
+ ROTATE_180(image32)
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+
+Imaging
+ImagingRotate270(Imaging imOut, Imaging imIn)
+{
+ ImagingSectionCookie cookie;
+ int x, y, yr;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+ if (imIn->xsize != imOut->ysize || imIn->ysize != imOut->xsize)
+ return (Imaging) ImagingError_Mismatch();
+
+ ImagingCopyInfo(imOut, imIn);
+
+ yr = imIn->ysize - 1;
+
+#define ROTATE_270(image)\
+ for (y = 0; y < imIn->ysize; y++, yr--)\
+ for (x = 0; x < imIn->xsize; x++)\
+ imOut->image[x][y] = imIn->image[yr][x];
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8)
+ ROTATE_270(image8)
+ else
+ ROTATE_270(image32)
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+
+/* -------------------------------------------------------------------- */
+/* Transforms */
+
+/* transform primitives (ImagingTransformMap) */
+
+static int
+affine_transform(double* xin, double* yin, int x, int y, void* data)
+{
+ /* full moon tonight. your compiler will generate bogus code
+ for simple expressions, unless you reorganize the code, or
+ install Service Pack 3 */
+
+ double* a = (double*) data;
+ double a0 = a[0]; double a1 = a[1]; double a2 = a[2];
+ double a3 = a[3]; double a4 = a[4]; double a5 = a[5];
+
+ xin[0] = a0 + a1*x + a2*y;
+ yin[0] = a3 + a4*x + a5*y;
+
+ return 1;
+}
+
+static int
+perspective_transform(double* xin, double* yin, int x, int y, void* data)
+{
+ double* a = (double*) data;
+ double a0 = a[0]; double a1 = a[1]; double a2 = a[2];
+ double a3 = a[3]; double a4 = a[4]; double a5 = a[5];
+ double a6 = a[6]; double a7 = a[7];
+
+ xin[0] = (a0 + a1*x + a2*y) / (a6*x + a7*y + 1);
+ yin[0] = (a3 + a4*x + a5*y) / (a6*x + a7*y + 1);
+
+ return 1;
+}
+
+static int
+quadratic_transform(double* xin, double* yin, int x, int y, void* data)
+{
+ double* a = (double*) data;
+
+ double a0 = a[0]; double a1 = a[1]; double a2 = a[2]; double a3 = a[3];
+ double a4 = a[4]; double a5 = a[5]; double a6 = a[6]; double a7 = a[7];
+ double a8 = a[8]; double a9 = a[9]; double a10 = a[10]; double a11 = a[11];
+
+ xin[0] = a0 + a1*x + a2*y + a3*x*x + a4*x*y + a5*y*y;
+ yin[0] = a6 + a7*x + a8*y + a9*x*x + a10*x*y + a11*y*y;
+
+ return 1;
+}
+
+static int
+quad_transform(double* xin, double* yin, int x, int y, void* data)
+{
+ /* quad warp: map quadrilateral to rectangle */
+
+ double* a = (double*) data;
+ double a0 = a[0]; double a1 = a[1]; double a2 = a[2]; double a3 = a[3];
+ double a4 = a[4]; double a5 = a[5]; double a6 = a[6]; double a7 = a[7];
+
+ xin[0] = a0 + a1*x + a2*y + a3*x*y;
+ yin[0] = a4 + a5*x + a6*y + a7*x*y;
+
+ return 1;
+}
+
+/* transform filters (ImagingTransformFilter) */
+
+#ifdef WITH_FILTERS
+
+static int
+nearest_filter8(void* out, Imaging im, double xin, double yin, void* data)
+{
+ int x = COORD(xin);
+ int y = COORD(yin);
+ if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize)
+ return 0;
+ ((UINT8*)out)[0] = im->image8[y][x];
+ return 1;
+}
+
+static int
+nearest_filter16(void* out, Imaging im, double xin, double yin, void* data)
+{
+ int x = COORD(xin);
+ int y = COORD(yin);
+ if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize)
+ return 0;
+ ((INT16*)out)[0] = ((INT16*)(im->image8[y]))[x];
+ return 1;
+}
+
+static int
+nearest_filter32(void* out, Imaging im, double xin, double yin, void* data)
+{
+ int x = COORD(xin);
+ int y = COORD(yin);
+ if (x < 0 || x >= im->xsize || y < 0 || y >= im->ysize)
+ return 0;
+ ((INT32*)out)[0] = im->image32[y][x];
+ return 1;
+}
+
+#define XCLIP(im, x) ( ((x) < 0) ? 0 : ((x) < im->xsize) ? (x) : im->xsize-1 )
+#define YCLIP(im, y) ( ((y) < 0) ? 0 : ((y) < im->ysize) ? (y) : im->ysize-1 )
+
+#define BILINEAR(v, a, b, d)\
+ (v = (a) + ( (b) - (a) ) * (d))
+
+#define BILINEAR_HEAD(type)\
+ int x, y;\
+ int x0, x1;\
+ double v1, v2;\
+ double dx, dy;\
+ type* in;\
+ if (xin < 0.0 || xin >= im->xsize || yin < 0.0 || yin >= im->ysize)\
+ return 0;\
+ xin -= 0.5;\
+ yin -= 0.5;\
+ x = FLOOR(xin);\
+ y = FLOOR(yin);\
+ dx = xin - x;\
+ dy = yin - y;
+
+#define BILINEAR_BODY(type, image, step, offset) {\
+ in = (type*) ((image)[YCLIP(im, y)] + offset);\
+ x0 = XCLIP(im, x+0)*step;\
+ x1 = XCLIP(im, x+1)*step;\
+ BILINEAR(v1, in[x0], in[x1], dx);\
+ if (y+1 >= 0 && y+1 < im->ysize) {\
+ in = (type*) ((image)[y+1] + offset);\
+ BILINEAR(v2, in[x0], in[x1], dx);\
+ } else\
+ v2 = v1;\
+ BILINEAR(v1, v1, v2, dy);\
+}
+
+static int
+bilinear_filter8(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BILINEAR_HEAD(UINT8);
+ BILINEAR_BODY(UINT8, im->image8, 1, 0);
+ ((UINT8*)out)[0] = (UINT8) v1;
+ return 1;
+}
+
+static int
+bilinear_filter32I(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BILINEAR_HEAD(INT32);
+ BILINEAR_BODY(INT32, im->image32, 1, 0);
+ ((INT32*)out)[0] = (INT32) v1;
+ return 1;
+}
+
+static int
+bilinear_filter32F(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BILINEAR_HEAD(FLOAT32);
+ BILINEAR_BODY(FLOAT32, im->image32, 1, 0);
+ ((FLOAT32*)out)[0] = (FLOAT32) v1;
+ return 1;
+}
+
+static int
+bilinear_filter32LA(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BILINEAR_HEAD(UINT8);
+ BILINEAR_BODY(UINT8, im->image, 4, 0);
+ ((UINT8*)out)[0] = (UINT8) v1;
+ ((UINT8*)out)[1] = (UINT8) v1;
+ ((UINT8*)out)[2] = (UINT8) v1;
+ BILINEAR_BODY(UINT8, im->image, 4, 3);
+ ((UINT8*)out)[3] = (UINT8) v1;
+ return 1;
+}
+
+static int
+bilinear_filter32RGB(void* out, Imaging im, double xin, double yin, void* data)
+{
+ int b;
+ BILINEAR_HEAD(UINT8);
+ for (b = 0; b < im->bands; b++) {
+ BILINEAR_BODY(UINT8, im->image, 4, b);
+ ((UINT8*)out)[b] = (UINT8) v1;
+ }
+ return 1;
+}
+
+#define BICUBIC(v, v1, v2, v3, v4, d) {\
+ double p1 = v2;\
+ double p2 = -v1 + v3;\
+ double p3 = 2*(v1 - v2) + v3 - v4;\
+ double p4 = -v1 + v2 - v3 + v4;\
+ v = p1 + (d)*(p2 + (d)*(p3 + (d)*p4));\
+}
+
+#define BICUBIC_HEAD(type)\
+ int x = FLOOR(xin);\
+ int y = FLOOR(yin);\
+ int x0, x1, x2, x3;\
+ double v1, v2, v3, v4;\
+ double dx, dy;\
+ type* in;\
+ if (xin < 0.0 || xin >= im->xsize || yin < 0.0 || yin >= im->ysize)\
+ return 0;\
+ xin -= 0.5;\
+ yin -= 0.5;\
+ x = FLOOR(xin);\
+ y = FLOOR(yin);\
+ dx = xin - x;\
+ dy = yin - y;\
+ x--; y--;
+
+#define BICUBIC_BODY(type, image, step, offset) {\
+ in = (type*) ((image)[YCLIP(im, y)] + offset);\
+ x0 = XCLIP(im, x+0)*step;\
+ x1 = XCLIP(im, x+1)*step;\
+ x2 = XCLIP(im, x+2)*step;\
+ x3 = XCLIP(im, x+3)*step;\
+ BICUBIC(v1, in[x0], in[x1], in[x2], in[x3], dx);\
+ if (y+1 >= 0 && y+1 < im->ysize) {\
+ in = (type*) ((image)[y+1] + offset);\
+ BICUBIC(v2, in[x0], in[x1], in[x2], in[x3], dx);\
+ } else\
+ v2 = v1;\
+ if (y+2 >= 0 && y+2 < im->ysize) {\
+ in = (type*) ((image)[y+2] + offset);\
+ BICUBIC(v3, in[x0], in[x1], in[x2], in[x3], dx);\
+ } else\
+ v3 = v2;\
+ if (y+3 >= 0 && y+3 < im->ysize) {\
+ in = (type*) ((image)[y+3] + offset);\
+ BICUBIC(v4, in[x0], in[x1], in[x2], in[x3], dx);\
+ } else\
+ v4 = v3;\
+ BICUBIC(v1, v1, v2, v3, v4, dy);\
+}
+
+
+static int
+bicubic_filter8(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BICUBIC_HEAD(UINT8);
+ BICUBIC_BODY(UINT8, im->image8, 1, 0);
+ if (v1 <= 0.0)
+ ((UINT8*)out)[0] = 0;
+ else if (v1 >= 255.0)
+ ((UINT8*)out)[0] = 255;
+ else
+ ((UINT8*)out)[0] = (UINT8) v1;
+ return 1;
+}
+
+static int
+bicubic_filter32I(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BICUBIC_HEAD(INT32);
+ BICUBIC_BODY(INT32, im->image32, 1, 0);
+ ((INT32*)out)[0] = (INT32) v1;
+ return 1;
+}
+
+static int
+bicubic_filter32F(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BICUBIC_HEAD(FLOAT32);
+ BICUBIC_BODY(FLOAT32, im->image32, 1, 0);
+ ((FLOAT32*)out)[0] = (FLOAT32) v1;
+ return 1;
+}
+
+static int
+bicubic_filter32LA(void* out, Imaging im, double xin, double yin, void* data)
+{
+ BICUBIC_HEAD(UINT8);
+ BICUBIC_BODY(UINT8, im->image, 4, 0);
+ if (v1 <= 0.0) {
+ ((UINT8*)out)[0] = 0;
+ ((UINT8*)out)[1] = 0;
+ ((UINT8*)out)[2] = 0;
+ } else if (v1 >= 255.0) {
+ ((UINT8*)out)[0] = 255;
+ ((UINT8*)out)[1] = 255;
+ ((UINT8*)out)[2] = 255;
+ } else {
+ ((UINT8*)out)[0] = (UINT8) v1;
+ ((UINT8*)out)[1] = (UINT8) v1;
+ ((UINT8*)out)[2] = (UINT8) v1;
+ }
+ BICUBIC_BODY(UINT8, im->image, 4, 3);
+ if (v1 <= 0.0)
+ ((UINT8*)out)[3] = 0;
+ else if (v1 >= 255.0)
+ ((UINT8*)out)[3] = 255;
+ else
+ ((UINT8*)out)[3] = (UINT8) v1;
+ return 1;
+}
+
+static int
+bicubic_filter32RGB(void* out, Imaging im, double xin, double yin, void* data)
+{
+ int b;
+ BICUBIC_HEAD(UINT8);
+ for (b = 0; b < im->bands; b++) {
+ BICUBIC_BODY(UINT8, im->image, 4, b);
+ if (v1 <= 0.0)
+ ((UINT8*)out)[b] = 0;
+ else if (v1 >= 255.0)
+ ((UINT8*)out)[b] = 255;
+ else
+ ((UINT8*)out)[b] = (UINT8) v1;
+ }
+ return 1;
+}
+
+static ImagingTransformFilter
+getfilter(Imaging im, int filterid)
+{
+ switch (filterid) {
+ case IMAGING_TRANSFORM_NEAREST:
+ if (im->image8)
+ switch (im->type) {
+ case IMAGING_TYPE_UINT8:
+ return (ImagingTransformFilter) nearest_filter8;
+ case IMAGING_TYPE_SPECIAL:
+ switch (im->pixelsize) {
+ case 1:
+ return (ImagingTransformFilter) nearest_filter8;
+ case 2:
+ return (ImagingTransformFilter) nearest_filter16;
+ case 4:
+ return (ImagingTransformFilter) nearest_filter32;
+ }
+ }
+ else
+ return (ImagingTransformFilter) nearest_filter32;
+ break;
+ case IMAGING_TRANSFORM_BILINEAR:
+ if (im->image8)
+ return (ImagingTransformFilter) bilinear_filter8;
+ else if (im->image32) {
+ switch (im->type) {
+ case IMAGING_TYPE_UINT8:
+ if (im->bands == 2)
+ return (ImagingTransformFilter) bilinear_filter32LA;
+ else
+ return (ImagingTransformFilter) bilinear_filter32RGB;
+ case IMAGING_TYPE_INT32:
+ return (ImagingTransformFilter) bilinear_filter32I;
+ case IMAGING_TYPE_FLOAT32:
+ return (ImagingTransformFilter) bilinear_filter32F;
+ }
+ }
+ break;
+ case IMAGING_TRANSFORM_BICUBIC:
+ if (im->image8)
+ return (ImagingTransformFilter) bicubic_filter8;
+ else if (im->image32) {
+ switch (im->type) {
+ case IMAGING_TYPE_UINT8:
+ if (im->bands == 2)
+ return (ImagingTransformFilter) bicubic_filter32LA;
+ else
+ return (ImagingTransformFilter) bicubic_filter32RGB;
+ case IMAGING_TYPE_INT32:
+ return (ImagingTransformFilter) bicubic_filter32I;
+ case IMAGING_TYPE_FLOAT32:
+ return (ImagingTransformFilter) bicubic_filter32F;
+ }
+ }
+ break;
+ }
+ /* no such filter */
+ return NULL;
+}
+
+#else
+#define getfilter(im, id) NULL
+#endif
+
+/* transformation engines */
+
+Imaging
+ImagingTransform(
+ Imaging imOut, Imaging imIn, int x0, int y0, int x1, int y1,
+ ImagingTransformMap transform, void* transform_data,
+ ImagingTransformFilter filter, void* filter_data,
+ int fill)
+{
+ /* slow generic transformation. use ImagingTransformAffine or
+ ImagingScaleAffine where possible. */
+
+ ImagingSectionCookie cookie;
+ int x, y;
+ char *out;
+ double xx, yy;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+
+ ImagingCopyInfo(imOut, imIn);
+
+ ImagingSectionEnter(&cookie);
+
+ if (x0 < 0)
+ x0 = 0;
+ if (y0 < 0)
+ y0 = 0;
+ if (x1 > imOut->xsize)
+ x1 = imOut->xsize;
+ if (y1 > imOut->ysize)
+ y1 = imOut->ysize;
+
+ for (y = y0; y < y1; y++) {
+ out = imOut->image[y] + x0*imOut->pixelsize;
+ for (x = x0; x < x1; x++) {
+ if (!transform(&xx, &yy, x-x0, y-y0, transform_data) ||
+ !filter(out, imIn, xx, yy, filter_data)) {
+ if (fill)
+ memset(out, 0, imOut->pixelsize);
+ }
+ out += imOut->pixelsize;
+ }
+ }
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+static Imaging
+ImagingScaleAffine(Imaging imOut, Imaging imIn,
+ int x0, int y0, int x1, int y1,
+ double a[6], int fill)
+{
+ /* scale, nearest neighbour resampling */
+
+ ImagingSectionCookie cookie;
+ int x, y;
+ int xin;
+ double xo, yo;
+ int xmin, xmax;
+ int *xintab;
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+
+ ImagingCopyInfo(imOut, imIn);
+
+ if (x0 < 0)
+ x0 = 0;
+ if (y0 < 0)
+ y0 = 0;
+ if (x1 > imOut->xsize)
+ x1 = imOut->xsize;
+ if (y1 > imOut->ysize)
+ y1 = imOut->ysize;
+
+ xintab = (int*) malloc(imOut->xsize * sizeof(int));
+ if (!xintab) {
+ ImagingDelete(imOut);
+ return (Imaging) ImagingError_MemoryError();
+ }
+
+ xo = a[0];
+ yo = a[3];
+
+ xmin = x1;
+ xmax = x0;
+
+ /* Pretabulate horizontal pixel positions */
+ for (x = x0; x < x1; x++) {
+ xin = COORD(xo);
+ if (xin >= 0 && xin < (int) imIn->xsize) {
+ xmax = x+1;
+ if (x < xmin)
+ xmin = x;
+ xintab[x] = xin;
+ }
+ xo += a[1];
+ }
+
+#define AFFINE_SCALE(pixel, image)\
+ for (y = y0; y < y1; y++) {\
+ int yi = COORD(yo);\
+ pixel *in, *out;\
+ out = imOut->image[y];\
+ if (fill && x1 > x0)\
+ memset(out+x0, 0, (x1-x0)*sizeof(pixel));\
+ if (yi >= 0 && yi < imIn->ysize) {\
+ in = imIn->image[yi];\
+ for (x = xmin; x < xmax; x++)\
+ out[x] = in[xintab[x]];\
+ }\
+ yo += a[5];\
+ }
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8) {
+ AFFINE_SCALE(UINT8, image8);
+ } else {
+ AFFINE_SCALE(INT32, image32);
+ }
+
+ ImagingSectionLeave(&cookie);
+
+ free(xintab);
+
+ return imOut;
+}
+
+static inline int
+check_fixed(double a[6], int x, int y)
+{
+ return (fabs(a[0] + x*a[1] + y*a[2]) < 32768.0 &&
+ fabs(a[3] + x*a[4] + y*a[5]) < 32768.0);
+}
+
+static inline Imaging
+affine_fixed(Imaging imOut, Imaging imIn,
+ int x0, int y0, int x1, int y1,
+ double a[6], int filterid, int fill)
+{
+ /* affine transform, nearest neighbour resampling, fixed point
+ arithmetics */
+
+ int x, y;
+ int xin, yin;
+ int xsize, ysize;
+ int xx, yy;
+ int a0, a1, a2, a3, a4, a5;
+
+ ImagingCopyInfo(imOut, imIn);
+
+ xsize = (int) imIn->xsize;
+ ysize = (int) imIn->ysize;
+
+/* use 16.16 fixed point arithmetics */
+#define FIX(v) FLOOR((v)*65536.0 + 0.5)
+
+ a0 = FIX(a[0]); a1 = FIX(a[1]); a2 = FIX(a[2]);
+ a3 = FIX(a[3]); a4 = FIX(a[4]); a5 = FIX(a[5]);
+
+#define AFFINE_TRANSFORM_FIXED(pixel, image)\
+ for (y = y0; y < y1; y++) {\
+ pixel *out;\
+ xx = a0;\
+ yy = a3;\
+ out = imOut->image[y];\
+ if (fill && x1 > x0)\
+ memset(out+x0, 0, (x1-x0)*sizeof(pixel));\
+ for (x = x0; x < x1; x++, out++) {\
+ xin = xx >> 16;\
+ if (xin >= 0 && xin < xsize) {\
+ yin = yy >> 16;\
+ if (yin >= 0 && yin < ysize)\
+ *out = imIn->image[yin][xin];\
+ }\
+ xx += a1;\
+ yy += a4;\
+ }\
+ a0 += a2;\
+ a3 += a5;\
+ }
+
+ if (imIn->image8)
+ AFFINE_TRANSFORM_FIXED(UINT8, image8)
+ else
+ AFFINE_TRANSFORM_FIXED(INT32, image32)
+
+ return imOut;
+}
+
+Imaging
+ImagingTransformAffine(Imaging imOut, Imaging imIn,
+ int x0, int y0, int x1, int y1,
+ double a[6], int filterid, int fill)
+{
+ /* affine transform, nearest neighbour resampling, floating point
+ arithmetics*/
+
+ ImagingSectionCookie cookie;
+ int x, y;
+ int xin, yin;
+ int xsize, ysize;
+ double xx, yy;
+ double xo, yo;
+
+ if (filterid || imIn->type == IMAGING_TYPE_SPECIAL) {
+ /* Filtered transform */
+ ImagingTransformFilter filter = getfilter(imIn, filterid);
+ if (!filter)
+ return (Imaging) ImagingError_ValueError("unknown filter");
+ return ImagingTransform(
+ imOut, imIn,
+ x0, y0, x1, y1,
+ affine_transform, a,
+ filter, NULL, fill);
+ }
+
+ if (a[2] == 0 && a[4] == 0)
+ /* Scaling */
+ return ImagingScaleAffine(imOut, imIn, x0, y0, x1, y1, a, fill);
+
+ if (!imOut || !imIn || strcmp(imIn->mode, imOut->mode) != 0)
+ return (Imaging) ImagingError_ModeError();
+
+ if (x0 < 0)
+ x0 = 0;
+ if (y0 < 0)
+ y0 = 0;
+ if (x1 > imOut->xsize)
+ x1 = imOut->xsize;
+ if (y1 > imOut->ysize)
+ y1 = imOut->ysize;
+
+ ImagingCopyInfo(imOut, imIn);
+
+ /* translate all four corners to check if they are within the
+ range that can be represented by the fixed point arithmetics */
+
+ if (check_fixed(a, 0, 0) && check_fixed(a, x1-x0, y1-y0) &&
+ check_fixed(a, 0, y1-y0) && check_fixed(a, x1-x0, 0))
+ return affine_fixed(imOut, imIn, x0, y0, x1, y1, a, filterid, fill);
+
+ /* FIXME: cannot really think of any reasonable case when the
+ following code is used. maybe we should fall back on the slow
+ generic transform engine in this case? */
+
+ xsize = (int) imIn->xsize;
+ ysize = (int) imIn->ysize;
+
+ xo = a[0];
+ yo = a[3];
+
+#define AFFINE_TRANSFORM(pixel, image)\
+ for (y = y0; y < y1; y++) {\
+ pixel *out;\
+ xx = xo;\
+ yy = yo;\
+ out = imOut->image[y];\
+ if (fill && x1 > x0)\
+ memset(out+x0, 0, (x1-x0)*sizeof(pixel));\
+ for (x = x0; x < x1; x++, out++) {\
+ xin = COORD(xx);\
+ if (xin >= 0 && xin < xsize) {\
+ yin = COORD(yy);\
+ if (yin >= 0 && yin < ysize)\
+ *out = imIn->image[yin][xin];\
+ }\
+ xx += a[1];\
+ yy += a[4];\
+ }\
+ xo += a[2];\
+ yo += a[5];\
+ }
+
+ ImagingSectionEnter(&cookie);
+
+ if (imIn->image8)
+ AFFINE_TRANSFORM(UINT8, image8)
+ else
+ AFFINE_TRANSFORM(INT32, image32)
+
+ ImagingSectionLeave(&cookie);
+
+ return imOut;
+}
+
+Imaging
+ImagingTransformPerspective(Imaging imOut, Imaging imIn,
+ int x0, int y0, int x1, int y1,
+ double a[8], int filterid, int fill)
+{
+ ImagingTransformFilter filter = getfilter(imIn, filterid);
+ if (!filter)
+ return (Imaging) ImagingError_ValueError("bad filter number");
+
+ return ImagingTransform(
+ imOut, imIn,
+ x0, y0, x1, y1,
+ perspective_transform, a,
+ filter, NULL,
+ fill);
+}
+
+Imaging
+ImagingTransformQuad(Imaging imOut, Imaging imIn,
+ int x0, int y0, int x1, int y1,
+ double a[8], int filterid, int fill)
+{
+ ImagingTransformFilter filter = getfilter(imIn, filterid);
+ if (!filter)
+ return (Imaging) ImagingError_ValueError("bad filter number");
+
+ return ImagingTransform(
+ imOut, imIn,
+ x0, y0, x1, y1,
+ quad_transform, a,
+ filter, NULL,
+ fill);
+}
+
+/* -------------------------------------------------------------------- */
+/* Convenience functions */
+
+Imaging
+ImagingResize(Imaging imOut, Imaging imIn, int filterid)
+{
+ double a[6];
+
+ if (imOut->xsize == imIn->xsize && imOut->ysize == imIn->ysize)
+ return ImagingCopy2(imOut, imIn);
+
+ memset(a, 0, sizeof a);
+ a[1] = (double) imIn->xsize / imOut->xsize;
+ a[5] = (double) imIn->ysize / imOut->ysize;
+
+ if (!filterid && imIn->type != IMAGING_TYPE_SPECIAL)
+ return ImagingScaleAffine(
+ imOut, imIn,
+ 0, 0, imOut->xsize, imOut->ysize,
+ a, 1);
+
+ return ImagingTransformAffine(
+ imOut, imIn,
+ 0, 0, imOut->xsize, imOut->ysize,
+ a, filterid, 1);
+}
+
+Imaging
+ImagingRotate(Imaging imOut, Imaging imIn, double theta, int filterid)
+{
+ int xsize, ysize;
+ double sintheta, costheta;
+ double a[6];
+
+ /* Setup an affine transform to rotate around the image center */
+ theta = -theta * M_PI / 180.0;
+ sintheta = sin(theta);
+ costheta = cos(theta);
+
+ xsize = imOut->xsize;
+ ysize = imOut->ysize;
+
+ a[0] = -costheta * xsize/2 - sintheta * ysize/2 + xsize/2;
+ a[1] = costheta;
+ a[2] = sintheta;
+ a[3] = sintheta * xsize/2 - costheta * ysize/2 + ysize/2;
+ a[4] = -sintheta;
+ a[5] = costheta;
+
+ return ImagingTransformAffine(
+ imOut, imIn,
+ 0, 0, imOut->xsize, imOut->ysize,
+ a, filterid, 1);
+}