Merge branch 'gaussian-refactor' into fast-box-blur

This commit is contained in:
homm 2014-10-25 15:50:57 +04:00
commit d89c9ab750
7 changed files with 237 additions and 275 deletions

View File

@ -149,11 +149,12 @@ class GaussianBlur(Filter):
"""
name = "GaussianBlur"
def __init__(self, radius=2):
def __init__(self, radius=2, effective_scale=None):
self.radius = radius
self.effective_scale = effective_scale
def filter(self, image):
return image.gaussian_blur(self.radius)
return image.gaussian_blur(self.radius, self.effective_scale or 2.6)
class UnsharpMask(Filter):

View File

@ -414,15 +414,18 @@ def solarize(image, threshold=128):
# --------------------------------------------------------------------
# PIL USM components, from Kevin Cazabon.
def gaussian_blur(im, radius=None):
""" PIL_usm.gblur(im, [radius])"""
def gaussian_blur(im, radius=None, effective_scale=None):
""" PIL_usm.gblur(im, [radius], [effective_scale])"""
if radius is None:
radius = 5.0
if effective_scale is None:
effective_scale = 2.6
im.load()
return im.im.gaussian_blur(radius)
return im.im.gaussian_blur(radius, effective_scale)
gblur = gaussian_blur

Binary file not shown.

After

Width:  |  Height:  |  Size: 1.3 KiB

View File

@ -5,6 +5,7 @@ from PIL import ImageOps
from PIL import ImageFilter
im = Image.open("Tests/images/hopper.ppm")
snakes = Image.open("Tests/images/color_snakes.png")
class TestImageOpsUsm(PillowTestCase):
@ -16,7 +17,7 @@ class TestImageOpsUsm(PillowTestCase):
self.assertEqual(i.size, (128, 128))
# i.save("blur.bmp")
i = ImageOps.usm(im, 2.0, 125, 8)
i = ImageOps.unsharp_mask(im, 2.0, 125, 8)
self.assertEqual(i.mode, "RGB")
self.assertEqual(i.size, (128, 128))
# i.save("usm.bmp")
@ -33,7 +34,7 @@ class TestImageOpsUsm(PillowTestCase):
self.assertEqual(i.mode, "RGB")
self.assertEqual(i.size, (128, 128))
def test_usm(self):
def test_usm_formats(self):
usm = ImageOps.unsharp_mask
self.assertRaises(ValueError, lambda: usm(im.convert("1")))
@ -45,7 +46,7 @@ class TestImageOpsUsm(PillowTestCase):
usm(im.convert("CMYK"))
self.assertRaises(ValueError, lambda: usm(im.convert("YCbCr")))
def test_blur(self):
def test_blur_formats(self):
blur = ImageOps.gaussian_blur
self.assertRaises(ValueError, lambda: blur(im.convert("1")))
@ -57,6 +58,33 @@ class TestImageOpsUsm(PillowTestCase):
blur(im.convert("CMYK"))
self.assertRaises(ValueError, lambda: blur(im.convert("YCbCr")))
def test_usm_accuracy(self):
i = snakes._new(ImageOps.unsharp_mask(snakes, 5, 1024, 0))
# Image should not be changed because it have only 0 and 255 levels.
self.assertEqual(i.tobytes(), snakes.tobytes())
def test_blur_accuracy(self):
i = snakes._new(ImageOps.gaussian_blur(snakes, .7))
# Alpha channel must match whole.
self.assertEqual(i.split()[3], snakes.split()[3])
# These pixels surrounded with pixels with 255 intensity.
# They must be very close to 255.
for x, y, c in [(1, 0, 1), (2, 0, 1), (7, 8, 1), (8, 8, 1), (2, 9, 1),
(7, 3, 0), (8, 3, 0), (5, 8, 0), (5, 9, 0), (1, 3, 0),
(4, 3, 2), (4, 2, 2)]:
self.assertGreaterEqual(i.im.getpixel((x, y))[c], 250)
# Fuzzy match.
gp = lambda x, y: i.im.getpixel((x, y))
self.assertTrue(211 <= gp(7, 4)[0] <= 213)
self.assertTrue(211 <= gp(7, 5)[2] <= 213)
self.assertTrue(211 <= gp(7, 6)[2] <= 213)
self.assertTrue(211 <= gp(7, 7)[1] <= 213)
self.assertTrue(211 <= gp(8, 4)[0] <= 213)
self.assertTrue(211 <= gp(8, 5)[2] <= 213)
self.assertTrue(211 <= gp(8, 6)[2] <= 213)
self.assertTrue(211 <= gp(8, 7)[1] <= 213)
if __name__ == '__main__':
unittest.main()

View File

@ -863,7 +863,8 @@ _gaussian_blur(ImagingObject* self, PyObject* args)
Imaging imOut;
float radius = 0;
if (!PyArg_ParseTuple(args, "f", &radius))
float effectiveScale = 2.6;
if (!PyArg_ParseTuple(args, "f|f", &radius, &effectiveScale))
return NULL;
imIn = self->image;
@ -871,7 +872,7 @@ _gaussian_blur(ImagingObject* self, PyObject* args)
if (!imOut)
return NULL;
if (!ImagingGaussianBlur(imIn, imOut, radius))
if (!ImagingGaussianBlur(imIn, imOut, radius, effectiveScale))
return NULL;
return PyImagingNew(imOut);

View File

@ -263,7 +263,8 @@ extern Imaging ImagingFilter(
FLOAT32 offset, FLOAT32 divisor);
extern Imaging ImagingFlipLeftRight(Imaging imOut, Imaging imIn);
extern Imaging ImagingFlipTopBottom(Imaging imOut, Imaging imIn);
extern Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius);
extern Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius,
float effectiveScale);
extern Imaging ImagingGetBand(Imaging im, int band);
extern int ImagingGetBBox(Imaging im, int bbox[4]);
typedef struct { int x, y; INT32 count; INT32 pixel; } ImagingColorItem;

View File

@ -9,62 +9,16 @@
#include "Python.h"
#include "Imaging.h"
#define PILUSMVERSION "0.6.1"
/* version history
0.6.1 converted to C and added to PIL 1.1.7
0.6.0 fixed/improved float radius support (oops!)
now that radius can be a float (properly), changed radius value to
be an actual radius (instead of diameter). So, you should get
similar results from PIL_usm as from other paint programs when
using the SAME values (no doubling of radius required any more).
Be careful, this may "break" software if you had it set for 2x
or 5x the radius as was recommended with earlier versions.
made PILusm thread-friendly (release GIL before lengthly operations,
and re-acquire it before returning to Python). This makes a huge
difference with multi-threaded applications on dual-processor
or "Hyperthreading"-enabled systems (Pentium4, Xeon, etc.)
0.5.0 added support for float radius values!
0.4.0 tweaked gaussian curve calculation to be closer to consistent shape
across a wide range of radius values
0.3.0 changed deviation calculation in gausian algorithm to be dynamic
_gblur now adds 1 to the user-supplied radius before using it so
that a value of "0" returns the original image instead of a
black one.
fixed handling of alpha channel in RGBX, RGBA images
improved speed of gblur by reducing unnecessary checks and assignments
0.2.0 fixed L-mode image support
0.1.0 initial release
*/
static inline UINT8 clip(double in)
{
if (in >= 255.0)
return (UINT8) 255;
if (in <= 0.0)
return (UINT8) 0;
return (UINT8) in;
}
static Imaging
gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
gblur(Imaging im, Imaging imOut, float radius, float effectiveScale, int channels)
{
ImagingSectionCookie cookie;
float *maskData = NULL;
int y = 0;
int x = 0;
float z = 0;
float sum = 0.0;
float dev = 0.0;
float *buffer = NULL;
@ -75,12 +29,10 @@ gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
float newPixel[4];
int channel = 0;
int offset = 0;
INT32 newPixelFinals;
int radius = 0;
float remainder = 0.0;
int i;
int effectiveRadius = 0;
int window = 0;
int hasAlpha = 0;
/* Do the gaussian blur */
@ -91,156 +43,130 @@ gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
radius of 5 instead of 25 lookups). So, we blur the lines first,
then we blur the resulting columns. */
/* first, round radius off to the next higher integer and hold the
remainder this is used so we can support float radius values
properly. */
remainder = floatRadius - ((int) floatRadius);
floatRadius = ceil(floatRadius);
/* Next, double the radius and offset by 2.0... that way "0" returns
the original image instead of a black one. We multiply it by 2.0
so that it is a true "radius", not a diameter (the results match
other paint programs closer that way too). */
radius = (int) ((floatRadius * 2.0) + 2.0);
/* Only pixels in effective radius from source pixel are accounted.
The Gaussian values outside 3 x radius is near zero. */
effectiveRadius = (int) ceil(radius * effectiveScale);
/* Window is number of pixels forming the result pixel on one axis.
It is source pixel and effective radius in both directions. */
window = effectiveRadius * 2 + 1;
/* create the maskData for the gaussian curve */
maskData = malloc(radius * sizeof(float));
/* FIXME: error checking */
for (x = 0; x < radius; x++) {
z = ((float) (x + 2) / ((float) radius));
dev = 0.5 + (((float) (radius * radius)) * 0.001);
/* you can adjust this factor to change the shape/center-weighting
of the gaussian */
maskData[x] = (float) pow((1.0 / sqrt(2.0 * 3.14159265359 * dev)),
((-(z - 1.0) * -(x - 1.0)) /
(2.0 * dev)));
maskData = malloc(window * sizeof(float));
for (pix = 0; pix < window; pix++) {
offset = pix - effectiveRadius;
if (radius) {
/* http://en.wikipedia.org/wiki/Gaussian_blur
"1 / sqrt(2 * pi * dev)" is constant and will be eliminated
by normalization. */
maskData[pix] = pow(2.718281828459,
-offset * offset / (2 * radius * radius));
} else {
maskData[pix] = 1;
}
sum += maskData[pix];
}
/* if there's any remainder, multiply the first/last values in
MaskData it. this allows us to support float radius values. */
if (remainder > 0.0) {
maskData[0] *= remainder;
maskData[radius - 1] *= remainder;
}
for (x = 0; x < radius; x++) {
/* this is done separately now due to the correction for float
radius values above */
sum += maskData[x];
}
for (i = 0; i < radius; i++) {
maskData[i] *= (1.0 / sum);
/* printf("%f\n", maskData[i]); */
for (pix = 0; pix < window; pix++) {
maskData[pix] *= (1.0 / sum);
// printf("%d %f\n", pix, maskData[pix]);
}
// printf("\n");
/* create a temporary memory buffer for the data for the first pass
memset the buffer to 0 so we can use it directly with += */
/* don't bother about alpha/padding */
/* don't bother about alpha */
buffer = calloc((size_t) (im->xsize * im->ysize * channels),
sizeof(float));
sizeof(float));
if (buffer == NULL)
return ImagingError_MemoryError();
return ImagingError_MemoryError();
/* be nice to other threads while you go off to lala land */
ImagingSectionEnter(&cookie);
/* memset(buffer, 0, sizeof(buffer)); */
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
/* perform a blur on each line, and place in the temporary storage buffer */
for (y = 0; y < im->ysize; y++) {
if (channels == 1 && im->image8 != NULL) {
line8 = (UINT8 *) im->image8[y];
} else {
line = im->image32[y];
}
for (x = 0; x < im->xsize; x++) {
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
/* for each neighbor pixel, factor in its value/weighting to the
current pixel */
for (pix = 0; pix < radius; pix++) {
/* figure the offset of this neighbor pixel */
offset =
(int) ((-((float) radius / 2.0) + (float) pix) + 0.5);
if (x + offset < 0)
offset = -x;
else if (x + offset >= im->xsize)
offset = im->xsize - x - 1;
if (channels == 1 && im->image8 != NULL) {
line8 = (UINT8 *) im->image8[y];
} else {
line = im->image32[y];
}
for (x = 0; x < im->xsize; x++) {
/* for each neighbor pixel, factor in its value/weighting to the
current pixel */
for (pix = 0; pix < window; pix++) {
/* figure the offset of this neighbor pixel */
offset = pix - effectiveRadius;
if (x + offset < 0)
offset = -x;
else if (x + offset >= im->xsize)
offset = im->xsize - x - 1;
/* add (neighbor pixel value * maskData[pix]) to the current
pixel value */
if (channels == 1) {
buffer[(y * im->xsize) + x] +=
((float) ((UINT8 *) & line8[x + offset])[0]) *
(maskData[pix]);
} else {
for (channel = 0; channel < channels; channel++) {
buffer[(y * im->xsize * channels) +
(x * channels) + channel] +=
((float) ((UINT8 *) & line[x + offset])
[channel]) * (maskData[pix]);
}
}
}
}
/* add (neighbor pixel value * maskData[pix]) to the current
pixel value */
if (channels == 1) {
buffer[(y * im->xsize) + x] +=
((float) ((UINT8 *) & line8[x + offset])[0]) *
(maskData[pix]);
} else {
for (channel = 0; channel < channels; channel++) {
buffer[(y * im->xsize * channels) +
(x * channels) + channel] +=
((float) ((UINT8 *) & line[x + offset])
[channel]) * (maskData[pix]);
}
}
}
}
}
if (strcmp(im->mode, "RGBX") == 0 || strcmp(im->mode, "RGBA") == 0) {
hasAlpha = 1;
}
/* perform a blur on each column in the buffer, and place in the
output image */
for (x = 0; x < im->xsize; x++) {
for (y = 0; y < im->ysize; y++) {
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = 0;
/* for each neighbor pixel, factor in its value/weighting to the
current pixel */
for (pix = 0; pix < radius; pix++) {
/* figure the offset of this neighbor pixel */
offset =
(int) (-((float) radius / 2.0) + (float) pix + 0.5);
if (y + offset < 0)
offset = -y;
else if (y + offset >= im->ysize)
offset = im->ysize - y - 1;
/* add (neighbor pixel value * maskData[pix]) to the current
pixel value */
for (channel = 0; channel < channels; channel++) {
newPixel[channel] +=
(buffer
[((y + offset) * im->xsize * channels) +
(x * channels) + channel]) * (maskData[pix]);
}
}
/* if the image is RGBX or RGBA, copy the 4th channel data to
newPixel, so it gets put in imOut */
if (strcmp(im->mode, "RGBX") == 0
|| strcmp(im->mode, "RGBA") == 0) {
newPixel[3] = (float) ((UINT8 *) & line[x + offset])[3];
}
for (y = 0; y < im->ysize; y++) {
newPixel[0] = newPixel[1] = newPixel[2] = newPixel[3] = .5;
/* for each neighbor pixel, factor in its value/weighting to the
current pixel */
for (pix = 0; pix < window; pix++) {
/* figure the offset of this neighbor pixel */
offset = pix - effectiveRadius;
if (y + offset < 0)
offset = -y;
else if (y + offset >= im->ysize)
offset = im->ysize - y - 1;
/* pack the channels into an INT32 so we can put them back in
the PIL image */
newPixelFinals = 0;
if (channels == 1) {
newPixelFinals = clip(newPixel[0]);
} else {
/* for RGB, the fourth channel isn't used anyways, so just
pack a 0 in there, this saves checking the mode for each
pixel. */
/* this doesn't work on little-endian machines... fix it! */
newPixelFinals =
clip(newPixel[0]) | clip(newPixel[1]) << 8 |
clip(newPixel[2]) << 16 | clip(newPixel[3]) << 24;
}
/* set the resulting pixel in imOut */
if (channels == 1) {
imOut->image8[y][x] = (UINT8) newPixelFinals;
} else {
imOut->image32[y][x] = newPixelFinals;
}
}
/* add (neighbor pixel value * maskData[pix]) to the current
pixel value */
for (channel = 0; channel < channels; channel++) {
newPixel[channel] +=
(buffer
[((y + offset) * im->xsize * channels) +
(x * channels) + channel]) * (maskData[pix]);
}
}
if (channels == 1) {
imOut->image8[y][x] = (UINT8)(newPixel[0]);
} else {
/* if the image is RGBX or RGBA, copy the 4th channel data to
newPixel, so it gets put in imOut */
if (hasAlpha) {
newPixel[3] = (float) ((UINT8 *) & im->image32[y][x])[3];
}
/* for RGB, the fourth channel isn't used anyways, so just
pack a 0 in there, this saves checking the mode for each
pixel. */
/* this might don't work on little-endian machines... fix it! */
imOut->image32[y][x] =
(UINT8)(newPixel[0]) | (UINT8)(newPixel[1]) << 8 |
(UINT8)(newPixel[2]) << 16 | (UINT8)(newPixel[3]) << 24;
}
}
}
/* free the buffer */
@ -252,42 +178,46 @@ gblur(Imaging im, Imaging imOut, float floatRadius, int channels, int padding)
return imOut;
}
Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius)
static inline UINT8 clip(double in)
{
if (in >= 255.0)
return (UINT8) 255;
if (in <= 0.0)
return (UINT8) 0;
return (UINT8) (in + 0.5);
}
Imaging ImagingGaussianBlur(Imaging im, Imaging imOut, float radius,
float effectiveScale)
{
int channels = 0;
int padding = 0;
if (strcmp(im->mode, "RGB") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "RGBA") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "RGBX") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "CMYK") == 0) {
channels = 4;
padding = 0;
channels = 4;
} else if (strcmp(im->mode, "L") == 0) {
channels = 1;
padding = 0;
channels = 1;
} else
return ImagingError_ModeError();
return ImagingError_ModeError();
return gblur(im, imOut, radius, channels, padding);
return gblur(im, imOut, radius, effectiveScale, channels);
}
Imaging
ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
int threshold)
int threshold)
{
ImagingSectionCookie cookie;
Imaging result;
int channel = 0;
int channels = 0;
int padding = 0;
int hasAlpha = 0;
int x = 0;
int y = 0;
@ -302,28 +232,23 @@ ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
INT32 newPixel = 0;
if (strcmp(im->mode, "RGB") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "RGBA") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "RGBX") == 0) {
channels = 3;
padding = 1;
channels = 3;
} else if (strcmp(im->mode, "CMYK") == 0) {
channels = 4;
padding = 0;
channels = 4;
} else if (strcmp(im->mode, "L") == 0) {
channels = 1;
padding = 0;
channels = 1;
} else
return ImagingError_ModeError();
return ImagingError_ModeError();
/* first, do a gaussian blur on the image, putting results in imOut
temporarily */
result = gblur(im, imOut, radius, channels, padding);
result = gblur(im, imOut, radius, 2.6, channels);
if (!result)
return NULL;
return NULL;
/* now, go through each pixel, compare "normal" pixel to blurred
pixel. if the difference is more than threshold values, apply
@ -332,64 +257,67 @@ ImagingUnsharpMask(Imaging im, Imaging imOut, float radius, int percent,
ImagingSectionEnter(&cookie);
for (y = 0; y < im->ysize; y++) {
if (channels == 1) {
lineIn8 = im->image8[y];
lineOut8 = imOut->image8[y];
} else {
lineIn = im->image32[y];
lineOut = imOut->image32[y];
}
for (x = 0; x < im->xsize; x++) {
newPixel = 0;
/* compare in/out pixels, apply sharpening */
if (channels == 1) {
diff =
((UINT8 *) & lineIn8[x])[0] -
((UINT8 *) & lineOut8[x])[0];
if (abs(diff) > threshold) {
/* add the diff*percent to the original pixel */
imOut->image8[y][x] =
clip((((UINT8 *) & lineIn8[x])[0]) +
(diff * ((float) percent) / 100.0));
} else {
/* newPixel is the same as imIn */
imOut->image8[y][x] = ((UINT8 *) & lineIn8[x])[0];
}
}
if (strcmp(im->mode, "RGBX") == 0 || strcmp(im->mode, "RGBA") == 0) {
hasAlpha = 1;
}
else {
for (channel = 0; channel < channels; channel++) {
diff = (int) ((((UINT8 *) & lineIn[x])[channel]) -
(((UINT8 *) & lineOut[x])[channel]));
if (abs(diff) > threshold) {
/* add the diff*percent to the original pixel
this may not work for little-endian systems, fix it! */
newPixel =
newPixel |
clip((float) (((UINT8 *) & lineIn[x])[channel])
+
(diff *
(((float) percent /
100.0)))) << (channel * 8);
} else {
/* newPixel is the same as imIn
this may not work for little-endian systems, fix it! */
newPixel =
newPixel | ((UINT8 *) & lineIn[x])[channel] <<
(channel * 8);
}
}
if (strcmp(im->mode, "RGBX") == 0
|| strcmp(im->mode, "RGBA") == 0) {
/* preserve the alpha channel
this may not work for little-endian systems, fix it! */
newPixel =
newPixel | ((UINT8 *) & lineIn[x])[channel] << 24;
}
imOut->image32[y][x] = newPixel;
}
}
for (y = 0; y < im->ysize; y++) {
if (channels == 1) {
lineIn8 = im->image8[y];
lineOut8 = imOut->image8[y];
} else {
lineIn = im->image32[y];
lineOut = imOut->image32[y];
}
for (x = 0; x < im->xsize; x++) {
newPixel = 0;
/* compare in/out pixels, apply sharpening */
if (channels == 1) {
diff =
((UINT8 *) & lineIn8[x])[0] -
((UINT8 *) & lineOut8[x])[0];
if (abs(diff) > threshold) {
/* add the diff*percent to the original pixel */
imOut->image8[y][x] =
clip((((UINT8 *) & lineIn8[x])[0]) +
(diff * ((float) percent) / 100.0));
} else {
/* newPixel is the same as imIn */
imOut->image8[y][x] = ((UINT8 *) & lineIn8[x])[0];
}
}
else {
for (channel = 0; channel < channels; channel++) {
diff = (int) ((((UINT8 *) & lineIn[x])[channel]) -
(((UINT8 *) & lineOut[x])[channel]));
if (abs(diff) > threshold) {
/* add the diff*percent to the original pixel
this may not work for little-endian systems, fix it! */
newPixel =
newPixel |
clip((float) (((UINT8 *) & lineIn[x])[channel])
+
(diff *
(((float) percent /
100.0)))) << (channel * 8);
} else {
/* newPixel is the same as imIn
this may not work for little-endian systems, fix it! */
newPixel =
newPixel | ((UINT8 *) & lineIn[x])[channel] <<
(channel * 8);
}
}
if (hasAlpha) {
/* preserve the alpha channel
this may not work for little-endian systems, fix it! */
newPixel =
newPixel | ((UINT8 *) & lineIn[x])[channel] << 24;
}
imOut->image32[y][x] = newPixel;
}
}
}
ImagingSectionLeave(&cookie);