pro thres, a
;   ========
; last modification: November,19,1998
;
; Read an image of spicules (0-14) in .bin format
; Compute and plot the histogram
; of daub4 (Daubechies's wavelet with 2 vanishing moments) wavelet coefficients
; on finest level of decomposition
; (and for smoothed image and details)
; Write the histogram in .jpg and in .ps format
 
   name1 = 'spi_c242'

   print,'Estimation of wavelet threshold using HISTOGRAM'

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

   n1   = 512
   n2   = n1
   npix = long(n1)*n2
   s    = n1
   is   = s-1
   s2   = s/2
   is2  = s2-1

   filter = 'daub4 (Daubechies, with 4 coefficients)'

   print,'Using ', filter, ' filter'

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

   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,*)
   
   set_plot,'x'
   window,3,xsize=512,ysize=512,title='Original image'

   tv, im

; finest level: n = 1
; multiscale decomposition
;
 coeffs = 4
 a = float( wtn( im, coeffs) )
 z = abs( reform (a, npix) )
 ma= max(z)
 print,'Wavelet coefficients on finest level: min = ',min(z),' max = ',ma

 coeff = 100./npix
 hz = coeff*histogram( z, min=0., max=30.)

 set_plot,'x'
 window,0,xsize=512,ysize=512

 st     = 'Image ' + ss
 smax   = string( format='(F8.2)', ma )
 stitre = 'Histogram of Daub(echies)4 wavelet coefficients(0-30-' + smax + ') in %'

 plot, hz, title = stitre, subtitle = st, xrange = [0.,30.]

; D for Daubechies
 file2 = 'Dhisto' + ss + '.jpg'
 write_jpeg, file2, tvrd()

 set_plot,'ps'
 file2 = 'Dhisto' + ss + '.ps'
 device,filename=file2,xsize=15.,ysize=15.

 plot, hz, title = stitre, subtitle = st, xrange = [0.,30.]
 device,/close

 s2s2=long(s2)*s2

 ; smoothed image (resume)
 d0  = abs( a(  0:is2,  0:is2 ) )

 z0 = reform( d0,s2s2 )
 coeff = 100./s2s2
 hz = coeff*histogram( z0, min=0.,max=30.)

 set_plot,'x'
 window,1,xsize=512,ysize=512
 st1 = st + '(Smoothed image)'

 plot, hz, title = stitre, subtitle = st1, xrange = [0.,30.]

;S for Smoothed image
 file2 = 'DhistoS' + ss + '.jpg'
 write_jpeg, file2, tvrd()

 set_plot,'ps'
;S for Smoothed image
 file2 = 'DhistoS' + ss + '.ps'
 device,filename=file2,xsize=15.,ysize=15.

 plot, hz, title = stitre, subtitle = st1, xrange = [0.,30.]
 device,/close

 ; details
 d1  = abs( a( s2:is,  s2:is  ) )
 d2  = abs( a(  0:is2, s2:is  ) )
 d3  = abs( a(  0:is2, s2:is  ) )

 z1  = reform( d1, s2s2 )
 z2  = reform( d2, s2s2 )
 z3  = reform( d3, s2s2 )

 z   = [[z1],[z2],[z3]]

 coeff = 100./(3.*s2s2)
 hz = coeff*histogram( z, min=0.,max=30.)

 set_plot,'x'
 window,2,xsize=512,ysize=512
 st2 = st + '(Details)'

 plot, hz, title = stitre, subtitle = st2, xrange = [0.,30.]

;D for Details
 file2 = 'DhistoD' + ss + '.jpg'
 write_jpeg, file2, tvrd()

 set_plot,'ps'
;D for Details
 file2 = 'DhistoD' + ss + '.ps'
 device,filename=file2,xsize=15.,ysize=15.

 plot, hz, title = stitre, subtitle = st2, xrange = [0.,30.]
 device,/close
 set_plot,'x'
 end