Pillow/libImaging/Resample.c

319 lines
9.5 KiB
C
Raw Normal View History

2010-07-31 06:52:47 +04:00
/*
* The Python Imaging Library
* $Id$
*
2015-06-02 15:50:22 +03:00
* Pillow image resampling support
2010-07-31 06:52:47 +04:00
*
* history:
* 2002-03-09 fl Created (for PIL 1.1.3)
* 2002-03-10 fl Added support for mode "F"
*
* Copyright (c) 1997-2002 by Secret Labs AB
*
* See the README file for information on usage and redistribution.
*/
#include "Imaging.h"
#include <math.h>
struct filter {
float (*filter)(float x);
float support;
};
static inline float sinc_filter(float x)
{
if (x == 0.0)
return 1.0;
x = x * M_PI;
return sin(x) / x;
}
static inline float lanczos_filter(float x)
2010-07-31 06:52:47 +04:00
{
/* truncated sinc */
2010-07-31 06:52:47 +04:00
if (-3.0 <= x && x < 3.0)
return sinc_filter(x) * sinc_filter(x/3);
return 0.0;
}
static struct filter LANCZOS = { lanczos_filter, 3.0 };
2010-07-31 06:52:47 +04:00
static inline float bilinear_filter(float x)
{
if (x < 0.0)
x = -x;
if (x < 1.0)
return 1.0-x;
return 0.0;
}
static struct filter BILINEAR = { bilinear_filter, 1.0 };
static inline float bicubic_filter(float x)
{
2014-10-24 12:46:51 +04:00
/* http://en.wikipedia.org/wiki/Bicubic_interpolation#Bicubic_convolution_algorithm */
2014-10-25 12:39:03 +04:00
#define a -0.5
2010-07-31 06:52:47 +04:00
if (x < 0.0)
x = -x;
if (x < 1.0)
2014-10-24 12:46:51 +04:00
return ((a + 2.0) * x - (a + 3.0)) * x*x + 1;
2010-07-31 06:52:47 +04:00
if (x < 2.0)
2014-10-24 12:46:51 +04:00
return (((x - 5) * x + 8) * x - 4) * a;
2010-07-31 06:52:47 +04:00
return 0.0;
#undef a
}
static struct filter BICUBIC = { bicubic_filter, 2.0 };
2014-10-25 06:04:04 +04:00
static inline UINT8 clip8(float in)
{
int out = (int) in;
if (out >= 255)
return 255;
if (out <= 0)
return 0;
return (UINT8) out;
}
/* This is work around bug in GCC prior 4.9 in 64 bit mode.
GCC generates code with partial dependency which 3 times slower.
See: http://stackoverflow.com/a/26588074/253146 */
2014-11-14 08:35:06 +03:00
#if defined(__x86_64__) && defined(__SSE__) && ! defined(__NO_INLINE__) && \
! defined(__clang__) && defined(GCC_VERSION) && (GCC_VERSION < 40900)
static float __attribute__((always_inline)) i2f(int v) {
float x;
__asm__("xorps %0, %0; cvtsi2ss %1, %0" : "=X"(x) : "r"(v) );
return x;
}
#else
static float inline i2f(int v) { return (float) v; }
#endif
2010-07-31 06:52:47 +04:00
Imaging
ImagingResampleHorizontal(Imaging imIn, int xsize, int filter)
2010-07-31 06:52:47 +04:00
{
ImagingSectionCookie cookie;
2014-10-26 00:11:39 +04:00
Imaging imOut;
2010-07-31 06:52:47 +04:00
struct filter *filterp;
float support, scale, filterscale;
float center, ww, ss, ss0, ss1, ss2, ss3;
2014-10-26 01:20:24 +04:00
int xx, yy, x, kmax, xmin, xmax;
2014-10-25 05:57:08 +04:00
int *xbounds;
float *k, *kk;
2010-07-31 06:52:47 +04:00
/* check filter */
switch (filter) {
case IMAGING_TRANSFORM_LANCZOS:
filterp = &LANCZOS;
2010-07-31 06:52:47 +04:00
break;
case IMAGING_TRANSFORM_BILINEAR:
filterp = &BILINEAR;
break;
case IMAGING_TRANSFORM_BICUBIC:
filterp = &BICUBIC;
break;
default:
return (Imaging) ImagingError_ValueError(
"unsupported resampling filter"
);
}
/* prepare for horizontal stretch */
2014-10-26 00:11:39 +04:00
filterscale = scale = (float) imIn->xsize / xsize;
2010-07-31 06:52:47 +04:00
/* determine support size (length of resampling filter) */
support = filterp->support;
if (filterscale < 1.0) {
filterscale = 1.0;
}
2010-07-31 06:52:47 +04:00
support = support * filterscale;
2014-10-25 05:37:33 +04:00
/* maximum number of coofs */
kmax = (int) ceil(support) * 2 + 1;
/* coefficient buffer */
2014-10-26 00:11:39 +04:00
kk = malloc(xsize * kmax * sizeof(float));
2014-10-25 05:37:33 +04:00
if ( ! kk)
return (Imaging) ImagingError_MemoryError();
2014-10-26 00:11:39 +04:00
xbounds = malloc(xsize * 2 * sizeof(int));
2014-10-25 05:37:33 +04:00
if ( ! xbounds) {
free(kk);
2010-07-31 06:52:47 +04:00
return (Imaging) ImagingError_MemoryError();
2014-10-25 05:37:33 +04:00
}
2014-10-26 00:11:39 +04:00
for (xx = 0; xx < xsize; xx++) {
2014-10-25 05:37:33 +04:00
k = &kk[xx * kmax];
center = (xx + 0.5) * scale;
ww = 0.0;
ss = 1.0 / filterscale;
2014-10-25 05:57:08 +04:00
xmin = (int) floor(center - support);
if (xmin < 0)
xmin = 0;
xmax = (int) ceil(center + support);
if (xmax > imIn->xsize)
xmax = imIn->xsize;
for (x = xmin; x < xmax; x++) {
2014-10-25 05:37:33 +04:00
float w = filterp->filter((x - center + 0.5) * ss) * ss;
2014-10-25 05:57:08 +04:00
k[x - xmin] = w;
2014-10-25 05:50:54 +04:00
ww += w;
2014-10-25 05:37:33 +04:00
}
2014-10-25 05:57:08 +04:00
for (x = 0; x < xmax - xmin; x++) {
2014-10-25 05:50:54 +04:00
if (ww != 0.0)
k[x] /= ww;
}
xbounds[xx * 2 + 0] = xmin;
xbounds[xx * 2 + 1] = xmax;
2014-10-25 05:37:33 +04:00
}
2010-07-31 06:52:47 +04:00
2014-10-26 00:11:39 +04:00
imOut = ImagingNew(imIn->mode, xsize, imIn->ysize);
if ( ! imOut) {
free(kk);
free(xbounds);
return NULL;
}
2010-07-31 06:52:47 +04:00
ImagingSectionEnter(&cookie);
/* horizontal stretch */
2014-10-25 05:41:38 +04:00
for (yy = 0; yy < imOut->ysize; yy++) {
if (imIn->image8) {
/* 8-bit grayscale */
2014-10-26 00:11:39 +04:00
for (xx = 0; xx < xsize; xx++) {
2014-10-25 05:50:54 +04:00
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
2014-10-25 05:41:38 +04:00
k = &kk[xx * kmax];
2014-10-25 05:50:54 +04:00
ss = 0.5;
2014-10-25 05:57:08 +04:00
for (x = xmin; x < xmax; x++)
ss += i2f(imIn->image8[yy][x]) * k[x - xmin];
2014-10-25 06:04:04 +04:00
imOut->image8[yy][xx] = clip8(ss);
2010-07-31 06:52:47 +04:00
}
} else {
switch(imIn->type) {
case IMAGING_TYPE_UINT8:
/* n-bit grayscale */
if (imIn->bands == 2) {
for (xx = 0; xx < xsize; xx++) {
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
k = &kk[xx * kmax];
2014-11-09 01:05:50 +03:00
ss0 = ss1 = 0.5;
for (x = xmin; x < xmax; x++) {
ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
ss1 += i2f((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
}
imOut->image[yy][xx*4 + 0] = clip8(ss0);
imOut->image[yy][xx*4 + 3] = clip8(ss1);
}
} else if (imIn->bands == 3) {
for (xx = 0; xx < xsize; xx++) {
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
k = &kk[xx * kmax];
ss0 = ss1 = ss2 = 0.5;
2014-10-26 01:20:24 +04:00
for (x = xmin; x < xmax; x++) {
ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
ss1 += i2f((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
ss2 += i2f((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
2014-10-26 01:20:24 +04:00
}
imOut->image[yy][xx*4 + 0] = clip8(ss0);
imOut->image[yy][xx*4 + 1] = clip8(ss1);
imOut->image[yy][xx*4 + 2] = clip8(ss2);
}
} else {
for (xx = 0; xx < xsize; xx++) {
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
k = &kk[xx * kmax];
ss0 = ss1 = ss2 = ss3 = 0.5;
2014-10-26 01:20:24 +04:00
for (x = xmin; x < xmax; x++) {
ss0 += i2f((UINT8) imIn->image[yy][x*4 + 0]) * k[x - xmin];
ss1 += i2f((UINT8) imIn->image[yy][x*4 + 1]) * k[x - xmin];
ss2 += i2f((UINT8) imIn->image[yy][x*4 + 2]) * k[x - xmin];
ss3 += i2f((UINT8) imIn->image[yy][x*4 + 3]) * k[x - xmin];
2014-10-26 01:20:24 +04:00
}
imOut->image[yy][xx*4 + 0] = clip8(ss0);
imOut->image[yy][xx*4 + 1] = clip8(ss1);
imOut->image[yy][xx*4 + 2] = clip8(ss2);
imOut->image[yy][xx*4 + 3] = clip8(ss3);
2010-07-31 06:52:47 +04:00
}
}
break;
case IMAGING_TYPE_INT32:
/* 32-bit integer */
2014-10-26 00:11:39 +04:00
for (xx = 0; xx < xsize; xx++) {
2014-10-25 05:50:54 +04:00
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
2014-10-25 05:41:38 +04:00
k = &kk[xx * kmax];
2010-07-31 06:52:47 +04:00
ss = 0.0;
2014-10-25 05:57:08 +04:00
for (x = xmin; x < xmax; x++)
ss += i2f(IMAGING_PIXEL_I(imIn, x, yy)) * k[x - xmin];
2014-10-25 05:50:54 +04:00
IMAGING_PIXEL_I(imOut, xx, yy) = (int) ss;
2010-07-31 06:52:47 +04:00
}
break;
case IMAGING_TYPE_FLOAT32:
/* 32-bit float */
2014-10-26 00:11:39 +04:00
for (xx = 0; xx < xsize; xx++) {
2014-10-25 05:50:54 +04:00
xmin = xbounds[xx * 2 + 0];
xmax = xbounds[xx * 2 + 1];
2014-10-25 05:41:38 +04:00
k = &kk[xx * kmax];
ss = 0.0;
2014-10-25 05:57:08 +04:00
for (x = xmin; x < xmax; x++)
ss += IMAGING_PIXEL_F(imIn, x, yy) * k[x - xmin];
2014-10-25 05:50:54 +04:00
IMAGING_PIXEL_F(imOut, xx, yy) = ss;
2010-07-31 06:52:47 +04:00
}
break;
}
}
2010-07-31 06:52:47 +04:00
}
ImagingSectionLeave(&cookie);
2014-10-25 05:37:33 +04:00
free(kk);
free(xbounds);
2010-07-31 06:52:47 +04:00
return imOut;
}
2014-10-26 00:11:39 +04:00
Imaging
ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
{
Imaging imTemp1, imTemp2, imTemp3;
2014-10-26 00:11:39 +04:00
Imaging imOut;
if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
return (Imaging) ImagingError_ModeError();
if (imIn->type == IMAGING_TYPE_SPECIAL)
return (Imaging) ImagingError_ModeError();
2014-10-26 00:11:39 +04:00
/* two-pass resize, first pass */
imTemp1 = ImagingResampleHorizontal(imIn, xsize, filter);
if ( ! imTemp1)
return NULL;
/* transpose image once */
2014-10-26 00:11:39 +04:00
imTemp2 = ImagingTransposeToNew(imTemp1);
ImagingDelete(imTemp1);
2014-10-26 00:11:39 +04:00
if ( ! imTemp2)
return NULL;
/* second pass */
imTemp3 = ImagingResampleHorizontal(imTemp2, ysize, filter);
ImagingDelete(imTemp2);
2014-10-26 00:11:39 +04:00
if ( ! imTemp3)
return NULL;
/* transpose result */
2014-10-26 00:11:39 +04:00
imOut = ImagingTransposeToNew(imTemp3);
ImagingDelete(imTemp3);
2014-10-26 00:11:39 +04:00
if ( ! imOut)
return NULL;
return imOut;
}