pro denoise,thresh,p,n,n1,im, bestim

; version du 17 mars 1998

;Denoising by wavelet thresholding using Haar filter

; thresh='soft' or 'hard'
; p  - threshold
; n  - number of levels of decomposition
; n1 - image size
; im     is input  data and it must be bytarr(n1,n1)
; bestim is output data and it is also bytarr(n1,n1)

  print,'Denoising by wavelet ',thresh,' thresholding...'

 l  = 2
 h  = fltarr(l)
 th = fltarr(l)
 g  = fltarr(l)
 tg = fltarr(l)

 filtername='Haar'
 print,'using the ',filtername,' filter'

 h (0)=0.5
 h (1)=0.5
 th(0)=0.5
 th(1)=0.5

 k=1
 for i=0,l-1 do begin
 tg(i)=h(l-1-i)*k
 g(i)=th(l-1-i)*k
 k=-k
 endfor

soft = 'soft'
hard = 'hard'

 n2 = n1
 npix = long(n1)*n2

 r = fltarr(n1,n2)
 a = fltarr(n1,n2)

 s = n1

 a = float(im)

; multiscale decomposition

 for e=n,1,-1 do begin
 is=s-1
 s2=s/2
 is2=s2-1

 r(0:is,0:is)=0.0

 ind=2*indgen(s2)

 for k=0,l-1 do begin
 indk=(ind+k) mod s
 r(0:is,0:is2)=r(0:is,0:is2)+th(k)*a(0:is,indk)*2.0
 r(0:is,s2:is)=r(0:is,s2:is)+tg(k)*a(0:is,indk)*2.0
 endfor

a(0:is,0:is)=r(0:is,0:is)
r(0:is,0:is)=0.0

for k=0,l-1 do begin
indk=(ind+k) mod s
r(0:is2,0:is)=r(0:is2,0:is)+th(k)*a(indk,0:is)
r(s2:is,0:is)=r(s2:is,0:is)+tg(k)*a(indk,0:is)
endfor

a(0:is,0:is)=r(0:is,0:is)
 
s=s/2

; next level
endfor

a2=reform(a,npix)
indi=where(abs(a2) le p, k)
if(k ne 0) then begin
a2(indi)=0.0
print,100.*k/float(npix),' % of wavelet coefficients < ',p
endif else begin
print,'all coeffs < p!'
goto,endprog
endelse

if thresh eq hard  then goto,reconstr

;soft thresholding will better suppress noise artifacts
;hard thresholding will better preserve the visuaal appearance of peaks and jumps

ind=where(abs(a2) gt p)
a2(ind) = a2(ind)/abs(a2(ind))*(abs(a2(ind))-p)

reconstr:

a=reform(a2,n1,n2)

;reconstruction
for e=0,n-1 do begin

is=s-1
s2=s*2
is2=s2-1

r(0:is2,0:is2)=0.0
ind=2*indgen(s)

for k=0,l-1 do begin
indk=(ind+k) mod s2
r(indk,0:is2)=r(indk,0:is2)+a(0:is,0:is2)*h(k)+a(s:is2,0:is2)*g(k)
endfor

a(0:is2,0:is2)=r(0:is2,0:is2)
r(0:is2,0:is2)=0.0

for k=0,l-1 do begin
indk=(ind+k) mod s2
r(0:is2,indk)=r(0:is2,indk)+2.0*a(0:is2,0:is)*h(k)+2.0*a(0:is2,s:is2)*g(k)
endfor

a(0:is2,0:is2)=r(0:is2,0:is2)
s=s*2
; next level
endfor

endprog:

   bestim = bytscl(a)

end