5

I'm working on an algorithm that manipulates pictures. Basically I will implement a diffusion (each pixel will get the median value of the 8 surrounding pixels + its own value).

what I will do is to create a array of 9 integers with the value, sort the array and get the median value at array[4].

I still don't know what to use for the problem, what is the best sorting function to use for relatively small arrays? The sorting function will roughly be called x times, x being the number of pixels.

Heapsort seems a bit overkill. Quicksort will not perform that well. And I don't want to implement really complex things.

What do you guys think?

alain.janinm
  • 19,951
  • 10
  • 65
  • 112
Sword22
  • 266
  • 3
  • 15

3 Answers3

17

If all you need is the median, there's no need to do any sorting at all! (For long arrays, see http://en.wikipedia.org/wiki/Selection_algorithm for an O(n) algorithm; of course we're talking only about short arrays here).

For median of 9 numbers, a little googling reveals the article Fast median search: an ANSI C implementation by N. Devillard, which points to the article Implementing median filters in XC4000E FPGAs by J. L. Smith, which provides this self-explanatory "sorting network" to get the median using 19 comparisons:

enter image description here

In terms of C:

typedef int T;

void sort2(T* a, T* b);
void sort3(T* a, T* b, T* c);
T min3(T a, T b, T c);
T max3(T a, T b, T c);

T median9(T p1, T p2, T p3, T p4, T p5, T p6, T p7, T p8, T p9)
{
    sort3(&p1, &p2, &p3);
    sort3(&p4, &p5, &p6);
    sort3(&p7, &p8, &p9);

    p7 = max3(p1, p4, p7);
    p3 = min3(p3, p6, p9);

    sort3(&p2, &p5, &p8);
    sort3(&p3, &p5, &p7);

    return p5;
}

void sort2(T* a, T* b)
{
    if (*a > *b)
    {
        T tmp = *b;
        *b = *a;
        *a = tmp;
    }
}

void sort3(T* a, T* b, T* c)
{
    sort2(b, c);
    sort2(a, b);
    sort2(b, c);
}

T min3(T a, T b, T c)
{
    if (a < b)
        return a < c ? a : c;
    else
        return b < c ? b : c;
}

T max3(T a, T b, T c)
{
    if (a > b)
        return a > c ? a : c;
    else
        return b > c ? b : c;
}

Edit: this file also contains the code for getting the median of 3, 5, 6, 7, 9 and 25 numbers.

#define PIX_SORT(a,b) { if ((a)>(b)) PIX_SWAP((a),(b)); }
#define PIX_SWAP(a,b) { pixelvalue temp=(a);(a)=(b);(b)=temp; }

/*----------------------------------------------------------------------------
   Function :   opt_med9()
   In       :   pointer to an array of 9 pixelvalues
   Out      :   a pixelvalue
   Job      :   optimized search of the median of 9 pixelvalues
   Notice   :   in theory, cannot go faster without assumptions on the
                signal.
                Formula from:
                XILINX XCELL magazine, vol. 23 by John L. Smith

                The input array is modified in the process
                The result array is guaranteed to contain the median
                value
                in middle position, but other elements are NOT sorted.
 ---------------------------------------------------------------------------*/

