pro spind, bestim,iminit,w

ouinon='oui'
print,'execution sans visualisation?(oui/non)' 
read,ouinon
if(ouinon eq 'oui') then goto,saut1

   set_plot,'x'

saut1:
n1=640
n2=512
pwr=512

   im_init   = bytarr(n1,n2)
   im        = bytarr(pwr,pwr)


pwr1=pwr-1
nbpix = 512L*pwr

   print,'number of files to process ?'
   read, nwindow

   nb = nwindow - 1


   name1 = 'spi242'

   for i = 0, nb do begin

   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,0:511)


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

print,'coeffs(4,12,20)=?' 
read,coeffs

print,'thres=?'
read,thres

image_1=im

suit:
wtn_image=wtn(image_1,coeffs)

w = wtn_image

wmin=min(wtn_image)
wmax=max(wtn_image)
print,'wmin=',wmin,' wmax=',wmax
binsize=0.01*(wmax-wmin)
print,'binsize=1% (wmax-wmin)',binsize
y=reform(wtn_image,nbpix)
print,'binsize?'
read,binsize
hw=histogram(y,min=wmin,max=wmax,binsize=binsize,reverse_indices=r)
help,hw
if ouinon eq 'non' then plot, hw,title='histogram'  
lhw=n_elements(hw)
if(lhw gt 100) then lhw=100
for i=0,lhw-1 do begin
print,i,hw(i)
endfor

hwmax=max(hw,loc)
print,'histogram:',hwmax,loc
thres=abs(loc*binsize+wmin)
print,r(loc+1)-r(loc),'elements ds bin loc'
print,'thres calcule=',thres

;openw,1,'original.dat'
;writeu,1,wtn_image
;status=fstat(1)
;close,1
;print,'the size of the file is ',status.size,' bytes.'

sprs_image=sprsin(wtn_image,thres=thres)

;write_spr,sprs_image,'sparse.dat'

;openr,1,'sparse.dat'
;status=fstat(1)
;close,1
;print,'the size of the file is ',status.size,' bytes.'

;sprs_image=read_spr('sparse.dat')

wtn_image=fulstr(sprs_image)

image_2=wtn(wtn_image,coeffs,/inverse)


suite:
m=image_2
print,100.*n_elements(where(abs(w) lt thres,count))/n_elements(image_1)
print,'The image is reconstructed from:'
print,100.-(100.*count/n_elements(image_1)),'% of original image data.'

if(ouinon eq 'oui') then goto,saut3
window,1,xsize=pwr*2,ysize=pwr,title='Initial and denoised images'
tvscl,im,0,0
tvscl,m,pwr1,0

saut3:
endfor

bestim=m

iminit=im
end