pro spinhisto
;   =========
; last modification: January 25,1999;
; Read an image of spicules (0-14) in .bin format
; Compute and plot the histogram of
; (biortho31,biortho35,biortho44,haar,ortho4,ortho6,orthosym8)
; wavelet coefficients
; on finest level of decomposition
; (and the histograms for smooth image and details, if test ne 'non')
; Compute S.Mallat's threshold estimation using the MEDIAN function
; and put its value in the histogram subtitle
; Write the histogram in .gif and in .ps format
 
   name1 = 'spi_c242'

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

   n1   = 512
   n2   = n1
   pwr   = 512
   pwr1  = pwr - 1
   npix = long(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)

   r = fltarr(n1,n2)

   autre:
;  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,'Estimation of wavelet threshold using HISTOGRAM'
   print,'of ', filter, ' wavelet coeeficients'

   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,*)
   imcut = im( k1:k2, l1:l2 )
   im    = 0
   im    = bytscl( congrid( imcut, pwr, pwr, /cubic ))
   imcut = 0
   
   set_plot,'x'
   window,3,xsize=512,ysize=512,title='Original image'
   tv, im
;-----------------------------------------------------
; finest level: n = 1
; multiscale decomposition
;
 a=float(im)

; zero mean image
 mom = moment( a )
 a   = a - mom(0)

 s = n1

 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)
;----------------------------------------------
 z = abs(reform(a,npix))
 ma= max(z)
 print,'Wavelet coefficients on finest level: min = ',min(z),' max = ',ma

; threshold estimation ( S.Mallat )
 sigma = 3.0* median( z )/0.6745
 print,'Threshold estimation: ',sigma
 ssig  = string( format='(F8.2)', sigma )

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

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

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

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

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

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

 plot, hz, title = stitre, subtitle = st, xrange = [0.,30.]
 device,/close
;----------------------------------------------------------
 test = 'non'
 if test eq 'non' then goto, suite

 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 = filter + '_histoS' + ss + '.gif'
 write_gif, file2, tvrd()

 set_plot,'ps'
;S for Smoothed image
 file2 = filter + '_histoS' + 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 = filter + '_histoD' + ss + '.gif'
 write_gif, file2, tvrd()

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

 plot, hz, title = stitre, subtitle = st2, xrange = [0.,30.]
 device,/close
;-----------------------------------------------------------
  suite:
  set_plot,'x'
  print,'end?(y/n)'
  rep='y'
  read,rep
  if rep eq 'n' then goto, autre
  wdelete,0,1,2,3

 end