Thread Subject: 3D Median Filter

Hi folks,


I need a median filter which works on a 3D array. The bare bones
version would replace each element in the array with the median of
the 3x3x3 "cube" that surrounds that element.


Here's my SLOOOOOOOOW non vectorized version:


%******
function VM = medfilt3(V);


if ndims(V)~=3
  error('V must be a 3D array.');
end


% Pad all around with zeros
sz = size(V);
szpad = sz + [2 2 2];
Vw = zeros(szpad);
VM = Vw;
Vw(2:end-1,2:end-1,2:end-1)=V;


% Slowly loop through all the elements, performing the filter
for i = 2:(sz(1)+1),
    for j = 2:(sz(2)+1),
        for k = 2:(sz(3)+1),
            data = Vw((i-1):(i+1),(j-1):(j+1),(k-1:k+1));
            VM(i,j,k)=median(data(:));
        end
    end
end


% Downsize the output
VM = VM(2:end-1,2:end-1,2:end-1);
%********


Anyone have a suggestion for speeding this up?


Thanks,


Larry Schultz

"Larry Schultz" <schultzl@ece.pdx.edu> wrote in message
news:eebf5c2.-1@WebX.raydaftYaTP...
> Hi folks,
>
>
> I need a median filter which works on a 3D array. The bare bones
> version would replace each element in the array with the median of
> the 3x3x3 "cube" that surrounds that element.
>
>
> Here's my SLOOOOOOOOW non vectorized version:
>
>
> %******
> function VM = medfilt3(V);
>
>
> if ndims(V)~=3
> error('V must be a 3D array.');
> end
>
>
> % Pad all around with zeros
> sz = size(V);
> szpad = sz + [2 2 2];
> Vw = zeros(szpad);
> VM = Vw;
> Vw(2:end-1,2:end-1,2:end-1)=V;
>
>
> % Slowly loop through all the elements, performing the filter
> for i = 2:(sz(1)+1),
> for j = 2:(sz(2)+1),
> for k = 2:(sz(3)+1),
> data = Vw((i-1):(i+1),(j-1):(j+1),(k-1:k+1));
> VM(i,j,k)=median(data(:));
> end
> end
> end
>
>
> % Downsize the output
> VM = VM(2:end-1,2:end-1,2:end-1);
> %********
>
>
> Anyone have a suggestion for speeding this up?
>
>
> Thanks,
>
>
> Larry Schultz

I don't have a solution, but here's how I might start...

I'd observe that the arrays are stored linearly so that, except at
boundaries, the linear indices of subsequent blocks are offset by 1. I'd use
this fact to set up an indexing vector to which I'd increment an offset as I
step through the array. I've provided a simple example. You would have to
loop through all the elements, but at least you wouldn't have to do a
complex index calculation each time... It might help.

You might even find a way to vectorize the whole thing by computing all the
indices in one giant array, use it to make a giant array with all your voxel
values and then use median on all the columns at once. Of course, the tricky
part will be to handle the boundary conditions. You may want to ignore the
issue in your indexing and compute the wrapped-around medians with all the
others and then remove them later.

In any case, here is a very simplified example.

>> x = reshape(1:25,5,5)

x =

     1 6 11 16 21
     2 7 12 17 22
     3 8 13 18 23
     4 9 14 19 24
     5 10 15 20 25

>> [J,I]=meshgrid(1:3,1:3)

J =

     1 2 3
     1 2 3
     1 2 3


I =

     1 1 1
     2 2 2
     3 3 3

>> K=sub2ind(size(x),I,J)

K =

     1 6 11
     2 7 12
     3 8 13

>> K = K(:)

K =

     1
     2
     3
     6
     7
     8
    11
    12
    13

>> x(K)

ans =

     1
     2
     3
     6
     7
     8
    11
    12
    13

>> x(K+1)

ans =

     2
     3
     4
     7
     8
     9
    12
    13
    14

>> x(K+2)

ans =

     3
     4
     5
     8
     9
    10
    13
    14
    15

>> x

x =

     1 6 11 16 21
     2 7 12 17 22
     3 8 13 18 23
     4 9 14 19 24
     5 10 15 20 25

>> median(x(K+2))

ans =

     9

>> K = cumsum([K,ones(9,12)],2)

K =

     1 2 3 4 5 6 7 8 9 10 11 12
13
     2 3 4 5 6 7 8 9 10 11 12 13
14
     3 4 5 6 7 8 9 10 11 12 13 14
15
     6 7 8 9 10 11 12 13 14 15 16 17
18
     7 8 9 10 11 12 13 14 15 16 17 18
