Hampel filter function

Share icon

Tue, 10/10/2023 - 11:34 pm

I wrote a function for Hampel filter. If anyone can help me check whether it is correct or nor? Thanks a lot.

function hampel(imgraw,k,nd)
wave imgraw
variable k,nd
duplicate/o imgraw,imgraw_smth
variable lx = dimsize(imgraw,1)

make/o/n=(lx) iLo=p-k
make/o/n=(lx) iHi=p+k
variable i,j
for(i=0;i<lx;i++)
  if(iLo[i]<0)
    iLo[i] = 0
  endif
endfor
for(i=0;i<lx;i++)
  if(iHi[i]>(lx-1))
    iHi[i] = lx-1
  endif
endfor
for(i=0;i<dimsize(imgraw,0);i++)
  for(j = 0;j<lx;j++)
    duplicate/o/rmd=[i,i][iLo(j),iHi(j)] imgraw,temp
    matrixtranspose temp
    make/n=(dimsize(temp,0))/d/o w
    w=temp
    variable medj = median(w)
    make/d/o/n=(lx) mmed,mmad
    mmed[j] = medj
    make/d/o/n=(iHi(j)-iLo(j)+1) middle
    middle=abs(w-medj)
    mmad[j] = median(middle)

    make/d/o/n=(lx) sd
    sd[j]= mmad[j]/sqrt(2)/inverseErf(0.5)
    make/d/o/n=(lx) ki
    ki[j]= abs(imgraw[i][j]-mmed[j]) > nd*sd[j]

    make/d/o/n=(lx) yu
 
    if(ki[j])
      yu[j] = mmed[j]
    else
      yu[j]= imgraw[i][j]
    endif
    imgraw_smth[i][j]=yu[j]
  endfor
endfor
killwaves w,temp,mmed,mmad,middle,sd,ki,yu
end

 

A suggestion for the validity test is to compare results from your algorithm against those using posted code for other analysis packages (e.g. python).

As to the code itself, I offer these suggestions

--> replace for loops by implicit statements

// not this
for(i=0;i<lx;i++)
  if(iLo[i]<0)
    iLo[i] = 0
  endif
endfor

// instead this
iLo = iLo < 0 ? 0 : iLo

--> use FREE waves, set waves outside loop when possible, and avoid making a wave just as a placeholder for a one step processing

// not this
    make/n=(dimsize(temp,0))/d/o w
    w=temp
    variable medj = median(w)
    make/d/o/n=(lx) mmed,mmad
    mmed[j] = medj
    make/d/o/n=(iHi(j)-iLo(j)+1) middle
    middle=abs(w-medj)
    mmad[j] = median(middle)

// instead this (or equivalent)

    make/N=(lx)/D/FREE mmed, mmad // this should be outside loop??
    med[j] = median(temp)

 

October 11, 2023 at 08:50 am - Permalink

Disclaimer: I have not tested your code and I have no direct experience with this filter.

Suggestions:

1.  it is not a good idea to have 6 Make statements inside an interior loop.  As long as the dimensions of the waves remain fixed (with respect to the loop), you should create them (execute the Make statement) only once before getting into the loop.  Note also that some of these waves, e.g., mmad and sd may be replaced by local variables.

2.  it is not a good idea to keep computing in the inner loop

sqrt(2)/inverseErf(0.5)

You should simply compute one scale factor e.g., 1/sqrt(2)/inverseERF(0.5) once before executing the loops.

 

 

October 11, 2023 at 11:12 am - Permalink

Thank you! I have compared my procedure with Matlab code and found they agree well. I used gnoise in Igor and copy the generated data into Matlab as orginal data and obtained the same filtered data using my Igor procedure and Matlab.
095502.png

October 11, 2023 at 06:13 pm - Permalink

Just curious. Your procedure seems to be designed for image data. Is your test from a line profile on an image?

October 12, 2023 at 06:39 am - Permalink

FWIW, I believe your filter function will not work properly on 2-D images. When you input an MxN matrix, the function iterates over the M lines in a 1D manner, looking only along the N points in the single line for a median (average). To work with images, the function should instead iterate over each point location in the MxN image matrix, analyzing in both the x and the y directions to obtain a representative average (median) for the given point.

If all you want is a filter to work on 1-D waves, you might consider removing the loop on the additional dimension (and also for example renaming the input wave as vectorraw rather than imgraw).

Also, in Igor Pro matrix notations, a 1-D wave is a row vector and an (M x N) matrix has M rows and N columns. In MatLab matrix notations, a 1-D input is a column vector and an (M x N) matrix has M columns and N rows. See the discussion at the link below.

https://www.wavemetrics.com/forum/general/translating-matlab-matrix-mul…

October 14, 2023 at 08:16 am - Permalink