Median value selection algorithm

At work I'm currently working on tomographic reconstruction algorithms. I have to implement a Bayesian iterative algorithm that requires to select the median value in a set of the 27 float values from a cube of 3x3x3 voxels. This operation must be performed for each voxel and for each iteration. We have to expect 256 million voxels to process for each iteration, but "only" 60 to 100 iterations.

Trying to find the most efficient algorithm I came up with a new algorithm considering the one described on the select algorithm page of wikipedia. I then submitted a question on StackOverflow to get some feedback. And I did get valuable feedback. I was first pointed to the C++ nth_element function I didn't know at the time. It was also suggest to optimize by sharing intermediate information which is indeed a smart thing to do to get a better initial guess on the median value.

After multiple tests and changes to the code I finally reached what seem to be a very efficient algorithm. The ratio is so good that I'm still unsure about it, but I checked everything. It could be due to memory cache or particularly favorable parallelization opportunities with sse instructions. I don't know.

Here is the outline of the algorithm. We have 27 values and have to find the value that split the set in two with 13 values smaller or equal to it and 13 values bigger or equal to it. This value is called the median value.

The fundamental idea of the algorithm is to use a heap data structure which has its smallest or biggest value at the top. Such data structure is very efficient for adding an element or to extract its top most element. It can also be very easily mapped into an array.

The algorithm uses two heaps initially empty, with a common top value which is the median value. Each heap has a capacity of at most 14 elements, where one is common to the two heaps, their top value and also the median value.
The algorithm proceed in two phases. In the first phase, the algorithm picks a value as initial median value guess. Subsequent values are then compare with this median value and added to the corresponding heap until one heap becomes full and contains 14 elements.

At this point the median value is a value in the full heap or in the remaining set of values to process. The second phase of the algorithm then starts where the remaining values are processed. Values that would not be inserted in the full heap are ignored. The other values are inserted in the full heap after deleting its top most value. The heap then gets a new top most value and thus also a new median value. When all the remaining elements have been processed, the median of the 27 value set is the top most value of the full heap.

Here are some benchmark results. See StackOverflow for more detailed information.

HeapSort        :2.287 
QuickSort       :2.297 
QuickMedian1 :0.967 
HeapMedian1  :0.858 
NthElement     :0.616 
QuickMedian2 :1.178 
HeapMedian2  :0.597 
HeapMedian3  :0.015  <-- best

It thus seem that HeapMedian3 is 33 times faster than NthElement. I used a 3GHz Intel E8400 processor and the Intel C++ compiler with options -03 and -xS for benchmarking.

Here is the code :

EDIT: The code has been removed because it was invalid. I tested it on a single random float value sequence were it gave a valid result by coincidence. A new blog post provides the correct code.