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 (see Median filter)
(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>
C
Support defines and functions:
<c>#define XGET(X) (((X)<0)?0:(((X)>=m->width)?m->width-1:(X)))
- define YGET(Y) (((Y)<0)?0:(((Y)>=m->height)?m->height-1:(Y)))
struct _pm {
pixel p; luminance lum;
};
static int _cmp(const void *a, const void *b) {
struct _pm *ap, *bp; ap = (struct _pm *)a; bp = (struct _pm *)b; if ( ap->lum > bp->lum ) return 1; if ( ap->lum < bp->lum ) return -1; return 0;
}
static void _get_median(image m,
unsigned int x, unsigned int y,
int r, pixel *p)
{
struct _pm *l;
int i, j;
unsigned int a,k;
l = malloc((2*r+1)*(2*r+1)*sizeof(struct _pm));
if ( l != NULL )
{
a = 0;
for(i=-r; i <= r; i++)
{
for(j=-r; j <= r; j++)
{
for(k=0; k < 3; k++)
{
l[a].p[k] = GET_PIXEL(m,XGET(x+i),YGET(y+j))[k];
}
l[a].lum = (2126*l[a].p[0] + 7152*l[a].p[1] +
722*l[a].p[2]) / 10000;
a++;
}
}
qsort(l, (2*r+1)*(2*r+1), sizeof(struct _pm), _cmp);
for(k=0;k<3;k++)
{
(*p)[k] = l[r].p[k];
}
free(l);
}
}</c>
The median filter:
<c>image median_filter(image si, int r) {
unsigned int x, y;
image dst;
pixel op;
dst = alloc_img(si->width, si->height);
if ( dst != NULL )
{
for(x=0; x < si->width; x++)
{
for(y=0; y < si->height; y++)
{
_get_median(si, x, y, r, &op);
put_pixel_unsafe(dst, x, y, op[0], op[1], op[2]);
}
}
}
return dst;
}</c>
Usage example:
<c>#include <stdio.h> /* #include "imglib.h" */
int main() {
image input, ei;
input = get_ppm(stdin);
if ( input == NULL ) return 1;
ei = median_filter(input, 2);
free_img(input);
if ( ei != NULL )
{
output_ppm(stdout, ei);
free_img(ei);
}
}</c>
This reads from stdin and writes to stdout.
The median filter, if implemented this way, can be very slow for big radius values or big images. This link taken from Wikipedia shows an algorithm that is O(1)!
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>