diff --git a/Tests/test_image_filter.py b/Tests/test_image_filter.py index 3a970520ccb..a097881844e 100644 --- a/Tests/test_image_filter.py +++ b/Tests/test_image_filter.py @@ -192,6 +192,24 @@ def test_rankfilter_overflow() -> None: im.filter(rankfilter) +@pytest.mark.parametrize( + "mode, fill", + ( + ("L", 0x12), + ("I", 0x12345678), + ("I;16", 0x1234), + ("I;16B", 0x1234), + ("F", 1.5), + ("RGB", (0x12, 0x34, 0x56)), + ), +) +@pytest.mark.parametrize("margin", (0, 1, 3)) +def test_expand(mode: str, fill: float | tuple[int, ...], margin: int) -> None: + im = Image.new(mode, (5, 4), fill) + out = im._new(im.im.expand(margin)) # As done internally by RankFilter + assert_image_equal(out, Image.new(mode, (5 + 2 * margin, 4 + 2 * margin), fill)) + + def test_builtinfilter_p() -> None: builtin_filter = ImageFilter.BuiltinFilter() diff --git a/src/_imaging.c b/src/_imaging.c index cc782698245..9ace9a6cd8b 100644 --- a/src/_imaging.c +++ b/src/_imaging.c @@ -1756,7 +1756,7 @@ _putdata(ImagingObject *self, PyObject *args) { for (i = x = y = 0; i < n; i++) { double value; set_value_to_item(seq, i); - IMAGING_PIXEL_INT32(image, x, y) = (INT32)(value * scale + offset); + image->image32[y][x] = (INT32)(value * scale + offset); if (++x >= (int)image->xsize) { x = 0, y++; } @@ -1766,7 +1766,7 @@ _putdata(ImagingObject *self, PyObject *args) { for (i = x = y = 0; i < n; i++) { double value; set_value_to_item(seq, i); - IMAGING_PIXEL_FLOAT32(image, x, y) = + ((FLOAT32 *)image->image32[y])[x] = (FLOAT32)(value * scale + offset); if (++x >= (int)image->xsize) { x = 0, y++; diff --git a/src/libImaging/Filter.c b/src/libImaging/Filter.c index 7982419ff16..619c0a83be8 100644 --- a/src/libImaging/Filter.c +++ b/src/libImaging/Filter.c @@ -76,38 +76,43 @@ ImagingExpand(Imaging imIn, int margin) { return NULL; } -#define EXPAND_LINE(type, image, yin, yout) \ - { \ - for (x = 0; x < margin; x++) { \ - imOut->image[yout][x] = imIn->image[yin][0]; \ - } \ - for (x = 0; x < imIn->xsize; x++) { \ - imOut->image[yout][x + margin] = imIn->image[yin][x]; \ - } \ - for (x = 0; x < margin; x++) { \ - imOut->image[yout][margin + imIn->xsize + x] = \ - imIn->image[yin][imIn->xsize - 1]; \ - } \ + int xsize = imIn->xsize, ysize = imIn->ysize; + +#define EXPAND_LINE(type, yin, yout) \ + { \ + const type *in = (const type *)imIn->image[yin]; \ + type *out = (type *)imOut->image[yout]; \ + for (x = 0; x < margin; x++) { \ + out[x] = in[0]; \ + } \ + for (x = 0; x < xsize; x++) { \ + out[x + margin] = in[x]; \ + } \ + for (x = 0; x < margin; x++) { \ + out[margin + xsize + x] = in[xsize - 1]; \ + } \ } -#define EXPAND(type, image) \ - { \ - for (y = 0; y < margin; y++) { \ - EXPAND_LINE(type, image, 0, y); \ - } \ - for (y = 0; y < imIn->ysize; y++) { \ - EXPAND_LINE(type, image, y, y + margin); \ - } \ - for (y = 0; y < margin; y++) { \ - EXPAND_LINE(type, image, imIn->ysize - 1, margin + imIn->ysize + y); \ - } \ +#define EXPAND(type) \ + { \ + for (y = 0; y < margin; y++) { \ + EXPAND_LINE(type, 0, y); \ + } \ + for (y = 0; y < ysize; y++) { \ + EXPAND_LINE(type, y, y + margin); \ + } \ + for (y = 0; y < margin; y++) { \ + EXPAND_LINE(type, ysize - 1, margin + ysize + y); \ + } \ } ImagingSectionEnter(&cookie); - if (imIn->image8) { - EXPAND(UINT8, image8); + if (imIn->type == IMAGING_TYPE_I16) { + EXPAND(UINT16); + } else if (imIn->image8) { + EXPAND(UINT8); } else { - EXPAND(INT32, image32); + EXPAND(INT32); } ImagingSectionLeave(&cookie); diff --git a/src/libImaging/Imaging.h b/src/libImaging/Imaging.h index 20f74b519a3..47196654f32 100644 --- a/src/libImaging/Imaging.h +++ b/src/libImaging/Imaging.h @@ -132,10 +132,6 @@ struct ImagingMemoryInstance { #define IMAGING_PIXEL_CMYK(im, x, y) ((im)->image[(y)][(x) * 4]) #define IMAGING_PIXEL_YCbCr(im, x, y) ((im)->image[(y)][(x) * 4]) -#define IMAGING_PIXEL_UINT8(im, x, y) ((im)->image8[(y)][(x)]) -#define IMAGING_PIXEL_INT32(im, x, y) ((im)->image32[(y)][(x)]) -#define IMAGING_PIXEL_FLOAT32(im, x, y) (((FLOAT32 *)(im)->image32[y])[x]) - struct ImagingAccessInstance { ModeID mode; void (*get_pixel)(Imaging im, int x, int y, void *pixel); diff --git a/src/libImaging/RankFilter.c b/src/libImaging/RankFilter.c index fc8ce2db96d..d8aedbd9225 100644 --- a/src/libImaging/RankFilter.c +++ b/src/libImaging/RankFilter.c @@ -17,52 +17,232 @@ /* Fast rank algorithm (due to Wirth), based on public domain code by Nicolas Devillard, available at http://ndevilla.free.fr */ -#define SWAP(type, a, b) \ - { \ - register type t = (a); \ - (a) = (b); \ - (b) = t; \ - } +#define RANK_INNER_BODY(type) \ + int i, j, l, m; \ + type x; \ + l = 0; \ + m = n - 1; \ + while (l < m) { \ + x = a[k]; \ + i = l; \ + j = m; \ + do { \ + while (a[i] < x) { \ + i++; \ + } \ + while (x < a[j]) { \ + j--; \ + } \ + if (i <= j) { \ + type t = a[i]; \ + a[i] = a[j]; \ + a[j] = t; \ + i++; \ + j--; \ + } \ + } while (i <= j); \ + if (j < k) { \ + l = i; \ + } \ + if (k < i) { \ + m = j; \ + } \ + } \ + return a[k] -#define MakeRankFunction(type) \ - static type Rank##type(type a[], int n, int k) { \ - register int i, j, l, m; \ - register type x; \ - l = 0; \ - m = n - 1; \ - while (l < m) { \ - x = a[k]; \ - i = l; \ - j = m; \ - do { \ - while (a[i] < x) { \ - i++; \ - } \ - while (x < a[j]) { \ - j--; \ - } \ - if (i <= j) { \ - SWAP(type, a[i], a[j]); \ - i++; \ - j--; \ - } \ - } while (i <= j); \ - if (j < k) { \ - l = i; \ - } \ - if (k < i) { \ - m = j; \ - } \ - } \ - return a[k]; \ - } +static UINT8 +RankUINT8(UINT8 a[], int n, int k) { + RANK_INNER_BODY(UINT8); +} + +static INT32 +RankINT32(INT32 a[], int n, int k) { + RANK_INNER_BODY(INT32); +} + +static FLOAT32 +RankFLOAT32(FLOAT32 a[], int n, int k) { + RANK_INNER_BODY(FLOAT32); +} + +#define RANK_MIN(a, b) ((a) < (b) ? (a) : (b)) +#define RANK_MAX(a, b) ((a) < (b) ? (b) : (a)) + +#define MINMAX_BODY(type, op) \ + do { \ + for (int y = 0; y < ysize; y++) { \ + type *out = (type *)imOut->image[y]; \ + memcpy(out, im->image[y], xsize * sizeof(type)); \ + for (int i = 0; i < size; i++) { \ + const type *row = (const type *)im->image[y + i]; \ + for (int k = (i == 0); k < size; k++) { \ + const type *in = row + k; \ + for (int x = 0; x < xsize; x++) { \ + out[x] = op(out[x], in[x]); \ + } \ + } \ + } \ + } \ + } while (0) + +// Median selection networks for 3x3 and 5x5 windows, +// original public-domain code via http://ndevilla.free.fr/median/median/src/optmed.c + +#define CSWAP(type, p, i, j) \ + do { \ + const type lo_ = p[i] < p[j] ? p[i] : p[j]; \ + p[j] = p[i] < p[j] ? p[j] : p[i]; \ + p[i] = lo_; \ + } while (0) + +#define MEDIAN_NETWORK_3(type, p) \ + CSWAP(type, p, 1, 2); \ + CSWAP(type, p, 4, 5); \ + CSWAP(type, p, 7, 8); \ + CSWAP(type, p, 0, 1); \ + CSWAP(type, p, 3, 4); \ + CSWAP(type, p, 6, 7); \ + CSWAP(type, p, 1, 2); \ + CSWAP(type, p, 4, 5); \ + CSWAP(type, p, 7, 8); \ + CSWAP(type, p, 0, 3); \ + CSWAP(type, p, 5, 8); \ + CSWAP(type, p, 4, 7); \ + CSWAP(type, p, 3, 6); \ + CSWAP(type, p, 1, 4); \ + CSWAP(type, p, 2, 5); \ + CSWAP(type, p, 4, 7); \ + CSWAP(type, p, 4, 2); \ + CSWAP(type, p, 6, 4); \ + CSWAP(type, p, 4, 2) -MakeRankFunction(UINT8) MakeRankFunction(INT32) MakeRankFunction(FLOAT32) +#define MEDIAN_NETWORK_5(type, p) \ + CSWAP(type, p, 0, 1); \ + CSWAP(type, p, 3, 4); \ + CSWAP(type, p, 2, 4); \ + CSWAP(type, p, 2, 3); \ + CSWAP(type, p, 6, 7); \ + CSWAP(type, p, 5, 7); \ + CSWAP(type, p, 5, 6); \ + CSWAP(type, p, 9, 10); \ + CSWAP(type, p, 8, 10); \ + CSWAP(type, p, 8, 9); \ + CSWAP(type, p, 12, 13); \ + CSWAP(type, p, 11, 13); \ + CSWAP(type, p, 11, 12); \ + CSWAP(type, p, 15, 16); \ + CSWAP(type, p, 14, 16); \ + CSWAP(type, p, 14, 15); \ + CSWAP(type, p, 18, 19); \ + CSWAP(type, p, 17, 19); \ + CSWAP(type, p, 17, 18); \ + CSWAP(type, p, 21, 22); \ + CSWAP(type, p, 20, 22); \ + CSWAP(type, p, 20, 21); \ + CSWAP(type, p, 23, 24); \ + CSWAP(type, p, 2, 5); \ + CSWAP(type, p, 3, 6); \ + CSWAP(type, p, 0, 6); \ + CSWAP(type, p, 0, 3); \ + CSWAP(type, p, 4, 7); \ + CSWAP(type, p, 1, 7); \ + CSWAP(type, p, 1, 4); \ + CSWAP(type, p, 11, 14); \ + CSWAP(type, p, 8, 14); \ + CSWAP(type, p, 8, 11); \ + CSWAP(type, p, 12, 15); \ + CSWAP(type, p, 9, 15); \ + CSWAP(type, p, 9, 12); \ + CSWAP(type, p, 13, 16); \ + CSWAP(type, p, 10, 16); \ + CSWAP(type, p, 10, 13); \ + CSWAP(type, p, 20, 23); \ + CSWAP(type, p, 17, 23); \ + CSWAP(type, p, 17, 20); \ + CSWAP(type, p, 21, 24); \ + CSWAP(type, p, 18, 24); \ + CSWAP(type, p, 18, 21); \ + CSWAP(type, p, 19, 22); \ + CSWAP(type, p, 8, 17); \ + CSWAP(type, p, 9, 18); \ + CSWAP(type, p, 0, 18); \ + CSWAP(type, p, 0, 9); \ + CSWAP(type, p, 10, 19); \ + CSWAP(type, p, 1, 19); \ + CSWAP(type, p, 1, 10); \ + CSWAP(type, p, 11, 20); \ + CSWAP(type, p, 2, 20); \ + CSWAP(type, p, 2, 11); \ + CSWAP(type, p, 12, 21); \ + CSWAP(type, p, 3, 21); \ + CSWAP(type, p, 3, 12); \ + CSWAP(type, p, 13, 22); \ + CSWAP(type, p, 4, 22); \ + CSWAP(type, p, 4, 13); \ + CSWAP(type, p, 14, 23); \ + CSWAP(type, p, 5, 23); \ + CSWAP(type, p, 5, 14); \ + CSWAP(type, p, 15, 24); \ + CSWAP(type, p, 6, 24); \ + CSWAP(type, p, 6, 15); \ + CSWAP(type, p, 7, 16); \ + CSWAP(type, p, 7, 19); \ + CSWAP(type, p, 13, 21); \ + CSWAP(type, p, 15, 23); \ + CSWAP(type, p, 7, 13); \ + CSWAP(type, p, 7, 15); \ + CSWAP(type, p, 1, 9); \ + CSWAP(type, p, 3, 11); \ + CSWAP(type, p, 5, 17); \ + CSWAP(type, p, 11, 17); \ + CSWAP(type, p, 9, 17); \ + CSWAP(type, p, 4, 10); \ + CSWAP(type, p, 6, 12); \ + CSWAP(type, p, 7, 14); \ + CSWAP(type, p, 4, 6); \ + CSWAP(type, p, 4, 7); \ + CSWAP(type, p, 12, 14); \ + CSWAP(type, p, 10, 14); \ + CSWAP(type, p, 6, 7); \ + CSWAP(type, p, 10, 12); \ + CSWAP(type, p, 6, 10); \ + CSWAP(type, p, 6, 17); \ + CSWAP(type, p, 12, 17); \ + CSWAP(type, p, 7, 17); \ + CSWAP(type, p, 7, 10); \ + CSWAP(type, p, 12, 18); \ + CSWAP(type, p, 7, 12); \ + CSWAP(type, p, 10, 18); \ + CSWAP(type, p, 12, 20); \ + CSWAP(type, p, 10, 20); \ + CSWAP(type, p, 10, 12) - Imaging ImagingRankFilter(Imaging im, int size, int rank) { +// restrict safe: imOut is a fresh allocation. +#define MEDIAN_BODY(type, size) \ + do { \ + for (int y = 0; y < ysize; y++) { \ + const type *rows[size]; \ + for (int i = 0; i < size; i++) { \ + rows[i] = (const type *)im->image[y + i]; \ + } \ + type *restrict out = (type *)imOut->image[y]; \ + for (int x = 0; x < xsize; x++) { \ + type p[size * size]; \ + for (int i = 0; i < size; i++) { \ + for (int k = 0; k < size; k++) { \ + p[i * size + k] = rows[i][x + k]; \ + } \ + } \ + MEDIAN_NETWORK_##size(type, p); \ + out[x] = p[size * size / 2]; \ + } \ + } \ + } while (0) + +Imaging +ImagingRankFilter(Imaging im, int size, int rank) { Imaging imOut = NULL; - int x, y; - int i, margin, size2; + int margin, size2; if (!im || im->bands != 1 || im->type == IMAGING_TYPE_I16) { return (Imaging)ImagingError_ModeError(); @@ -89,35 +269,53 @@ MakeRankFunction(UINT8) MakeRankFunction(INT32) MakeRankFunction(FLOAT32) if (!imOut) { return NULL; } + int xsize = imOut->xsize, ysize = imOut->ysize; /* malloc check ok, checked above */ -#define RANK_BODY(type) \ - do { \ - type *buf = malloc(size2 * sizeof(type)); \ - if (!buf) { \ - goto nomemory; \ - } \ - for (y = 0; y < imOut->ysize; y++) { \ - for (x = 0; x < imOut->xsize; x++) { \ - for (i = 0; i < size; i++) { \ - memcpy( \ - buf + i * size, \ - &IMAGING_PIXEL_##type(im, x, y + i), \ - size * sizeof(type) \ - ); \ - } \ - IMAGING_PIXEL_##type(imOut, x, y) = Rank##type(buf, size2, rank); \ - } \ - } \ - free(buf); \ + // restrict safe: buf is a private allocation, imOut is a fresh allocation. +#define RANK_BODY(type, rank_fn) \ + do { \ + type *restrict buf = malloc(size2 * sizeof(type)); \ + if (!buf) { \ + goto nomemory; \ + } \ + for (int y = 0; y < ysize; y++) { \ + type *restrict out = (type *)imOut->image[y]; \ + for (int x = 0; x < xsize; x++) { \ + type *restrict p = buf; \ + for (int i = 0; i < size; i++) { \ + const type *row = (const type *)im->image[y + i] + x; \ + for (int k = 0; k < size; k++) { \ + *p++ = row[k]; \ + } \ + } \ + out[x] = rank_fn(buf, size2, rank); \ + } \ + } \ + free(buf); \ + } while (0) + +#define RANK_DISPATCH(type) \ + do { \ + if (rank == 0) { \ + MINMAX_BODY(type, RANK_MIN); \ + } else if (rank == size2 - 1) { \ + MINMAX_BODY(type, RANK_MAX); \ + } else if (rank == size2 / 2 && size == 3) { \ + MEDIAN_BODY(type, 3); \ + } else if (rank == size2 / 2 && size == 5) { \ + MEDIAN_BODY(type, 5); \ + } else { \ + RANK_BODY(type, Rank##type); \ + } \ } while (0) if (im->image8) { - RANK_BODY(UINT8); + RANK_DISPATCH(UINT8); } else if (im->type == IMAGING_TYPE_INT32) { - RANK_BODY(INT32); + RANK_DISPATCH(INT32); } else if (im->type == IMAGING_TYPE_FLOAT32) { - RANK_BODY(FLOAT32); + RANK_DISPATCH(FLOAT32); } else { /* safety net (we shouldn't end up here) */ ImagingDelete(imOut);