de Rigo, D. (2012). Multi-dimensional weighted median: the module "wmedian" of the Mastrave modelling library. Mastrave project technical report. http://mastrave.org/doc/mtv_m/wmedian

## Multi-dimensional weighted median: the module "wmedian" of the Mastrave modelling library

Daniele de Rigo

Abstract: Weighted median (WM) filtering is a well known technique for dealing with noisy images and a variety of WM-based algorithms have been proposed as effective ways for reducing uncertainties or reconstructing degraded signals by means of available information with heterogeneous reliability. Here a generalized module for applying weighted median filtering to multi-dimensional arrays of information with associated multi-dimensional arrays of corresponding weights is presented. Weights may be associated to single elements or to groups of elements along given dimensions of the multi-dimensional arrays. The filtered information derives from a reduction operator applied along a custom dimension.

#### Copyright and license notice of the function wmedian

The file wmedian.m is part of Mastrave.

Mastrave is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

Mastrave is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with Mastrave. If not, see http://www.gnu.org/licenses/.

#### Function declaration

```answer = wmedian( values        ,
dim     = []  ,
weights = []  )

```

#### Description

Utility to extend the function median(.) providing weighted medians of the elements of the array values along a given dimension dim .

The weighted median of a column vector v with integer weights w is equivalent to the median of the vector [ v(1)◇w(1) ; v(2)◇w(2) ; ... ], where the operator ◇ denotes duplications (Yin et al., 1996), i.e.

v(1)◇w(1) = repmat( v(1), w(1), 1 )

The weighted median wm of v with nonnegative weights w2 is defined as

wm = arg min( w2' * abs( v - wm ) )

If one or more weights are Inf, the corresponding elements of values aligned along the dimension dim are weighted as if the Inf weights were al having the same weight and all remaining weights were zeros. This does not affect elements of values which are not aligned with Inf-weighted elements. If one or more weights are NaN, they are considered as zero weights. In case one or more all-zeros sequences of weights are aligned along the dimension dim of weights , the corresponding elements of values are weighted as if all weights were uniform.

##### References

Yin, L., Yang, R., Gabbouj, M., Neuvo, M. (1996): Weighted Median
Filters: A Tutorial. IEEE Transactions on Circuits and Systems II:
Analog and Digital Signal Processing, Vol. 43, No. 3, pp. 157-192,
March 1996.
DOI: 10.1109/82.486465.
Free access version:
http://www.cs.tut.fi/~moncef/publications/weighted-median-filters.pdf

#### Input arguments

```
values            ::numeric::
Vector, matrix or multi-dimensional array of numbers.

dim               ::scalar_index|empty::
Scalar positive integer representing the dimension along
which the weighted medians have to be computed.
If  dim  is an empty array  [] , the dimension is the
first non-singleton dimension.  In case  values  is a
vector, this definition means that the default
dimension is the one along which the elements of the
vector  values  are aligned.
If omitted, the default value is  [].

weights           ::nonnegative::
Vector, matrix or multi-dimensional array of nonnegative
numbers representing the weights to be associeted to the
corresponding elements of  values .
If  weights  and  values  do not have the same size, they
are expected to be instances of flats (linear manifolds)
suitable to be combined within a bsxfun(...) call.
If  weights  is an empty array  [], then it is considered
as an array of ones with the same size as  values .
If omitted, the default value is  [].

```

#### Example of usage

```
% Basic usage

% Vectors:
v   = ceil( rand(1,7)* 100 )
wm  = wmedian( v )
assert( wm == median(v) )
w   = bsxfun( @power, 1:7, [0:4].' );
w   = [w(end:-1:2,end:-1:1); w]
wm  = wmedian( v,2,w )

% Verifying the definition of weighted median
def = @(v,w,x)abs( bsxfun( @minus, v(:).', x(:) ) )*w(:)
x   = [1:100].';
hold off
for i=1:size(w,1)
wi = w(i,:).';
mi = abs( v - wm(i) ) * wi;
semilogy( x, def(v,wi,x), wm(i), mi , 'or');
text( wm(i), mi*.8, sprintf( 'w( %d, : )', i ) )
hold on;
end; hold off

% Matrices:
v   = ceil( rand(5,7)* 100 )
wm  = wmedian( v )
assert( wm == median(v) )

% Passing a custom dimension
wm  = wmedian( v , 2 )
assert( wm == median(v,2) )

% Dealing with multi-dimensional arrays
v   = ceil( rand(5,7,3)* 100 )
wm  = wmedian( v , 1 )
assert( wm == median(v,1) )

wm  = wmedian( v , 2 )
assert( wm == median(v,2) )

```
```Memory requirements:
O( numel( bsxfun( @plus,  values ,  weights  ) ) )

cumstd, cumvar, cumsumsq, groupfun

Keywords:
weighted operators, reduction

Version: 0.5.5```

#### Support

The Mastrave modelling library is committed to provide reusable and general - but also robust and scalable - modules for research modellers dealing with computational science.  You can help the Mastrave project by providing feedbacks on unexpected behaviours of this module.  Despite all efforts, all of us - either developers or users - (should) know that errors are unavoidable.  However, the free software paradigm successfully highlights that scientific knowledge freedom also implies an impressive opportunity for collectively evolve the tools and ideas upon which our daily work is based.  Reporting a problem that you found using Mastrave may help the developer team to find a possible bug.  Please, be aware that Mastrave is entirely based on voluntary efforts: in order for your help to be as effective as possible, please read carefully the section on reporting problems.  Thank you for your collaboration.

Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Daniele de Rigo