[R] Asking Favor For the Script of Median Filter

[David Winsemius]

On Mar 27, 2011, at 10:56 AM, chuan_zl wrote:

> Hello,everybody. My name is Chuan Zun Liang. I come from Malaysia. I  
> am just
> a beginner for R. Kindly to ask favor about median filter. The  
> problem I
> facing as below:
>
>
>> x<-matrix(sample(1:30,25),5,5)
>> x
>     [,1] [,2] [,3] [,4] [,5]
> [1,]    7    8   30   29   13
> [2,]    4    6   12    5    9
> [3,]   25    3   22   14   24
> [4,]    2   15   26   23   19
> [5,]   28   18   10   11   20
>
> This is example original matrices of an image. I want apply with  
> median
> filter with window size 3X# to remove salt and pepper noise in my  
> matric.
> Here are the script I attend to writing.The script and output shown as
> below:
>
>> MedFilter<-function(mat,sz)
> + {out<-matrix(0,nrow(mat),ncol(mat))
> +  for(p in 1:(nrow(mat)-(sz-1)))
> + {for(q in 1:(ncol(mat)-(sz-1)))
> + {outrow<-median(as.vector(mat[p:(p+(sz-1)),q:(q+(sz-1))]))
> + out[(p+p+(sz-1))/2,(q+q+(sz-1))/2]<-outrow}}
> + out}
>

Noting that median is probably the "rate-limiting" factor, I looked  
for another way to get the "middle" of 9 items. Using order seem faster:

 > system.time( replicate(100000, x2s <- sort(x2)))
    user  system elapsed
   9.829   0.212  10.029
 > system.time( replicate(100000, x2m <- median(x2)))
    user  system elapsed
   7.169   0.126   7.272
 > system.time( replicate(100000, x2s <-x2[order(x2)[5] ]))
    user  system elapsed
   1.907   0.051   1.960

So see if this is any faster. On my system it's about three times  
faster:

x <- matrix(sample(364*364), 364,364)
out <- matrix(0, 364,364)
for(xi in 1:(nrow(x)-2)) {
       for(yi in 1:(ncol(x)-2) ) {
          xm <- x[xi+0:2, yi+0:2]
           d[xi+1, yi+1] <-xm[order(xm)[5] ]}}

#---------tests ------
 > system.time(for(xi in 1:(nrow(x)-2)) {
+       for(yi in 1:(ncol(x)-2) ) {
+          xm <- x[xi+0:2, yi+0:2]
+           d[xi+1, yi+1] <-xm[order(xm)[5] ]}} )
    user  system elapsed
   3.806   0.083   3.887

 > system.time(MedFilter(x,3) )
    user  system elapsed
  11.242   0.202  11.427


>> MedFilter(x,3)
>     [,1] [,2] [,3] [,4] [,5]
> [1,]    0    0    0    0    0
> [2,]    0    8   12   14    0
> [3,]    0   12   14   19    0
> [4,]    0   18   15   20    0
> [5,]    0    0    0    0    0
>
> Example to getting value 8 and 12 as below:
>
> 7  8 30                         8 30 29
> 4  6 12 (median=8)         6 12   5 (median=12)
> 25 3 22                         3 22 14
>
> Even the script can give output. However, it is too slow. My image  
> size is
> 364*364. It is time consumption. Is it get other ways to improving it?

David Winsemius, MD
West Hartford, CT