~ubuntu-branches/debian/squeeze/python-imaging/squeeze

« back to all changes in this revision

Viewing changes to libImaging/UnsharpMask.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2009-11-20 19:22:59 UTC
  • mfrom: (1.1.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20091120192259-n3iy0f17n5akogom
Tags: 1.1.7-1
New upstream version.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/* PILusm, a gaussian blur and unsharp masking library for PIL
 
2
   By Kevin Cazabon, copyright 2003
 
3
   kevin_cazabon@hotmail.com
 
4
   kevin@cazabon.com */
 
5
 
 
6
/* Originally released under LGPL.  Graciously donated to PIL
 
7
   for distribution under the standard PIL license in 2009." */
 
8
 
 
9
#include "Python.h"
 
10
#include "Imaging.h"
 
11
 
 
12
#define PILUSMVERSION "0.6.1"
 
13
 
 
14
/* version history
 
15
 
 
16
0.6.1   converted to C and added to PIL 1.1.7
 
17
 
 
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.)
 
29
 
 
30
0.5.0   added support for float radius values!
 
31
 
 
32
0.4.0   tweaked gaussian curve calculation to be closer to consistent shape
 
33
            across a wide range of radius values
 
34
 
 
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
 
38
            black one.
 
39
        fixed handling of alpha channel in RGBX, RGBA images
 
40
        improved speed of gblur by reducing unnecessary checks and assignments
 
41
 
 
42
0.2.0   fixed L-mode image support
 
43
 
 
44
0.1.0   initial release
 
45
 
 
46
*/
 
47
 
 
48
static inline UINT8 clip(double in)
 
49
{
 
50
    if (in >= 255.0)
 
51
        return (UINT8) 255;
 
52
    if (in <= 0.0)
 
53
        return (UINT8) 0;
 
54
    return (UINT8) in;
 
55
}
 
56
 
 
57
static Imaging
 
58
gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
 
59
{
 
60
    ImagingSectionCookie cookie;
 
61
 
 
62
    float *maskData = NULL;
 
63
    int y = 0;
 
64
    int x = 0;
 
65
    float z = 0;
 
66
    float sum = 0.0;
 
67
    float dev = 0.0;
 
68
 
 
69
    float *buffer = NULL;
 
70
 
 
71
    int *line = NULL;
 
72
    UINT8 *line8 = NULL;
 
73
 
 
74
    int pix = 0;
 
75
    float newPixel[4];
 
76
    int channel = 0;
 
77
    int offset = 0;
 
78
    INT32 newPixelFinals;
 
79
 
 
80
    int radius = 0;
 
81
    float remainder = 0.0;
 
82
 
 
83
    int i;
 
84
 
 
85
    /* Do the gaussian blur */
 
86
 
 
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. */
 
93
 
 
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
 
96
       properly. */
 
97
 
 
98
    remainder = floatRadius - ((int) floatRadius);
 
99
    floatRadius = ceil(floatRadius);
 
100
 
 
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);
 
106
 
 
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
 
114
           of the gaussian */
 
115
        maskData[x] = (float) pow((1.0 / sqrt(2.0 * 3.14159265359 * dev)),
 
116
                                  ((-(z - 1.0) * -(x - 1.0)) /
 
117
                                   (2.0 * dev)));
 
118
    }
 
119
 
 
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;
 
125
    }
 
126
 
 
127
    for (x = 0; x < radius; x++) {
 
128
        /* this is done separately now due to the correction for float
 
129
           radius values above */
 
130
        sum += maskData[x];
 
131
    }
 
132
 
 
133
    for (i = 0; i < radius; i++) {
 
134
        maskData[i] *= (1.0 / sum);
 
135
        /* printf("%f\n", maskData[i]); */
 
136
    }
 
137
 
 
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 += */
 
140
 
 
141
    /* don't bother about alpha/padding */
 
142
    buffer = calloc((size_t) (im->xsize * im->ysize * channels),
 
143
                    sizeof(float));
 
144
    if (buffer == NULL)
 
145
        return ImagingError_MemoryError();
 
