Median filter
You are encouraged to solve this task according to the task description, using any language you may know.
The median filter takes in the neighbourhood the median color.
(to test the function below, you can use these input and output solutions)
Ada
<ada> function Median (Picture : Image; Radius : Positive) return Image is
type Extended_Luminance is range 0..10_000_000; type VRGB is record Color : Pixel; Value : Luminance; end record; Width : constant Positive := 2*Radius*(Radius+1); type Window is array (-Width..Width) of VRGB; Sorted : Window; Next : Integer;
procedure Put (Color : Pixel) is -- Sort using binary search pragma Inline (Put); This : constant Luminance := Luminance ( ( 2_126 * Extended_Luminance (Color.R) + 7_152 * Extended_Luminance (Color.G) + 722 * Extended_Luminance (Color.B) ) / 10_000 ); That : Luminance; Low : Integer := Window'First; High : Integer := Next - 1; Middle : Integer := (Low + High) / 2; begin while Low <= High loop That := Sorted (Middle).Value; if That > This then High := Middle - 1; elsif That < This then Low := Middle + 1; else exit; end if; Middle := (Low + High) / 2; end loop; Sorted (Middle + 1..Next) := Sorted (Middle..Next - 1); Sorted (Middle) := (Color, This); Next := Next + 1; end Put; Result : Image (Picture'Range (1), Picture'Range (2));
begin
for I in Picture'Range (1) loop for J in Picture'Range (2) loop Next := Window'First; for X in I - Radius .. I + Radius loop for Y in J - Radius .. J + Radius loop Put ( Picture ( Integer'Min (Picture'Last (1), Integer'Max (Picture'First (1), X)), Integer'Min (Picture'Last (2), Integer'Max (Picture'First (2), Y)) ) ); end loop; end loop; Result (I, J) := Sorted (0).Color; end loop; end loop; return Result;
end Median; </ada> The implementation works with an arbitrary window width. The window is specified by its radius R>0. The resulting width is 2R+1. The filter uses the original pixels of the image from the median of the window sorted according to the luminance. The image edges are extrapolated using the nearest pixel on the border. Sorting uses binary search. (For practical use, note that median filter is extremely slow.)
The following sample code illustrates use: <ada>
F1, F2 : File_Type;
begin
Open (F1, In_File, "city.ppm"); Create (F2, Out_File, "city_median.ppm"); Put_PPM (F2, Median (Get_PPM (F1), 1)); -- Window 3x3 Close (F1); Close (F2);
</ada>
OCaml
<ocaml>let color_add (r1,g1,b1) (r2,g2,b2) =
( (r1 + r2), (g1 + g2), (b1 + b2) )
let color_div (r,g,b) d =
( (r / d), (g / d), (b / d) )
let compare_as_grayscale (r1,g1,b1) (r2,g2,b2) =
let v1 = (2_126 * r1 + 7_152 * g1 + 722 * b1) and v2 = (2_126 * r2 + 7_152 * g2 + 722 * b2) in (Pervasives.compare v1 v2)
let get_rgb img x y =
let _, r_channel,_,_ = img in let width = Bigarray.Array2.dim1 r_channel and height = Bigarray.Array2.dim2 r_channel in if (x < 0) || (x >= width) then (0,0,0) else if (y < 0) || (y >= height) then (0,0,0) else (* feed borders with black *) (get_pixel img x y)
let median_value img radius =
let samples = (radius*2+1) * (radius*2+1) in fun x y -> let sample = ref [] in
for _x = (x - radius) to (x + radius) do for _y = (y - radius) to (y + radius) do let v = get_rgb img _x _y in
sample := v :: !sample; done; done;
let ssample = List.sort compare_as_grayscale !sample in let mid = (samples / 2) in
if (samples mod 2) = 1 then List.nth ssample (mid+1) else let median1 = List.nth ssample (mid) and median2 = List.nth ssample (mid+1) in (color_div (color_add median1 median2) 2)
let median img radius =
let _, r_channel,_,_ = img in let width = Bigarray.Array2.dim1 r_channel and height = Bigarray.Array2.dim2 r_channel in
let _median_value = median_value img radius in
let res = new_img ~width ~height in for y = 0 to pred height do for x = 0 to pred width do let color = _median_value x y in put_pixel res color x y; done; done; (res)</ocaml>
an alternate version of the function median_value
using arrays instead of lists:
<ocaml>let median_value img radius =
let samples = (radius*2+1) * (radius*2+1) in let sample = Array.make samples (0,0,0) in fun x y -> let i = ref 0 in for _x = (x - radius) to (x + radius) do for _y = (y - radius) to (y + radius) do let v = get_rgb img _x _y in sample.(!i) <- v; incr i; done; done;
Array.sort compare_as_grayscale sample; let mid = (samples / 2) in
if (samples mod 2) = 1 then sample.(mid+1) else (color_div (color_add sample.(mid) sample.(mid+1)) 2)</ocaml>