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
<lang 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; </lang> 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: <lang 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);
</lang>
C
Interface:
<lang c>image median_filter(image si, int r);</lang>
Support defines and functions:
<lang c>#include "imglib.h"
- 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);
}
}</lang>
The median filter:
<lang 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;
}</lang>
Usage example:
<lang 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);
}
}</lang>
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)!
D
This uses structures defined on Basic bitmap storage problem page.
The implementation uses algorithm described in Median Filtering in Constant Time paper with some slight differences, that shouldn't have impact on complexity.
Unfortunatelly because Bitmap's structure overloaded opIndex operator (as defined in Basic bitmap storage), doesn't return reference (since you can't return a reference in D 1.0, but this handy feature is present in D 2.0), calculations are made directly on underlying data buffer, good side of that is this is probably faster.
<lang D> ubyte median(uint[] cummulative, int no) {
int localSum;
foreach (k, v; cummulative)
if (v) {
localSum += v;
if (localSum > no / 2)
return k;
}
return 0;
}
RgbBitmap filterChannel(int Radius = 10)(RgbBitmap rgb, uint idx) {
auto res = RgbBitmap(rgb.width, rgb.height); res.data = rgb.data.dup; int edge = 2*Radius + 1, rsqr = edge*edge; alias uint[256] Hist; Hist H; Hist[] hCol;
static assert ( Radius >= 1 );
hCol.length = rgb.width;
// create histogram columns
for (int i = 0; i < edge - 1; i++)
for (int j = 0; j < rgb.width; j++)
hCol[j][ rgb.data[i*rgb.width + j].value[idx] ]++;
for (int i = Radius; i < rgb.height - Radius; i++) {
// add to each histogram column lower pixel
for (int j = 0; j < rgb.width; j++) {
hCol[j][ rgb.data[(i + Radius)*rgb.width + j].value[idx] ]++;
}
// calculate main Histogram using first edge-1 columns
H[] = 0;
for (int x = 0; x < edge - 1; x++)
foreach(k,v; hCol[x]) if (v)
H[k] += v;
for (int j = Radius; j < rgb.width - Radius; j++) {
// add right-most column
foreach (k, v; hCol[j + Radius]) if (v)
H[k] += v;
res.data[i*res.width + j].value[idx] = median(H, rsqr);
// drop left-most column
foreach (k, v; hCol[j - Radius]) if (v)
H[k] -= v;
}
// substract the upper pixels
for (int j = 0; j < rgb.width; j++)
hCol[j][ rgb.data[(i - Radius)*rgb.width + j].value[idx] ]--;
}
return res;
} </lang>
Sample usage: <lang D> alias filterChannel!(10) median10; auto inFile = new P6Image(new FileConduit("input.ppm"));
auto c1 = inFile.bitmap; auto c2 = median10(c1, 0); // red auto c3 = median10(c2, 1); // green auto c4 = median10(c3, 2); // blue
// this requires constructor defined in Write ppm file auto outFile = new P6Image(c4, inFile.maxVal); write (outFile); </lang>
OCaml
<lang 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)</lang>
an alternate version of the function median_value using arrays instead of lists: <lang 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)</lang>
Ruby
<lang ruby>class Pixmap
def median_filter(radius=3)
radius += 1 if radius.even?
filtered = self.class.new(@width, @height)
pb = ProgressBar.new(@height) if $DEBUG
@height.times do |y|
@width.times do |x|
window = []
(x - radius).upto(x + radius).each do |win_x|
(y - radius).upto(y + radius).each do |win_y|
win_x = 0 if win_x < 0
win_y = 0 if win_y < 0
win_x = @width-1 if win_x >= @width
win_y = @height-1 if win_y >= @height
window << self[win_x, win_y]
end
end
# median
filtered[x, y] = window.sort[window.length / 2]
end
pb.update(y) if $DEBUG
end
pb.close if $DEBUG
filtered
end
end
class RGBColour
# refactoring def luminosity Integer(0.2126*@red + 0.7152*@green + 0.0722*@blue) end def to_grayscale l = luminosity self.class.new(l, l, l) end
# defines how to compare (and hence, sort) def <=>(other) self.luminosity <=> other.luminosity end
end
class ProgressBar
def initialize(max)
$stdout.sync = true
@progress_max = max
@progress_pos = 0
@progress_view = 68
$stdout.print "[#{'-'*@progress_view}]\r["
end
def update(n)
new_pos = n * @progress_view/@progress_max
if new_pos > @progress_pos
@progress_pos = new_pos
$stdout.print '='
end
end
def close $stdout.puts '=]' end
end
bitmap = Pixmap.open('file') filtered = bitmap.median_filter</lang>
Tcl
<lang tcl>package require Tk
- Set the color of a pixel
proc applyMedian {srcImage x y -> dstImage} {
set x0 [expr {$x==0 ? 0 : $x-1}]
set y0 [expr {$y==0 ? 0 : $y-1}]
set x1 $x
set y1 $y
set x2 [expr {$x+1==[image width $srcImage] ? $x : $x+1}]
set y2 [expr {$y+1==[image height $srcImage] ? $y : $y+1}]
set r [set g [set b {}]]
foreach X [list $x0 $x1 $x2] {
foreach Y [list $y0 $y1 $y2] { lassign [$srcImage get $X $Y] rPix gPix bPix lappend r $rPix lappend g $gPix lappend b $bPix }
} set r [lindex [lsort -integer $r] 4] set g [lindex [lsort -integer $g] 4] set b [lindex [lsort -integer $b] 4] $dstImage put [format "#%02x%02x%02x" $r $g $b] -to $x $y
}
- Apply the filter to the whole image
proc medianFilter {srcImage {dstImage ""}} {
if {$dstImage eq ""} {
set dstImage [image create photo]
}
set w [image width $srcImage]
set h [image height $srcImage]
for {set x 0} {$x < $w} {incr x} {
for {set y 0} {$y < $h} {incr y} { applyMedian $srcImage $x $y -> $dstImage }
} return $dstImage
}
- Demonstration code using the Tk widget demo's teapot image
image create photo teapot -file $tk_library/demos/images/teapot.ppm pack [labelframe .src -text Source] -side left pack [label .src.l -image teapot] update pack [labelframe .dst -text Median] -side left pack [label .dst.l -image [medianFilter teapot]]</lang>