pro dmaxpws, array_out
;    =======
; modified   November, 19, 1998 by Olga Koutchmy

; read maxpows image in jpg format
; plot histogram
; denoising by daub4 wavelet thresholding and cycle averaging
; save denoised image in jpg format

;  maxpows image is read from:
; 'spimaxpa.jpg' for all  images
; 'spimaxpe.jpg' for even images
; 'spimaxpo.jpg' for odd  images
; 'spimaxpf.jpg' for first 7 images
; 'spimaxpl.jpg' for last  8 images

; histogram plot of maxpows image is saved to:
; 'spimhista.ps'  for all  images
; 'spimhiste.ps'  for even images
; 'spimhisto.ps'  for odd  images
; 'spimhistf.jpg' for first 7 images
; 'spimhistl.jpg' for last  8 images

; denoising by daub4 wavelet thresholding and cycle averaging
; denoised image is saved to:
; 'spibesta.jpg' for  all images
; 'spibeste.jpg' for even images
; 'spibesto.jpg' for odd  images
; 'spibestf.jpg' for first 7 images
; 'spibestl.jpg' for last  8 images

; to read: read_jpeg,fname,image,r,g,b; tvlct,r,g,b; tv,image

 kmax =  27
 kkmax=  3
 lmax =  27
 llmax=  3

; cut boundary zeros
 n1   = 512 - (kmax+kkmax+2)
 n2   = 512 - (lmax+llmax+2)
 npix = long(n1)*n2
 npix0= long(512)*512

 choice = 0

 print,'0: even     files'
 print,'1: odd      files'
 print,'2: all      files'
 print,'3: 0 to  6  files'
 print,'4: 7 to 14  files'
 print,'choice = ?'
 read, choice

; read maxpows image
 case choice of
 0:fname= 'spimaxpe.jpg'
 1:fname= 'spimaxpo.jpg'
 2:fname= 'spimaxpa.jpg'
 3:fname= 'spimaxpf.jpg'
 4:fname= 'spimaxpl.jpg'
 else: goto,endprog
 endcase

 read_jpeg, fname, restored_image

 set_plot,'x'

 window,1,xsize=n1,ysize=n2,title='After MAXPOW'
 tvscl, restored_image

 coef = 100./npix
 coeffs = 4.

 w = abs( wtn(congrid(restored_image,512,512,/cubic),coeffs) )
; histogram
 hz =  histogram( reform(w,npix0),min=1.,max=30. )
; free old allocation for
 w  = 0

 window,2,xsize=n1,ysize=n2
 plot, coef* hz, title='Histogram of wavelet coeffs (in %) after MAXPOW',xrange=[1,30]

; save histogram plot of maxpows image
 case choice of
 0:fname= 'spimhiste.ps'
 1:fname= 'spimhisto.ps'
 2:fname= 'spimhista.ps'
 3:fname= 'spimhistf.ps'
 4:fname= 'spimhistl.ps'
 else: goto,endprog
 endcase

 set_plot, 'ps'
 device, filename = fname,xsize=15.,ysize=15.
 plot, coef* hz, title='Histogram of wavelet coeffs (in %) after MAXPOW',xrange=[1,30]
 device, /close
 set_plot, 'x'

; free old allocation for hz
 hz = 0

 pp = 4.
;pp = 7. seems better...
 print,'Threshold=?'
 read,pp

 imager = bytscl( congrid( restored_image, 512, 512, /cubic ))
; free old allocation for
 restored_image = 0
 image0 = bytarr(512,512)

; image is denoised by daub4 wavelet thresholding and cycle averaging
 denoises, pp, imager, image0

 array_out = image0( kmax+1:511-kkmax-1, lmax+1:511-llmax-1)

; free old allocation for image0
 image0 = 0

 window,3,xsize=n1,ysize=n2,title='Denoised image'
 tv, array_out

 mom = moment( array_out ,MDEV=adevi,SDEV=sdevi )
 print,'Denoised image:'
 case choice of
 0:print,'for even files'
 1:print,'for odd  files'
 2:print,'for all  files'
 3:print,'for 0 to 6 files'
 4:print,'for 7 to 14 files'
 else: goto,endprog
 endcase

 print,'Mean                   =', mom(0)
 print,'Variance               =', mom(1)
 print,'Mean Absolute Deviation=', adevi
 print,'Standard Deviation     =', sdevi
 rms = sqrt( total( (array_out - mom(0) )^2 )/npix )
 print,'RMS                    =', rms

; save the denoised image
 case choice of
 0:fname= 'spibeste.jpg'
 1:fname= 'spibesto.jpg'
 2:fname= 'spibesta.jpg'
 3:fname= 'spibestf.jpg'
 4:fname= 'spibestl.jpg'
 else: goto,endprog
 endcase

 write_jpeg,fname,tvrd()

 endprog:
 end