Pillow/libImaging/RankFilter.c

/*
 * The Python Imaging Library
 * $Id$
 *
 * min, max, median filters
 *
 * history:
 * 2002-06-08 fl    Created
 *
 * Copyright (c) Secret Labs AB 2002.  All rights reserved.
 *
 * See the README file for information on usage and redistribution.
 */

#include "Imaging.h"

/* 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 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];\
}

MakeRankFunction(UINT8)
MakeRankFunction(INT32)
MakeRankFunction(FLOAT32)

Imaging
ImagingRankFilter(Imaging im, int size, int rank)
{
    Imaging imOut = NULL;
    int x, y;
    int i, margin, size2;

    if (!im || im->bands != 1 || im->type == IMAGING_TYPE_SPECIAL)
        return (Imaging) ImagingError_ModeError();

    if (!(size & 1))
        return (Imaging) ImagingError_ValueError("bad filter size");
    
    /* malloc check ok, for overflow in the define below */
    if (size > INT_MAX / size || 
        size > INT_MAX / (size * sizeof(FLOAT32))) {
        return (Imaging) ImagingError_ValueError("filter size too large");
    }

    size2 = size * size;
    margin = (size-1) / 2;

    if (rank < 0 || rank >= size2)
        return (Imaging) ImagingError_ValueError("bad rank value");

    imOut = ImagingNew(im->mode, im->xsize - 2*margin, im->ysize - 2*margin);
    if (!imOut)
        return NULL;

    /* 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); \
} while (0)

    if (im->image8)
        RANK_BODY(UINT8);
    else if (im->type == IMAGING_TYPE_INT32)
        RANK_BODY(INT32);
    else if (im->type == IMAGING_TYPE_FLOAT32)
        RANK_BODY(FLOAT32);
    else {
        /* safety net (we shouldn't end up here) */
        ImagingDelete(imOut);
        return (Imaging) ImagingError_ModeError();
    }

    ImagingCopyInfo(imOut, im);

    return imOut;

nomemory:
    ImagingDelete(imOut);
    return (Imaging) ImagingError_MemoryError();
}