OpenEuphoria

Re: Error in median()

• Posted by Spock May 06, 2020
• 1173 views

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