diff options
Diffstat (limited to 'Imaging/libImaging/Effects.c')
| -rw-r--r-- | Imaging/libImaging/Effects.c | 370 |
1 files changed, 370 insertions, 0 deletions
diff --git a/Imaging/libImaging/Effects.c b/Imaging/libImaging/Effects.c new file mode 100644 index 0000000..2ec08c7 --- /dev/null +++ b/Imaging/libImaging/Effects.c @@ -0,0 +1,370 @@ +/* + * The Python Imaging Library + * $Id: Effects.c 2134 2004-10-06 08:55:20Z fredrik $ + * + * various special effects and image generators + * + * history: + * 1997-05-21 fl Just for fun + * 1997-06-05 fl Added mandelbrot generator + * 2003-05-24 fl Added perlin_turbulence generator (in progress) + * + * Copyright (c) 1997-2003 by Fredrik Lundh. + * Copyright (c) 1997 by Secret Labs AB. + * + * See the README file for information on usage and redistribution. + */ + + +#include "Imaging.h" + +#include <math.h> + +Imaging +ImagingEffectMandelbrot(int xsize, int ysize, double extent[4], int quality) +{ + /* Generate a Mandelbrot set covering the given extent */ + + Imaging im; + int x, y, k; + double width, height; + double x1, y1, xi2, yi2, cr, ci, radius; + double dr, di; + + /* Check arguments */ + width = extent[2] - extent[0]; + height = extent[3] - extent[1]; + if (width < 0.0 || height < 0.0 || quality < 2) + return (Imaging) ImagingError_ValueError(NULL); + + im = ImagingNew("L", xsize, ysize); + if (!im) + return NULL; + + dr = width/(xsize-1); + di = height/(ysize-1); + + radius = 100.0; + + for (y = 0; y < ysize; y++) { + UINT8* buf = im->image8[y]; + for (x = 0; x < xsize; x++) { + x1 = y1 = xi2 = yi2 = 0.0; + cr = x*dr + extent[0]; + ci = y*di + extent[1]; + for (k = 1;; k++) { + y1 = 2*x1*y1 + ci; + x1 = xi2 - yi2 + cr; + xi2 = x1*x1; + yi2 = y1*y1; + if ((xi2 + yi2) > radius) { + buf[x] = k*255/quality; + break; + } + if (k > quality) { + buf[x] = 0; + break; + } + } + } + } + return im; +} + +Imaging +ImagingEffectNoise(int xsize, int ysize, float sigma) +{ + /* Generate gaussian noise centered around 128 */ + + Imaging imOut; + int x, y; + int nextok; + double this, next; + + imOut = ImagingNew("L", xsize, ysize); + if (!imOut) + return NULL; + + next = 0.0; + nextok = 0; + + for (y = 0; y < imOut->ysize; y++) { + UINT8* out = imOut->image8[y]; + for (x = 0; x < imOut->xsize; x++) { + if (nextok) { + this = next; + nextok = 0; + } else { + /* after numerical recepies */ + double v1, v2, radius, factor; + do { + v1 = rand()*(2.0/32767.0) - 1.0; + v2 = rand()*(2.0/32767.0) - 1.0; + radius= v1*v1 + v2*v2; + } while (radius >= 1.0); + factor = sqrt(-2.0*log(radius)/radius); + this = factor * v1; + next = factor * v2; + } + out[x] = (unsigned char) (128 + sigma * this); + } + } + + return imOut; +} + +Imaging +ImagingEffectPerlinTurbulence(int xsize, int ysize) +{ + /* Perlin turbulence (In progress) */ + + return NULL; +} + +Imaging +ImagingEffectSpread(Imaging imIn, int distance) +{ + /* Randomly spread pixels in an image */ + + Imaging imOut; + int x, y; + + imOut = ImagingNew(imIn->mode, imIn->xsize, imIn->ysize); + + if (!imOut) + return NULL; + +#define SPREAD(type, image)\ + for (y = 0; y < imIn->ysize; y++)\ + for (x = 0; x < imIn->xsize; x++) {\ + int xx = x + (rand() % distance) - distance/2;\ + int yy = y + (rand() % distance) - distance/2;\ + if (xx >= 0 && xx < imIn->xsize && yy >= 0 && yy < imIn->ysize) {\ + imOut->image[yy][xx] = imIn->image[y][x];\ + imOut->image[y][x] = imIn->image[yy][xx];\ + } else\ + imOut->image[y][x] = imIn->image[y][x];\ + } + + if (imIn->image8) { + SPREAD(UINT8, image8); + } else { + SPREAD(INT32, image32); + } + + ImagingCopyInfo(imOut, imIn); + + return imOut; +} + +/* -------------------------------------------------------------------- */ +/* Taken from the "C" code in the W3C SVG specification. Translated + to C89 by Fredrik Lundh */ + +/* Produces results in the range [1, 2**31 - 2]. +Algorithm is: r = (a * r) mod m +where a = 16807 and m = 2**31 - 1 = 2147483647 +See [Park & Miller], CACM vol. 31 no. 10 p. 1195, Oct. 1988 +To test: the algorithm should produce the result 1043618065 +as the 10,000th generated number if the original seed is 1. +*/ +#define RAND_m 2147483647 /* 2**31 - 1 */ +#define RAND_a 16807 /* 7**5; primitive root of m */ +#define RAND_q 127773 /* m / a */ +#define RAND_r 2836 /* m % a */ + +static long +perlin_setup_seed(long lSeed) +{ + if (lSeed <= 0) lSeed = -(lSeed % (RAND_m - 1)) + 1; + if (lSeed > RAND_m - 1) lSeed = RAND_m - 1; + return lSeed; +} + +static long +perlin_random(long lSeed) +{ + long result; + result = RAND_a * (lSeed % RAND_q) - RAND_r * (lSeed / RAND_q); + if (result <= 0) result += RAND_m; + return result; +} + +#define BSize 0x100 +#define BM 0xff +#define PerlinN 0x1000 +#define NP 12 /* 2^PerlinN */ +#define NM 0xfff +static int perlin_uLatticeSelector[BSize + BSize + 2]; +static double perlin_fGradient[4][BSize + BSize + 2][2]; +typedef struct +{ + int nWidth; /* How much to subtract to wrap for stitching. */ + int nHeight; + int nWrapX; /* Minimum value to wrap. */ + int nWrapY; +} StitchInfo; + +static void +perlin_init(long lSeed) +{ + double s; + int i, j, k; + lSeed = perlin_setup_seed(lSeed); + for(k = 0; k < 4; k++) + { + for(i = 0; i < BSize; i++) + { + perlin_uLatticeSelector[i] = i; + for (j = 0; j < 2; j++) + perlin_fGradient[k][i][j] = (double)(((lSeed = perlin_random(lSeed)) % (BSize + BSize)) - BSize) / BSize; + s = (double) (sqrt(perlin_fGradient[k][i][0] * perlin_fGradient[k][i][0] + perlin_fGradient[k][i][1] * perlin_fGradient[k][i][1])); + perlin_fGradient[k][i][0] /= s; + perlin_fGradient[k][i][1] /= s; + } + } + while(--i) + { + k = perlin_uLatticeSelector[i]; + perlin_uLatticeSelector[i] = perlin_uLatticeSelector[j = (lSeed = perlin_random(lSeed)) % BSize]; + perlin_uLatticeSelector[j] = k; + } + for(i = 0; i < BSize + 2; i++) + { + perlin_uLatticeSelector[BSize + i] = perlin_uLatticeSelector[i]; + for(k = 0; k < 4; k++) + for(j = 0; j < 2; j++) + perlin_fGradient[k][BSize + i][j] = perlin_fGradient[k][i][j]; + } +} + +#define s_curve(t) ( t * t * (3. - 2. * t) ) +#define lerp(t, a, b) ( a + t * (b - a) ) +static double +perlin_noise2(int nColorChannel, double vec[2], StitchInfo *pStitchInfo) +{ + int bx0, bx1, by0, by1, b00, b10, b01, b11; + double rx0, rx1, ry0, ry1, *q, sx, sy, a, b, t, u, v; + register int i, j; + + t = vec[0] + (double) PerlinN; + bx0 = (int)t; + bx1 = bx0+1; + rx0 = t - (int)t; + rx1 = rx0 - 1.0f; + t = vec[1] + (double) PerlinN; + by0 = (int)t; + by1 = by0+1; + ry0 = t - (int)t; + ry1 = ry0 - 1.0f; + + /* If stitching, adjust lattice points accordingly. */ + if(pStitchInfo != NULL) + { + if(bx0 >= pStitchInfo->nWrapX) + bx0 -= pStitchInfo->nWidth; + if(bx1 >= pStitchInfo->nWrapX) + bx1 -= pStitchInfo->nWidth; + if(by0 >= pStitchInfo->nWrapY) + by0 -= pStitchInfo->nHeight; + if(by1 >= pStitchInfo->nWrapY) + by1 -= pStitchInfo->nHeight; + } + + bx0 &= BM; + bx1 &= BM; + by0 &= BM; + by1 &= BM; + + i = perlin_uLatticeSelector[bx0]; + j = perlin_uLatticeSelector[bx1]; + b00 = perlin_uLatticeSelector[i + by0]; + b10 = perlin_uLatticeSelector[j + by0]; + b01 = perlin_uLatticeSelector[i + by1]; + b11 = perlin_uLatticeSelector[j + by1]; + sx = (double) (s_curve(rx0)); + sy = (double) (s_curve(ry0)); + q = perlin_fGradient[nColorChannel][b00]; u = rx0 * q[0] + ry0 * q[1]; + q = perlin_fGradient[nColorChannel][b10]; v = rx1 * q[0] + ry0 * q[1]; + a = lerp(sx, u, v); + q = perlin_fGradient[nColorChannel][b01]; u = rx0 * q[0] + ry1 * q[1]; + q = perlin_fGradient[nColorChannel][b11]; v = rx1 * q[0] + ry1 * q[1]; + b = lerp(sx, u, v); + return lerp(sy, a, b); +} + +double +perlin_turbulence( + int nColorChannel, double *point, double fBaseFreqX, double fBaseFreqY, + int nNumOctaves, int bFractalSum, int bDoStitching, + double fTileX, double fTileY, double fTileWidth, double fTileHeight) +{ + StitchInfo stitch; + StitchInfo *pStitchInfo = NULL; /* Not stitching when NULL. */ + + double fSum = 0.0f; + double vec[2]; + double ratio = 1; + + int nOctave; + + vec[0] = point[0] * fBaseFreqX; + vec[1] = point[1] * fBaseFreqY; + + /* Adjust the base frequencies if necessary for stitching. */ + if(bDoStitching) + { + /* When stitching tiled turbulence, the frequencies must be adjusted */ + /* so that the tile borders will be continuous. */ + if(fBaseFreqX != 0.0) + { + double fLoFreq = (double) (floor(fTileWidth * fBaseFreqX)) / fTileWidth; + double fHiFreq = (double) (ceil(fTileWidth * fBaseFreqX)) / fTileWidth; + if(fBaseFreqX / fLoFreq < fHiFreq / fBaseFreqX) + fBaseFreqX = fLoFreq; + else + fBaseFreqX = fHiFreq; + } + + if(fBaseFreqY != 0.0) + { + double fLoFreq = (double) (floor(fTileHeight * fBaseFreqY)) / fTileHeight; + double fHiFreq = (double) (ceil(fTileHeight * fBaseFreqY)) / fTileHeight; + if(fBaseFreqY / fLoFreq < fHiFreq / fBaseFreqY) + fBaseFreqY = fLoFreq; + else + fBaseFreqY = fHiFreq; + } + + /* Set up initial stitch values. */ + pStitchInfo = &stitch; + stitch.nWidth = (int) (fTileWidth * fBaseFreqX + 0.5f); + stitch.nWrapX = (int) (fTileX * fBaseFreqX + PerlinN + stitch.nWidth); + stitch.nHeight = (int) (fTileHeight * fBaseFreqY + 0.5f); + stitch.nWrapY = (int) (fTileY * fBaseFreqY + PerlinN + stitch.nHeight); + } + + for(nOctave = 0; nOctave < nNumOctaves; nOctave++) + { + if(bFractalSum) + fSum += (double) (perlin_noise2(nColorChannel, vec, pStitchInfo) / ratio); + else + fSum += (double) (fabs(perlin_noise2(nColorChannel, vec, pStitchInfo)) / ratio); + + vec[0] *= 2; + vec[1] *= 2; + ratio *= 2; + + if(pStitchInfo != NULL) + { + /* Update stitch values. Subtracting PerlinN before the multiplication and */ + /* adding it afterward simplifies to subtracting it once. */ + stitch.nWidth *= 2; + stitch.nWrapX = 2 * stitch.nWrapX - PerlinN; + stitch.nHeight *= 2; + stitch.nWrapY = 2 * stitch.nWrapY - PerlinN; + } + } + return fSum; +} + |