146
 
 
147
    /* be nice to other threads while you go off to lala land */
 
148
    ImagingSectionEnter(&cookie);
 
149
 
 
150
    /* memset(buffer, 0, sizeof(buffer)); */
 
151
 
 
152
    newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
 
153
 
 
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];
 
158
        } else {
 
159
            line = im->image32[y];
 
160
        }
 
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
 
164
               current pixel */
 
165
            for (pix = 0; pix < radius; pix++) {
 
166
                /* figure the offset of this neighbor pixel */
 
167
                offset =
 
168
                    (int) ((-((float) radius / 2.0) + (float) pix) + 0.5);
 
169
                if (x + offset < 0)
 
170
                    offset = -x;
 
171
                else if (x + offset >= im->xsize)
 
172
                    offset = im->xsize - x - 1;
 
173
 
 
174
                /* add (neighbor pixel value * maskData[pix]) to the current
 
175
                   pixel value */
 
176
                if (channels == 1) {
 
177
                    buffer[(y * im->xsize) + x] +=
 
178
                        ((float) ((UINT8 *) & line8[x + offset])[0]) *
 
179
                        (maskData[pix]);
 
180
                } else {
 
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]);
 
186
                    }
 
187
                }
 
188
            }
 
189
        }
 
190
    }
 
191
 
 
192
    /* perform a blur on each column in the buffer, and place in the
 
193
       output image */
 
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
 
198
               current pixel */
 
199
            for (pix = 0; pix < radius; pix++) {
 
200
                /* figure the offset of this neighbor pixel */
 
201
                offset =
 
202
                    (int) (-((float) radius / 2.0) + (float) pix + 0.5);
 
203
                if (y + offset < 0)
 
204
                    offset = -y;
 
205
                else if (y + offset >= im->ysize)
 
206
                    offset = im->ysize - y - 1;
 
207
                /* add (neighbor pixel value * maskData[pix]) to the current
 
208
                   pixel value */
 
209
                for (channel = 0; channel < channels; channel++) {
 
210
                    newPixel[channel] +=
 
211
                        (buffer
 
212
                         [((y + offset) * im->xsize * channels) +
 
213
                          (x * channels) + channel]) * (maskData[pix]);
 
214
                }
 
215
            }
 
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];
 
221
            }
 
222
 
 
223
            /* pack the channels into an INT32 so we can put them back in
 
224
               the PIL image */
 
225
            newPixelFinals = 0;
 
226
            if (channels == 1) {
 
227
                newPixelFinals = clip(newPixel[0]);
 
228
            } else {
 
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
 
231
                   pixel. */
 
232
                /* this doesn't work on little-endian machines... fix it! */
 
233
                newPixelFinals =
 
234
                    clip(newPixel[0]) | clip(newPixel[1]) << 8 |
 
235
                    clip(newPixel[2]) << 16 | clip(newPixel[3]) << 24;
 
236
            }
 
237
            /* set the resulting pixel in imOut */
 
238
            if (channels == 1) {
 
239
                imOut->image8[y][x] = (UINT8) newPixelFinals;
 
240
            } else {
 
241
                imOut->image32[y][x] = newPixelFinals;
 
242
            }
 
243
        }
 
244
    }
 
245
 
 
246
    /* free the buffer */
 
247
    free(buffer);
 
248
 
 
249
    /* get the GIL back so Python knows who you are */
 
250
    ImagingSectionLeave(&cookie);
 
251
 
 
252
    return imOut;
 
253
}
 
254
 
 
255
Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius)
 
256
{
 
257
    int channels = 0;
 
258
    int padding = 0;
 
259
 
 
260
    if (strcmp(im->mode, "RGB") == 0) {
 
261
        channels = 3;
 
262
        padding = 1;
 
263
    } else if (strcmp(im->mode, "RGBA") == 0) {
 
264
        channels = 3;
 
265
        padding = 1;
 
266
    } else if (strcmp(im->mode, "RGBX") == 0) {
 
267
        channels = 3;
 
268
        padding = 1;
 
269
    } else if (strcmp(im->mode, "CMYK") == 0) {
 
270
        channels = 4;
 
271
        padding = 0;
 
272
    } else if (strcmp(im->mode, "L") == 0) {
 
273
        channels = 1;
 
274
        padding = 0;
 
275
    } else
 
276
        return ImagingError_ModeError();
 
277
 
 
278
    return gblur(im, imOut, radius, channels, padding);
 
279
}
 
