### Re: Error in median()

- Posted by Spock May 06, 2020
- 645 views

Quite so.

I've reverted to my own function, which I wrote way back before the stats package was available. To use this, I changed median() in the code I'm writing to Median(). Simple enough.

....

Hi Craig,

No doubt your use case differs to mine but there will be similarities. Recently I needed to use a **weighted** median function. This was for my graphics library to do basic image denoising. Maybe I didn't look hard enough but there didn't seem to be any simple but efficient algorithm on the internet. A naive way would be to take the list of each color channel and replicate each element based on its weight and then take the center element after sorting [you could do this as a quick 'n dirty way of weighting your median function]. But using this method, even for small kernels, the processing times soon become prohibitive. The largest kernel size I could go (using standard sorting) was a measly 7x7 pixels. This was sorting 4K elements per pixel. In a multi mega-pixel image that is not a sensible option. Fortunately, I can take advantage of counting sort using a 256 length sequence of accumulators for each of the 3 color channels. A single pass over the N x N kernel, adding the weights to each accumulator, followed by a linear search for the mid entry of each channel accumulator yields a very fast, but synthetic, median RGB value.

For highest quality I also added a function to select the centermost pixel in 3d space. There is no natural way to directly sort RGB pixels. The algorithm I use computes the accumulated distances from each pixel to all other pixels. Lowest accumulated value is deemed to be the median pixel. Unlike the other algorithm a real pixel value is actually returned. I use this code (in Orac) which is reminiscent of Bubble Sort in its n^2 time complexity. Although it's slow the quality seems to be the best. The trick to the weighting is to multiply the **other** pixel (see p1, p2 below) with the weights associated to a pixel:

global function weighted_median_rgb(seq s, seq wgts) -- get the median rgb value from the list of pixels -- this method is not efficient.. -- but unlike *synth() it picks a value rather than compute one .. -- the method is to compute the sum of euclidean distance between all points -- lowest overall score wins -- init seq d = repeat(0, ~s) -- distances -- loop for i = 1 to ~s-1 do int p1 = s.i int w1 = wgts.i for j = i+1 to ~s do int p2 = s.j int w2 = wgts.j atom dist = ColourDiff2(p1, p2) d.i += dist * w2 -- see how they are opposites here.. d.j += dist * w1 end for end for -- find minimum int n = 1 atom min = d.1 for i = 2 to ~d do if d.i < d.n then n = i end if end for -- exit return s.n end function

Spock