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,'nbpix=',nbpix


   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


; suppression des zeros aux bords

for i=0,n2-1 do begin
i1=i
v=where( im_init(*,i) eq 0B, count) 
if( count ne n1) then goto,sui1
endfor  

sui1:
for i=n2-1,0,-1 do begin
i2=i
v=where( im_init(*,i) eq 0B, count) 
if( count ne n1) then goto,sui2
endfor

sui2:
for i=0,n1-1 do begin
j1=i
v=where( im_init(i,*) eq 0B, count) 
if( count ne n2) then goto,sui3
endfor

sui3:
for i=n1-1,0,-1 do begin
j2=i
v=where( im_init(i,*) eq 0B, count) 
if( count ne n2) then goto,sui4
endfor

sui4:
print,'suppression des zeros aux bords :',i1,i2,j1,j2

   imsz= im_init(j1:j2,i1:i2)
 
   
help,imsz

;prolongement par approximation cubique

dim1=j2-j1+1
print,'dim1=',dim1
;on decoupe pwr pixels en x
if(dim1 gt pwr) then imc=imsz(0:pwr1,*)
help,imc
dim2=i2-i1+1
print,'dim2=',dim2
if((dim1 ne pwr) or (dim2 ne pwr)) then imc=congrid(imc,pwr,pwr,/cubic)

im=imc

help,im

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


if(ouinon eq 'oui') then goto,saut2
   window, 1, xsize=512, ysize=512,title='Initial image.'
   tvscl, im 

saut2:
; tests

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

print,'thres=?'
read,thres


;m=fltarr(pwr,pwr)
;m=replicate(0.0,pwr,pwr)

ns=8
ns1=ns-1
ns2=float(ns*ns)

;for sh2=0,ns1 do begin
;for sh1=0,ns1 do begin

image_1=im

;image_1=shift(im,sh1,sh2)

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.2*(wmax-wmin)
print,'binsize=',binsize
y=reform(wtn_image,nbpix)
wma=binsize/2.
wmi=-wma
print,'wmi=',wmi,' wma=',wma
print,'binsize?'
read,binsize
hw=histogram(y,min=wmi,max=wma,binsize=binsize,reverse_indices=r)
help,hw
lhw=n_elements(hw)
if(lhw gt 10) then lhw=10
for i=0,lhw-1 do begin
print,i,hw(i)
endfor

hwmax=max(hw,loc)
print,'histogram:',hwmax,loc
thres=loc*binsize+wmi
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)

;m=m+shift(image_2,-sh1,-sh2)/ns2

;endfor
;endfor

suite:
m=image_2
print,100.*n_elements(where(abs(wtn_image) 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,2,xsize=512,ysize=512,title='Denoised image.'
tvscl,m

saut3:
endfor

bestim=m

iminit=im
end