source: libpipi/trunk/pipi/filter/median.c

Last change on this file was 3564, checked in by Sam Hocevar, 10 years ago

That optimisation sucked. Reverted median filter to something that works
both on Linux and Windows.

File size: 4.0 KB
Line 
1/*
2 *  libpipi       Pathetic image processing interface library
3 *  Copyright (c) 2004-2008 Sam Hocevar <sam@zoy.org>
4 *                All Rights Reserved
5 *
6 *  $Id$
7 *
8 *  This library is free software. It comes without any warranty, to
9 *  the extent permitted by applicable law. You can redistribute it
10 *  and/or modify it under the terms of the Do What The Fuck You Want
11 *  To Public License, Version 2, as published by Sam Hocevar. See
12 *  http://sam.zoy.org/wtfpl/COPYING for more details.
13 */
14
15/*
16 * median.c: median filter functions
17 */
18
19#include "config.h"
20
21#include <stdlib.h>
22#include <stdio.h>
23#include <string.h>
24#include <math.h>
25
26#include "pipi.h"
27#include "pipi_internals.h"
28
29static int cmpdouble(void const *i1, void const *i2)
30{
31    double a = *(double const *)i1;
32    double b = *(double const *)i2;
33
34    return (a > b) - (a < b);
35}
36
37pipi_image_t *pipi_median(pipi_image_t *src, int radius)
38{
39    return pipi_median_ext(src, radius, radius);
40}
41
42/* FIXME: this is in desperate want of optimisation. Here is what could
43 * be done to improve the performance:
44 *  - prefetch the neighbourhood; most neighbours are the same as the
45 *    previous pixels.
46 *  - use a better sort algorithm; bubble sort is ridiculous
47 *  - even better, use state-of-the art median selection, ie. O(3n), for
48 *    most common combinations (9, 25, 49, 81). */
49pipi_image_t *pipi_median_ext(pipi_image_t *src, int rx, int ry)
50{
51    pipi_image_t *dst;
52    pipi_pixels_t *srcp, *dstp;
53    float *srcdata, *dstdata;
54    double *list;
55    int x, y, w, h, i, j, size, gray;
56
57    w = src->w;
58    h = src->h;
59    size = (2 * rx + 1) * (2 * ry + 1);
60
61    gray = (src->last_modified == PIPI_PIXELS_Y_F32);
62
63    srcp = gray ? pipi_get_pixels(src, PIPI_PIXELS_Y_F32)
64                : pipi_get_pixels(src, PIPI_PIXELS_RGBA_F32);
65    srcdata = (float *)srcp->pixels;
66
67    dst = pipi_new(w, h);
68    dstp = gray ? pipi_get_pixels(dst, PIPI_PIXELS_Y_F32)
69                : pipi_get_pixels(dst, PIPI_PIXELS_RGBA_F32);
70    dstdata = (float *)dstp->pixels;
71
72    list = malloc(size * (gray ? 1 : 4) * sizeof(double));
73
74    for(y = 0; y < h; y++)
75    {
76        for(x = 0; x < w; x++)
77        {
78            double *parser = list;
79
80            /* Make a list of neighbours */
81            for(j = -ry; j <= ry; j++)
82            {
83                int y2 = y + j;
84                if(y2 < 0) y2 = h - 1 - ((-y2 - 1) % h);
85                else if(y2 > 0) y2 = y2 % h;
86
87                for(i = -rx; i <= rx; i++)
88                {
89                    int x2 = x + i;
90                    if(x2 < 0) x2 = w - 1 - ((-x2 - 1) % w);
91                    else if(x2 > 0) x2 = x2 % w;
92
93                    if(gray)
94                    {
95                        *parser++ = srcdata[y2 * w + x2];
96                    }
97                    else
98                    {
99                        parser[0] = srcdata[4 * (y2 * w + x2)];
100                        parser[size * 1] = srcdata[4 * (y2 * w + x2) + 1];
101                        parser[size * 2] = srcdata[4 * (y2 * w + x2) + 2];
102                        parser[size * 3] = srcdata[4 * (y2 * w + x2) + 3];
103                        parser++;
104                    }
105                }
106            }
107
108            /* Sort the list */
109            qsort(list, size, sizeof(double), cmpdouble);
110            if(!gray)
111            {
112                qsort(list + size, size, sizeof(double), cmpdouble);
113                qsort(list + 2 * size, size, sizeof(double), cmpdouble);
114                qsort(list + 3 * size, size, sizeof(double), cmpdouble);
115            }
116
117            /* Store the median value */
118            if(gray)
119            {
120                dstdata[y * w + x] = list[size / 2];
121            }
122            else
123            {
124                dstdata[4 * (y * w + x)] = list[size / 2];
125                dstdata[4 * (y * w + x) + 1] = list[size / 2 + size * 1];
126                dstdata[4 * (y * w + x) + 2] = list[size / 2 + size * 2];
127                dstdata[4 * (y * w + x) + 3] = list[size / 2 + size * 3];
128            }
129        }
130    }
131
132    return dst;
133}
134
Note: See TracBrowser for help on using the repository browser.