19
     8 9 10 11 12 13 14 15 16 17 18 19
20
    11 12 13 14 15 16 17 18 19 20 21 22
23
    12 13 14 15 16 17 18 19 20 21 22 23
24
    13 14 15 16 17 18 19 20 21 22 23 24
25

>> median(x(K))

ans =

     7 8 9 10 11 12 13 14 15 16 17 18
19

"Larry Schultz" <schultzl@ece.pdx.edu> writes:

> Hi folks,
>
>
> I need a median filter which works on a 3D array. The bare bones
> version would replace each element in the array with the median of
> the 3x3x3 "cube" that surrounds that element.

A high-efficiency/high-effort is to use the source for the MEX
function ordfilt2 in the image processing toolbox. It is not a
simple function, but is probably extendible to 3D arrays.


--
Peter Boettcher <boettcher@ll.mit.edu>
MIT Lincoln Laboratory
MATLAB FAQ: http://www.mit.edu/~pwb/cssm/

"Ken Davis" <kendavis@alum.mit.edu> wrote in message
news:55F969AE19647F250D8213758BB01024@in.WebX.raydaftYaTP...
>
> "Larry Schultz" <schultzl@ece.pdx.edu> wrote in message
> news:eebf5c2.-1@WebX.raydaftYaTP...
> > Hi folks,
> >
> >
> > I need a median filter which works on a 3D array. The bare bones
> > version would replace each element in the array with the median of
> > the 3x3x3 "cube" that surrounds that element.
> >
> >
> > Here's my SLOOOOOOOOW non vectorized version:
> >
> >
> > %******
> > function VM = medfilt3(V);
> >
> >
> > if ndims(V)~=3
> > error('V must be a 3D array.');
> > end
> >
> >
> > % Pad all around with zeros
> > sz = size(V);
> > szpad = sz + [2 2 2];
> > Vw = zeros(szpad);
> > VM = Vw;
> > Vw(2:end-1,2:end-1,2:end-1)=V;
> >
> >
> > % Slowly loop through all the elements, performing the filter
> > for i = 2:(sz(1)+1),
> > for j = 2:(sz(2)+1),
> > for k = 2:(sz(3)+1),
> > data = Vw((i-1):(i+1),(j-1):(j+1),(k-1:k+1));
> > VM(i,j,k)=median(data(:));
> > end
> > end
> > end
> >
> >
> > % Downsize the output
> > VM = VM(2:end-1,2:end-1,2:end-1);
> > %********
> >
> >
> > Anyone have a suggestion for speeding this up?
> >
> >
> > Thanks,
> >
> >
> > Larry Schultz
>
> I don't have a solution, but here's how I might start...
>
> I'd observe that the arrays are stored linearly so that, except at
> boundaries, the linear indices of subsequent blocks are offset by 1. I'd
use
> this fact to set up an indexing vector to which I'd increment an offset as
I
> step through the array. I've provided a simple example. You would have to
> loop through all the elements, but at least you wouldn't have to do a
> complex index calculation each time... It might help.
>
> You might even find a way to vectorize the whole thing by computing all
the
> indices in one giant array, use it to make a giant array with all your
voxel
> values and then use median on all the columns at once. Of course, the
tricky
> part will be to handle the boundary conditions. You may want to ignore the
> issue in your indexing and compute the wrapped-around medians with all the
> others and then remove them later.
>
> In any case, here is a very simplified example.
>
[snip]

