pro spin
;    ==== 

;   Modified  November, 19, 1998

    print,'Denoising of 15 images of spicules...'
    print,'by hard thresholding and cycle averaging'
    print,'using Haar filter'

    set_plot,'x'

    n1    = 640
    n2    = 512
    pwr   = 512
    pwr1  = pwr - 1

    im_init   = bytarr( n1, n2  )

    print,'number of files to process = ?'
    read, nwindow
;         nwindow = 15

;   recentered images
;   -----------------
    name1 = 'spi_c242'

;   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)
    m1   = pwr - (kmax+kkmax+2)
    m2   = pwr - (lmax+llmax+2)

    imav = replicate( 0., m1, m2 )

    openw, 2, 'moments_h.dat'

    print,'threshold = ?'
    read, p
;         p = 7.

    nb = nwindow - 1

   ;Boucle sur les images
   ;---------------------
   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 (', n1, ',', n2, ') from file ', file

   im = im_init(0:pwr1,*)
   
   titre = 'Image (recentered) ' + no
   window,0,xsize=pwr,ysize=pwr,title=titre
   tv, im 

   imcut  = im( k1:k2, l1:l2 )
   im     = 0

;  denoising and cycle averaging
;  -----------------------------

   im    = bytscl( congrid( imcut, pwr, pwr, /cubic ))
   imcut = 0   

   print,'denoising and cycle averaging...'
   denoisec, p, im, bestim
   im = 0

   titre = 'Denoised image (recentered) ' + no
   window,1,xsize=pwr,ysize=pwr,title=titre
   tv, bestim

   file = name1 + ss + 'h.jpg'
   write_jpeg, file, tvrd()

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

 printf,2,'      '
 printf,2,'Denoised Image (recentered) ', i
 printf,2,'==========================='
 printf,2,'Mean                   =', mom(0)
 printf,2,'Variance               =', mom(1)
 printf,2,'Mean Absolute Deviation=', adevi
 printf,2,'Standard Deviation     =', sdevi

;  average image
   imav   = imav + bestim(k1:k2,l1:l2)
   bestim = 0

   ;next file
   endfor

   imav = imav/nwindow

   window,2,xsize=pwr,ysize=pwr,title='Average denoised image(recentered)'
   tvscl, congrid(imav,pwr,pwr,/cubic)

   file = name1  + 'hav.jpg
   write_jpeg, file, tvrd()

 mom = moment( imav ,MDEV=adevi,SDEV=sdevi )

 printf,2,'      '
 printf,2,'Average image (all denoised images (recentered))'
 printf,2,'================================================'
 printf,2,'Mean                   =', mom(0)
 printf,2,'Variance               =', mom(1)
 printf,2,'Mean Absolute Deviation=', adevi
 printf,2,'Standard Deviation     =', sdevi
 close, 2

 imav = 0

 endprog:
 end