pixelvalue opt_med9(pixelvalue * p)
{
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
    PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ;
    PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
    PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ;
    PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ;
    PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ;
    PIX_SORT(p[4], p[2]) ; return(p[4]) ;
}
kennytm
  • 510,854
  • 105
  • 1,084
  • 1,005
  • 3
    +1. Someone finally pointed out that you don't need a complete sort for this. Plus some very clean code. – Oliver Charlesworth Mar 15 '12 at 14:39
  • 1
    using a sorting network (as proposed by Jeff Foster's link), you'd end up with 25 comparisons, so ~30% overhead... – Christoph Mar 15 '12 at 14:42
  • indeed sorting networks are meant for n large, were overhead can be considered small in comparison to the total running time. However, this was overlooked. –  Mar 15 '12 at 14:58
  • 1
    Performance point: the last two invocations of `sort3` (within the `median9`) may be replaced by something that finds the median, rather than sorting all 3 elements – valdo Mar 15 '12 at 15:35
  • @valdo, nice: `T med3(a,b,c){ return a – AShelly Mar 15 '12 at 21:45
2

You know that the "standard" qsort of the c library is normally quite optimized? Often it's quicksort + insertion sort for small sub-arrays (when you subdivide very much the base array). For example read the comments of this version taken from the gnu c library

xanatos
  • 109,618
  • 12
  • 197
  • 280
  • 2
    `qsort` will probably be appalling for this, for two reasons. (1) 9 is a very small number. (2) You will incur the function call overhead for every single comparison operation (this is why C++'s `std::sort` is substantially faster than `qsort`). – Oliver Charlesworth Mar 15 '12 at 13:54
  • @OliCharlesworth When I have those problems I normally take the source from a BSD licensed version and replace the parts of the code that are useless to me :-) (so I would replace all the calls to the comparator with direct compares) – xanatos Mar 15 '12 at 13:57
  • The only reason that one has to use qsort over merge sort is space complexity, O(logn) vs O(n) for merge sort. If you are to use a library implementation use merge sort which is O(nlogn) time complexity, which make it a good choice for small arrays. –  Mar 15 '12 at 14:08
  • @g24l Very simply: unless I have VERY special problems, the single line I need to use `qsort` (the c library one) is much better than maintaining another source file with the sorting algorithm. Clearly if my business is "sorting" then I will spend my days finding how to sort in the best way... (I won't probably go as far as this http://stackoverflow.com/questions/2786899/fastest-sort-of-fixed-length-6-int-array but at least I will compare good implementations of the "standard" sort algorithms). My POV is that the OP will be more than happy with the default `qsort`. – xanatos Mar 15 '12 at 14:18
  • 1
    did you notice that he says "Quicksort will not perform that well" . On the rest I agree with you. –  Mar 15 '12 at 14:26
  • @g24l One thing is saying "quicksort will not perform that well", another is "the `qsort` of gnu c library takes x microseconds. I think it's too much. What can I do?". The difference between the two sentences is VERY VERY big. A naive quicksort is probably an order of time slower than the qsort of glibc. I gave a response to the first (implicit) question, someone else can give a response to the second (implicit) question. – xanatos Mar 15 '12 at 14:31
  • Sorry, but when I read questions I stick to the explicit part. –  Mar 15 '12 at 15:17
  • @g24l I have rechecked his sentence: `Heapsort seems a bit overkill. Quicksort will not perform that well.` Now... Does the c runtime have a heapsort sorting algo? No, so he isn't speaking of c library sorting algos, otherwise he would be comparing apples and oranges. – xanatos Mar 15 '12 at 15:21
  • Standard libraries implement bubblesort instead of quicksort() if they "sense" small arrays. At least linux ones did few years ago (for trully small arrays it's pretty fast and overhead is very low). Try some sample quicksort over small array and let gcc produce source, it will be 90% implemented as bubble sort. – Tomas Pruzina Mar 16 '12 at 00:40
1

Theoretically speaking, you want to find a median, which is a particular case of a selection problem. It may be done in linear time.

But this is in theory. Practically I'd use "partition-based general selection algorithm" (mentioned in the article). It works in linear time in average, plus it's practically simple and fast.

template <typename T>
T* SubDivide(T* pB, T* pE)
{
    ASSERT(pB < pE);

    T* pPivot = --pE;
    const T pivot = *pPivot;

    while (pB < pE)
    {
        if (*pB > pivot)
        {
            --pE;
            std::swap(*pB, *pE);
        } else
            ++pB;
    }

    std::swap(*pE, *pPivot);
    return pE;
}

template <typename T>
void SelectElement(T* pB, T* pE, size_t k)
{
    ASSERT((pE > pB) && (k < unsigned(pE - pB)));

    while (true)
    {
        T* pPivot = SubDivide(pB, pE);
        size_t n = pPivot - pB;

        if (n == k)
            break;

        if (n > k)
            pE = pPivot;
        else
        {
            pB = pPivot + 1;
            k -= (n + 1);
        }
    }
}

// Usage example
int pArr[9] = { /* fill with values*/ };

// Select the 4th element (which is a median)
SelectElement(pArr, pArr + 9, 4);

// pArr[4] now holds the median value
valdo
  • 12,632
  • 2
  • 37
  • 67