In my work with median filters (which can be wonderful for outlier
rejection!) the proper way to handle bound conditions is to reflect the data
about the edges. (I should add that I like have the size of the data coming
out of a filter equal to the size of the data going into the filter.)
e.g. is you input data looks like:
[a b c d e f g h]
and you have a median filter window size of, say L=5, then reflect (L-1)/2
points as shown:
[c b a b c d e f g h g f]
and now your filtered data will look like:
[median([c b a b c], ...
median([b a b c d], ...
median([b a b c d], ...
median([b a b c d], ...
median([b a b c d], ...
median([b a b c d], ...

"AJ "no z" Johnson" <aj.jozhnson@lmco.com> wrote in message
news:bcd2rg$iur2@cui1.lmms.lmco.com...
> "Ken Davis" <kendavis@alum.mit.edu> wrote in message
> news:55F969AE19647F250D8213758BB01024@in.WebX.raydaftYaTP...
> >
> > "Larry Schultz" <schultzl@ece.pdx.edu> wrote in message
> > news:eebf5c2.-1@WebX.raydaftYaTP...
> > > Hi folks,
> > >
> > >
> > > I need a median filter which works on a 3D array. The bare bones
> > > version would replace each element in the array with the median of
> > > the 3x3x3 "cube" that surrounds that element.
> > >
> > >
> > > Here's my SLOOOOOOOOW non vectorized version:
> > >
> > >
> > > %******
> > > function VM = medfilt3(V);
> > >
> > >
> > > if ndims(V)~=3
> > > error('V must be a 3D array.');
> > > end
> > >
> > >
> > > % Pad all around with zeros
> > > sz = size(V);
> > > szpad = sz + [2 2 2];
> > > Vw = zeros(szpad);
> > > VM = Vw;
> > > Vw(2:end-1,2:end-1,2:end-1)=V;
> > >
> > >
> > > % Slowly loop through all the elements, performing the filter
> > > for i = 2:(sz(1)+1),
> > > for j = 2:(sz(2)+1),
> > > for k = 2:(sz(3)+1),
> > > data = Vw((i-1):(i+1),(j-1):(j+1),(k-1:k+1));
> > > VM(i,j,k)=median(data(:));
> > > end
> > > end
> > > end
> > >
> > >
> > > % Downsize the output
> > > VM = VM(2:end-1,2:end-1,2:end-1);
> > > %********
> > >
> > >
> > > Anyone have a suggestion for speeding this up?
> > >
> > >
> > > Thanks,
> > >
> > >
> > > Larry Schultz
> >
> > I don't have a solution, but here's how I might start...
> >
> > I'd observe that the arrays are stored linearly so that, except at
> > boundaries, the linear indices of subsequent blocks are offset by 1. I'd
> use
> > this fact to set up an indexing vector to which I'd increment an offset
as
> I
> > step through the array. I've provided a simple example. You would have
to
> > loop through all the elements, but at least you wouldn't have to do a
> > complex index calculation each time... It might help.
> >
> > You might even find a way to vectorize the whole thing by computing all
> the
> > indices in one giant array, use it to make a giant array with all your
> voxel
> > values and then use median on all the columns at once. Of course, the
> tricky
> > part will be to handle the boundary conditions. You may want to ignore
the
> > issue in your indexing and compute the wrapped-around medians with all
the
> > others and then remove them later.
> >
> > In any case, here is a very simplified example.
> >
> [snip]
>
> In my work with median filters (which can be wonderful for outlier
> rejection!) the proper way to handle bound conditions is to reflect the
data
> about the edges. (I should add that I like have the size of the data
coming
> out of a filter equal to the size of the data going into the filter.)
> e.g. is you input data looks like:
> [a b c d e f g h]
> and you have a median filter window size of, say L=5, then reflect (L-1)/2
> points as shown:
> [c b a b c d e f g h g f]
> and now your filtered data will look like:
> [median([c b a b c], ...
> median([b a b c d], ...
> median([b a b c d], ...
> median([b a b c d], ...
> median([b a b c d], ...
> median([b a b c d], ...
>
>
If you do use wrap around (which is only good for stationary/isotropic
processes), then you might find this mod1 function helpful for indexing:

function z=mod1(x,y)
%MOD1 Modulo with respect to interval [1, y] instead of interval [0, (y-1)].
z = 1 + mod(x-1, y);

if you wanted to use it to wrap around a 5-point data set with an extra two
points at the boundary, you would use:

>> data = 1:5

data =

     1 2 3 4 5

>> data = data(mod1(-1:7,5))

data =

     4 5 1 2 3 4 5 1 2

Ken

Peter Boettcher <boettcher@ll.mit.edu> writes:

> "Larry Schultz" <schultzl@ece.pdx.edu> writes:
>
> > Hi folks,
> >
> >
> > I need a median filter which works on a 3D array. The bare bones
> > version would replace each element in the array with the median of
> > the 3x3x3 "cube" that surrounds that element.
>
> A high-efficiency/high-effort is to use the source for the MEX
> function ordfilt2 in the image processing toolbox. It is not a
> simple function, but is probably extendible to 3D arrays.

Funny, it was just a couple of weeks ago that someone here created an
enhancement request to extend the IPT median- and order-statistic filtering
functions to arbitrary dimensions. We'll take a look at when we might be
able to do this.

--
Steve Eddins
The MathWorks, Inc.
[Mail sent to news4@eddins.net is forwarded to my MathWorks address.]

"Ken Davis" <kendavis@alum.mit.edu> wrote in
news:BED725561618F9A08233BC15D20287BD@in.WebX.raydaftYaTP:

> "AJ "no z" Johnson" <aj.jozhnson@lmco.com> wrote in message
> news:bcd2rg$iur2@cui1.lmms.lmco.com...
>> "Ken Davis" <kendavis@alum.mit.edu> wrote in message
>> news:55F969AE19647F250D8213758BB01024@in.WebX.raydaftYaTP...
>> >
>> > "Larry Schultz" <schultzl@ece.pdx.edu> wrote in message
>> > news:eebf5c2.-1@WebX.raydaftYaTP...
>> > > Hi folks,
>> > >
>> > >
>> > > I need a median filter which works on a 3D array. The bare bones
>> > > version would replace each element in the array with the median
>> > > of the 3x3x3 "cube" that surrounds that element.
>> > >
>> > >
>> > > Here's my SLOOOOOOOOW non vectorized version:
>> > >
>> > >
>> > > %******
>> > > function VM = medfilt3(V);
>> > >
>> > >
>> > > if ndims(V)~=3
>> > > error('V must be a 3D array.');
>> > > end
>> > >
>> > >
>> > > % Pad all around with zeros
>> > > sz = size(V);
>> > > szpad = sz + [2 2 2];
>> > > Vw = zeros(szpad);
>> > > VM = Vw;
>> > > Vw(2:end-1,2:end-1,2:end-1)=V;
>> > >
>> > >
>> > > % Slowly loop through all the elements, performing the filter
>> > > for i = 2:(sz(1)+1),
>> > > for j = 2:(sz(2)+1),
>> > > for k = 2:(sz(3)+1),
>> > > data = Vw((i-1):(i+1),(j-1):(j+1),(k-1:k+1));
>> > > VM(i,j,k)=median(data(:));
>> > > end
>> > > end
>> > > end
>> > >
>> > >
>> > > % Downsize the output
>> > > VM = VM(2:end-1,2:end-1,2:end-1);
>> > > %********
>> > >
>> > >
>> > > Anyone have a suggestion for speeding this up?
>> > >
>> > >
>> > > Thanks,
>> > >
>> > >
>> > > Larry Schultz
>> >
>> > I don't have a solution, but here's how I might start...
>> >
>> > I'd observe that the arrays are stored linearly so that, except at
>> > boundaries, the linear indices of subsequent blocks are offset by
>> > 1. I'd
>> use
>> > this fact to set up an indexing vector to which I'd increment an
>> > offset
> as
>> I
>> > step through the array. I've provided a simple example. You would
>> > have
> to
>> > loop through all the elements, but at least you wouldn't have to do
>> > a complex index calculation each time... It might help.
>> >
>> > You might even find a way to vectorize the whole thing by computing
>> > all
>> the
>> > indices in one giant array, use it to make a giant array with all
>> > your
>> voxel
>> > values and then use median on all the columns at once. Of course,
>> > the
>> tricky
>> > part will be to handle the boundary conditions. You may want to
>> > ignore
> the
>> > issue in your indexing and compute the wrapped-around medians with
>> > all
> the
>> > others and then remove them later.
>> >
>> > In any case, here is a very simplified example.
>> >
>> [snip]
>>
>> In my work with median filters (which can be wonderful for outlier
>> rejection!) the proper way to handle bound conditions is to reflect
>> the
> data
>> about the edges. (I should add that I like have the size of the data
> coming
>> out of a filter equal to the size of the data going into the filter.)
>> e.g. is you input data looks like:
>> [a b c d e f g h]
>> and you have a median filter window size of, say L=5, then reflect
>> (L-1)/2 points as shown:
>> [c b a b c d e f g h g f]
>> and now your filtered data will look like:
>> [median([c b a b c], ...
>> median([b a b c d], ...
>> median([b a b c d], ...
>> median([b a b c d], ...
>> median([b a b c d], ...
>> median([b a b c d], ...
>>
>>
> If you do use wrap around (which is only good for stationary/isotropic
> processes), then you might find this mod1 function helpful for
> indexing:
>
> function z=mod1(x,y)
> %MOD1 Modulo with respect to interval [1, y] instead of interval [0,
> (y-1)]. z = 1 + mod(x-1, y);
>
> if you wanted to use it to wrap around a 5-point data set with an
> extra two points at the boundary, you would use:
>
>>> data = 1:5
>
> data =
>
> 1 2 3 4 5
>
>>> data = data(mod1(-1:7,5))
>
> data =
>
> 4 5 1 2 3 4 5 1 2
>
> Ken
>
>

The reflection described, though, is 3 2 1 2 3 4 5 4 3!

--
Scott
Reverse first field of address to reply

[snip]
> >>
> >> In my work with median filters (which can be wonderful for outlier
> >> rejection!) the proper way to handle bound conditions is to reflect
> >> the
> > data
> >> about the edges. (I should add that I like have the size of the data
> > coming
> >> out of a filter equal to the size of the data going into the filter.)
[snip
>
> The reflection described, though, is 3 2 1 2 3 4 5 4 3!
>
> --
> Scott
> Reverse first field of address to reply

Oops...

I guess I had better read what the postings actually say instead of
presuming (often incorrectly) that I already know.

Ken

Steve Eddins <news4@eddins.net> wrote in message
news:uel1xhmj1.fsf@eddins.net...
> Peter Boettcher <boettcher@ll.mit.edu> writes:
>
> > "Larry Schultz" <schultzl@ece.pdx.edu> writes:
>>
> > > Hi folks,
> > > I need a median filter which works on a 3D array. The bare bones
> > > version would replace each element in the array with the median of
> > > the 3x3x3 "cube" that surrounds that element.
> >
> > A high-efficiency/high-effort is to use the source for the MEX
> > function ordfilt2 in the image processing toolbox. It is not a
> > simple function, but is probably extendible to 3D arrays.
>
> Funny, it was just a couple of weeks ago that someone here created an
> enhancement request to extend the IPT median- and order-statistic
filtering
> functions to arbitrary dimensions. We'll take a look at when we might
be
> able to do this.
>
> Steve Eddins
> The MathWorks, Inc.
> [Mail sent to news4@eddins.net is forwarded to my MathWorks address.]

I have not evaluated the image processing toolbox, but I have used and
examined the median .m files in the signal processing toolbox. An m
point running median window in an array produces an array in the
workspace m times the size of the filtered array (yes it can enlarge
only a chunk at a time) and each output point makes a call to sort.

Creating your own MEX file to perform a running mean with a fixed size
work area and no complete sorts required can be much faster. Another
approach to the edge processing is to grow the window. For the running
calculation I suggest the following references:

A review of approaches:

M. Juhola, J. Katajainen, and T. Raita, "Comparison of algorithms for
standard median filtering," IEEE Trans. Signal Processing, vol. 39, pp.
204 - 208, January 1991.
ABSTRACT:
In standard median filtering we search repeatedly for a median from a
sample set which changes only slightly between the subsequent searches.
We review several well-known methods for solving this running median
problem, analyze the (asymptotical) time complexities of the methods,
and propose simple variants which are especially suited for small sample
sets, a frequent situation. Although we have restricted our discussion
to the one-dimensional case, the ideas are easily extended to higher
dimensions.

This one I've used for floating point data since the 1990s:

J. T. Astola and T. G. Campbell, "On computation of the running median,"
IEEE Trans. Acoust., Speech, Signal Processing, vol. 37, pp. 572 - 574,
April 1989.
ABSTRACT:
This correspondence presents a fast median filtering algorithm with
logarithmic time complexity based on a special data structure, a double
heap, which naturally supports the median. With slight modification, the
approach can be used to implement any rank order filter.

This one I've used for fixed point (11 bit) data in custom hardware in
the 1980s and 1990s:

E. Ataman, V. K. Aatre, and K. M. Wong, "A fast method for real-time
median filtering," IEEE Trans. Acoust., Speech, Signal Processing, vol.
28, pp. 415 - 421, August 1980.
ABSTRACT:
A fast real-time algorithm is presented for median filtering of signals
and images. The algorithm determines the kth bit of the median by
inspecting the k most significant bits of the samples. The total number
of full-word comparison steps is equal to the wordlength of the samples.
Speed and hardware complexity of the algorithm is compared with two
other fast methods for median filtering.

Dale B. Dalrymple
Senior Staff Systems Engineer
Signal Processing Systems Div. of ISL Inc.

"Dale B. Dalrymple" <dbd@abac.com> writes:

> I have not evaluated the image processing toolbox, but I have used and
> examined the median .m files in the signal processing toolbox. An m
> point running median window in an array produces an array in the
> workspace m times the size of the filtered array (yes it can enlarge
> only a chunk at a time) and each output point makes a call to sort.

For arbitrarily-shaped neighborhoods, ORDFILT2 uses an order(N) selection
algorithm based on quicksort-style partitioning.

For large rectangular neighborhoods, ORDFILT2 uses a histogram-based
method inspired by "A fast two-dimensional median filtering algorithm,"
T.S.Huang, G.J.Yang, G.Y.Tang, IEEE Transactions on Acoustics, Speech and
Signal Processing, Vol ASSP 27 No 1, February 1979.

Both algorithms are implemented as MEX-files. MEDFILT2 calls ORDFILT2.

--
Steve Eddins
The MathWorks, Inc.
[Mail sent to news4@eddins.net is forwarded to my MathWorks address.]