pro thresh
;   ======
;  modified: January 11, 1999

   print,'Denoising by wavelet soft or hard global thresholding' 
   print,'( without cyclic averaging )...'
;using a filter from (biortho31,biortho35,biortho44,qbiortho44,haar,ortho4,ortho6,orthosym8)

   set_plot,'x'

   im_init   = bytarr(640,512)
   im        = bytarr(512,512)

   n1 = 512
   n2 = n1
   pwr = 512
   pwr1 = pwr - 1
   npix = long(n1)*n2

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

;  test for boundary zeros
;  if name1 eq 'spi_c242' then testfor0, kmax,kkmax,lmax,llmax
   kmax =  27
   kkmax=  3  
   lmax =  27
   llmax=  3

   k1   = kmax+1
   k2   = pwr1-kkmax-1
   l1   = lmax+1
   l2   = pwr1-llmax-1
;  imcut = image( kmax+1:pwr1-kkmax-1, lmax+1:pwr1-llmax-1)

;  filtres
   filter = 'haar'
print, 'Filter name?(biortho31,biortho35,biortho44,qbiortho44,haar,ortho4,ortho6,orthosym8)?'
   read, filter

   openr, 3, filter
   readf, 3, l

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

   readf, 3, h, th
   close, 3

   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

   print,'File number (0-14) to process ?'
   read, i

   name1 = 'spi_c242'

   number = string( format='(2I)', i)
   no     = strcompress( number,/remove_all)
   ss =  no
   if (i lt 10) then ss = '0' + no

   file = name1 + ss + '.bin'

   openr, 1, file
   readu, 1, im_init
   close,1

   print, 'Read unformatted initial image (640,512) from file ', file

   im = im_init(0:511,*)
   imcut = im( k1:k2, l1:l2 )
   im    = bytscl( congrid( imcut, pwr, pwr, /cubic ))
   imcut = 0
   
   titre = 'Original Image ' + no
   window, 0, xsize=512, ysize=512, title=titre
   tv,im

   filename = 'moments_dImage' + no
   openw, 4, filename
   printf, 4, 'Image ' + no

   mom = moment( im( k1:k2, l1:l2 ),MDEV=adevi,SDEV=sdevi )

   printf,4,' '
   printf,4,'Mean                   =', mom(0)
   printf,4,'Variance               =', mom(1)
   printf,4,'Mean Absolute Deviation=', adevi
   printf,4,'Standard Deviation     =', sdevi

;--------------------------------------------------------
 soft   = 'soft'
 hard   = 'hard'
 thresh = soft

 autre:
 print,'what thresholding?(hard/soft)'
 read, thresh

 print,'threshold = ?'
 read, p

 print,'number of levels = ?'
 read,n

;--------------------------------------
; multiscale decomposition
 s = n1
 a = float(im)
 mom = moment( a )
 a = a - mom(0)

 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)

if( thresh eq soft ) then  p = 1.5*p

;hard thresholding will better preserve the visual appearance of peaks and jumps
indi= where(abs(a2) le p, k)
if(k ne 0) then begin
a2(indi) = 0.0
ltp = 100.*k/float(npix)
print,ltp,' % 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
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

   gtp = 100. - ltp
   sgtp = strcompress(string( format='(F8.2)', gtp ),/remove_all) 
   sp = strcompress(string( format='(F8.2)', p ),/remove_all) 
   titre = 'Image '+no+'; '+filter+' filter; '+thresh+' thresh='+sp+'; '+sgtp+'%'

   a = a + mom(0)
   window, 1, xsize=512, ysize=512, title= titre
   tvscl, a

   file2 = 'Denoised_Image' + ss + '.gif'
   write_gif, file2, tvrd()

   mom = moment( a( k1:k2, l1:l2 ),MDEV=adevi,SDEV=sdevi )

   printf, 4,  'Noise removed using '+filter + ' filter; ' + 'thresholding =' + sp
   printf, 4, 'Image is reconstructed from ' + sgtp +'% of original data'

   printf,4,' '
   printf,4,'Mean                   =', mom(0)
   printf,4,'Variance               =', mom(1)
   printf,4,'Mean Absolute Deviation=', adevi
   printf,4,'Standard Deviation     =', sdevi
   close, 4

   print,'end?(y/n)'
   rep='y'
   read,rep
   if rep eq 'n' then goto, autre

endprog:
   wdelete,0,1
end