1
/* PILusm, a gaussian blur and unsharp masking library for PIL
2
By Kevin Cazabon, copyright 2003
3
kevin_cazabon@hotmail.com
6
/* Originally released under LGPL. Graciously donated to PIL
7
for distribution under the standard PIL license in 2009." */
12
#define PILUSMVERSION "0.6.1"
16
0.6.1 converted to C and added to PIL 1.1.7
18
0.6.0 fixed/improved float radius support (oops!)
19
now that radius can be a float (properly), changed radius value to
20
be an actual radius (instead of diameter). So, you should get
21
similar results from PIL_usm as from other paint programs when
22
using the SAME values (no doubling of radius required any more).
23
Be careful, this may "break" software if you had it set for 2x
24
or 5x the radius as was recommended with earlier versions.
25
made PILusm thread-friendly (release GIL before lengthly operations,
26
and re-acquire it before returning to Python). This makes a huge
27
difference with multi-threaded applications on dual-processor
28
or "Hyperthreading"-enabled systems (Pentium4, Xeon, etc.)
30
0.5.0 added support for float radius values!
32
0.4.0 tweaked gaussian curve calculation to be closer to consistent shape
33
across a wide range of radius values
35
0.3.0 changed deviation calculation in gausian algorithm to be dynamic
36
_gblur now adds 1 to the user-supplied radius before using it so
37
that a value of "0" returns the original image instead of a
39
fixed handling of alpha channel in RGBX, RGBA images
40
improved speed of gblur by reducing unnecessary checks and assignments
42
0.2.0 fixed L-mode image support
48
static inline UINT8 clip(double in)
58
gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
60
ImagingSectionCookie cookie;
62
float *maskData = NULL;
81
float remainder = 0.0;
85
/* Do the gaussian blur */
87
/* For a symmetrical gaussian blur, instead of doing a radius*radius
88
matrix lookup, you get the EXACT same results by doing a radius*1
89
transform, followed by a 1*radius transform. This reduces the
90
number of lookups exponentially (10 lookups per pixel for a
91
radius of 5 instead of 25 lookups). So, we blur the lines first,
92
then we blur the resulting columns. */
94
/* first, round radius off to the next higher integer and hold the
95
remainder this is used so we can support float radius values
98
remainder = floatRadius - ((int) floatRadius);
99
floatRadius = ceil(floatRadius);
101
/* Next, double the radius and offset by 2.0... that way "0" returns
102
the original image instead of a black one. We multiply it by 2.0
103
so that it is a true "radius", not a diameter (the results match
104
other paint programs closer that way too). */
105
radius = (int) ((floatRadius * 2.0) + 2.0);
107
/* create the maskData for the gaussian curve */
108
maskData = malloc(radius * sizeof(float));
109
/* FIXME: error checking */
110
for (x = 0; x < radius; x++) {
111
z = ((float) (x + 2) / ((float) radius));
112
dev = 0.5 + (((float) (radius * radius)) * 0.001);
113
/* you can adjust this factor to change the shape/center-weighting
115
maskData[x] = (float) pow((1.0 / sqrt(2.0 * 3.14159265359 * dev)),
116
((-(z - 1.0) * -(x - 1.0)) /
120
/* if there's any remainder, multiply the first/last values in
121
MaskData it. this allows us to support float radius values. */
122
if (remainder > 0.0) {
123
maskData[0] *= remainder;
124
maskData[radius - 1] *= remainder;
127
for (x = 0; x < radius; x++) {
128
/* this is done separately now due to the correction for float
129
radius values above */
133
for (i = 0; i < radius; i++) {
134
maskData[i] *= (1.0 / sum);
135
/* printf("%f\n", maskData[i]); */
138
/* create a temporary memory buffer for the data for the first pass
139
memset the buffer to 0 so we can use it directly with += */
141
/* don't bother about alpha/padding */
142
buffer = calloc((size_t) (im->xsize * im->ysize * channels),
145
return ImagingError_MemoryError();
147
/* be nice to other threads while you go off to lala land */
148
ImagingSectionEnter(&cookie);
150
/* memset(buffer, 0, sizeof(buffer)); */
152
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
154
/* perform a blur on each line, and place in the temporary storage buffer */
155
for (y = 0; y < im->ysize; y++) {
156
if (channels == 1 && im->image8 != NULL) {
157
line8 = (UINT8 *) im->image8[y];
159
line = im->image32[y];
161
for (x = 0; x < im->xsize; x++) {
162
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
163
/* for each neighbor pixel, factor in its value/weighting to the
165
for (pix = 0; pix < radius; pix++) {
166
/* figure the offset of this neighbor pixel */
168
(int) ((-((float) radius / 2.0) + (float) pix) + 0.5);
171
else if (x + offset >= im->xsize)
172
offset = im->xsize - x - 1;
174
/* add (neighbor pixel value * maskData[pix]) to the current
177
buffer[(y * im->xsize) + x] +=
178
((float) ((UINT8 *) & line8[x + offset])[0]) *
181
for (channel = 0; channel < channels; channel++) {
182
buffer[(y * im->xsize * channels) +
183
(x * channels) + channel] +=
184
((float) ((UINT8 *) & line[x + offset])
185
[channel]) * (maskData[pix]);
192
/* perform a blur on each column in the buffer, and place in the
194
for (x = 0; x < im->xsize; x++) {
195
for (y = 0; y < im->ysize; y++) {
196
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
197
/* for each neighbor pixel, factor in its value/weighting to the
199
for (pix = 0; pix < radius; pix++) {
200
/* figure the offset of this neighbor pixel */
202
(int) (-((float) radius / 2.0) + (float) pix + 0.5);
205
else if (y + offset >= im->ysize)
206
offset = im->ysize - y - 1;
207
/* add (neighbor pixel value * maskData[pix]) to the current
209
for (channel = 0; channel < channels; channel++) {
212
[((y + offset) * im->xsize * channels) +
213
(x * channels) + channel]) * (maskData[pix]);
216
/* if the image is RGBX or RGBA, copy the 4th channel data to
217
newPixel, so it gets put in imOut */
218
if (strcmp(im->mode, "RGBX") == 0
219
|| strcmp(im->mode, "RGBA") == 0) {
220
newPixel[3] = (float) ((UINT8 *) & line[x + offset])[3];
223
/* pack the channels into an INT32 so we can put them back in
227
newPixelFinals = clip(newPixel[0]);
229
/* for RGB, the fourth channel isn't used anyways, so just
230
pack a 0 in there, this saves checking the mode for each
232
/* this doesn't work on little-endian machines... fix it! */
234
clip(newPixel[0]) | clip(newPixel[1]) << 8 |
235
clip(newPixel[2]) << 16 | clip(newPixel[3]) << 24;
237
/* set the resulting pixel in imOut */
239
imOut->image8[y][x] = (UINT8) newPixelFinals;
241
imOut->image32[y][x] = newPixelFinals;
246
/* free the buffer */
249
/* get the GIL back so Python knows who you are */
250
ImagingSectionLeave(&cookie);
255
Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius)
260
if (strcmp(im->mode, "RGB") == 0) {
263
} else if (strcmp(im->mode, "RGBA") == 0) {
266
} else if (strcmp(im->mode, "RGBX") == 0) {
269
} else if (strcmp(im->mode, "CMYK") == 0) {
272
} else if (strcmp(im->mode, "L") == 0) {
276
return ImagingError_ModeError();
278
return gblur(im, imOut, radius, channels, padding);
282
ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
285
ImagingSectionCookie cookie;
297
UINT8 *lineIn8 = NULL;
298
UINT8 *lineOut8 = NULL;
304
if (strcmp(im->mode, "RGB") == 0) {
307
} else if (strcmp(im->mode, "RGBA") == 0) {
310
} else if (strcmp(im->mode, "RGBX") == 0) {
313
} else if (strcmp(im->mode, "CMYK") == 0) {
316
} else if (strcmp(im->mode, "L") == 0) {
320
return ImagingError_ModeError();
322
/* first, do a gaussian blur on the image, putting results in imOut
324
result = gblur(im, imOut, radius, channels, padding);
328
/* now, go through each pixel, compare "normal" pixel to blurred
329
pixel. if the difference is more than threshold values, apply
330
the OPPOSITE correction to the amount of blur, multiplied by
333
ImagingSectionEnter(&cookie);
335
for (y = 0; y < im->ysize; y++) {
337
lineIn8 = im->image8[y];
338
lineOut8 = imOut->image8[y];
340
lineIn = im->image32[y];
341
lineOut = imOut->image32[y];
343
for (x = 0; x < im->xsize; x++) {
345
/* compare in/out pixels, apply sharpening */
348
((UINT8 *) & lineIn8[x])[0] -
349
((UINT8 *) & lineOut8[x])[0];
350
if (abs(diff) > threshold) {
351
/* add the diff*percent to the original pixel */
352
imOut->image8[y][x] =
353
clip((((UINT8 *) & lineIn8[x])[0]) +
354
(diff * ((float) percent) / 100.0));
356
/* newPixel is the same as imIn */
357
imOut->image8[y][x] = ((UINT8 *) & lineIn8[x])[0];
362
for (channel = 0; channel < channels; channel++) {
363
diff = (int) ((((UINT8 *) & lineIn[x])[channel]) -
364
(((UINT8 *) & lineOut[x])[channel]));
365
if (abs(diff) > threshold) {
366
/* add the diff*percent to the original pixel
367
this may not work for little-endian systems, fix it! */
370
clip((float) (((UINT8 *) & lineIn[x])[channel])
374
100.0)))) << (channel * 8);
376
/* newPixel is the same as imIn
377
this may not work for little-endian systems, fix it! */
379
newPixel | ((UINT8 *) & lineIn[x])[channel] <<
383
if (strcmp(im->mode, "RGBX") == 0
384
|| strcmp(im->mode, "RGBA") == 0) {
385
/* preserve the alpha channel
386
this may not work for little-endian systems, fix it! */
388
newPixel | ((UINT8 *) & lineIn[x])[channel] << 24;
390
imOut->image32[y][x] = newPixel;
395
ImagingSectionLeave(&cookie);