;
; Time-stamp: <2000-02-01 09:30:39 baron>
;
;This example demonstrates the use of IDL's discrete wavelet transform
; and sparse array storage format to compress and store an spectrum 
; (which is simply a 1-dimensional image
; First, a spectrum is read in from an input file
; and transformed into its wavelet representation
; only about 3% of the information is needed to perform the
; inverse wavelet transformation. 



;Finally, the transformed image is reconstructed from the storage file
; and displayed alongside the original.


; Begin by choosing the number of wavelet coefficients to use and a
; threshold value:
; Open thedata file, read a spectrum
; and close the file:
filename='/home/baron/scr/spec/sn87a/sn1987a.2'
filelength,filename,npts
;
readit,filename,npts,wl,flux
; Expand the image to the nearest power of two using cubic
; convolution, and transform the image into its wavelet
; representation using the WTN function:
pwr = 1024*2
flux = CONGRID(flux, pwr, /CUBIC)
coeffs = 12 & thres = 0.1*max(flux)
wtn_flux = WTN(flux, coeffs)
; Write the image to a file using the WRITEU procedure and check
; the size of the file (in bytes) using the FSTAT function:
OPENW, 1, 'original.dat'
WRITEU, 1, wtn_flux
status = FSTAT(1)
CLOSE, 1
PRINT, 'Size of the file is ', status.size, ' bytes.'
;
; now compress the wavelet transformed image
;
wtn_flux_1 = fltarr(N_ELEMENTS(flux))
index_good = where(abs(wtn_flux) ge thres) 
index_bad = where(abs(wtn_flux) lt thres) 
wtn_flux_1[index_good] = wtn_flux[index_good]
wtn_flux_1[index_bad] = 0.
print,n_elements(flux),n_elements(wtn_flux_1)
;
PRINT, 'Percentage of elements under threshold: ',$
   100.*N_ELEMENTS(WHERE(ABS(wtn_flux) LT thres, $
   count)) / N_ELEMENTS(flux)
;
; Apply the inverse wavelet transform to the image:
;
flux_2 = WTN(wtn_flux, COEFFS, /INVERSE)
flux_3 = WTN(wtn_flux_1, COEFFS, /INVERSE)
; Calculate and print the amount of data used in reconstruction of
; the image:
PRINT, 'The image on the right is reconstructed from:', $
   100.0 - (100.* count/N_ELEMENTS(flux)),$
   '% of original image data.'
; Finally, display the original and reconstructed images side by
; side:
WINDOW, 1, $ 
   TITLE = 'Wavelet Image Compression and File I/O'
loadct,14
plot,wl,flux
oplot,wl,flux_3,color=200
oplot,wl,flux_2,color=300