280
 
 
281
Imaging
 
282
ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
 
283
                   int threshold)
 
284
{
 
285
    ImagingSectionCookie cookie;
 
286
 
 
287
    Imaging result;
 
288
    int channel = 0;
 
289
    int channels = 0;
 
290
    int padding = 0;
 
291
 
 
292
    int x = 0;
 
293
    int y = 0;
 
294
 
 
295
    int *lineIn = NULL;
 
296
    int *lineOut = NULL;
 
297
    UINT8 *lineIn8 = NULL;
 
298
    UINT8 *lineOut8 = NULL;
 
299
 
 
300
    int diff = 0;
 
301
 
 
302
    INT32 newPixel = 0;
 
303
 
 
304
    if (strcmp(im->mode, "RGB") == 0) {
 
305
        channels = 3;
 
306
        padding = 1;
 
307
    } else if (strcmp(im->mode, "RGBA") == 0) {
 
308
        channels = 3;
 
309
        padding = 1;
 
310
    } else if (strcmp(im->mode, "RGBX") == 0) {
 
311
        channels = 3;
 
312
        padding = 1;
 
313
    } else if (strcmp(im->mode, "CMYK") == 0) {
 
314
        channels = 4;
 
315
        padding = 0;
 
316
    } else if (strcmp(im->mode, "L") == 0) {
 
317
        channels = 1;
 
318
        padding = 0;
 
319
    } else
 
320
        return ImagingError_ModeError();
 
321
 
 
322
    /* first, do a gaussian blur on the image, putting results in imOut
 
323
       temporarily */
 
324
    result = gblur(im, imOut, radius, channels, padding);
 
325
    if (!result)
 
326
        return NULL;
 
327
 
 
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
 
331
       percent. */
 
332
 
 
333
    ImagingSectionEnter(&cookie);
 
334
 
 
335
    for (y = 0; y < im->ysize; y++) {
 
336
        if (channels == 1) {
 
337
            lineIn8 = im->image8[y];
 
338
            lineOut8 = imOut->image8[y];
 
339
        } else {
 
340
            lineIn = im->image32[y];
 
341
            lineOut = imOut->image32[y];
 
342
        }
 
343
        for (x = 0; x < im->xsize; x++) {
 
344
            newPixel = 0;
 
345
            /* compare in/out pixels, apply sharpening */
 
346
            if (channels == 1) {
 
347
                diff =
 
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));
 
355
                } else {
 
356
                    /* newPixel is the same as imIn */
 
357
                    imOut->image8[y][x] = ((UINT8 *) & lineIn8[x])[0];
 
358
                }
 
359
            }
 
360
 
 
361
            else {
 
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! */
 
368
                        newPixel =
 
369
                            newPixel |
 
370
                            clip((float) (((UINT8 *) & lineIn[x])[channel])
 
371
                                 +
 
372
                                 (diff *
 
373
                                  (((float) percent /
 
374
                                    100.0)))) << (channel * 8);
 
375
                    } else {
 
376
                        /* newPixel is the same as imIn
 
377
                           this may not work for little-endian systems, fix it! */
 
378
                        newPixel =
 
379
                            newPixel | ((UINT8 *) & lineIn[x])[channel] <<
 
380
                            (channel * 8);
 
381
                    }
 
382
                }
 
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! */
 
387
                    newPixel =
 
388
                        newPixel | ((UINT8 *) & lineIn[x])[channel] << 24;
 
389
                }
 
390
                imOut->image32[y][x] = newPixel;
 
391
            }
 
392
        }
 
393
    }
 
394
 
 
395
    ImagingSectionLeave(&cookie);
 
396
 
 
397
    return imOut;
 
398
}