pro maxpowspinh, array_out
;    ===========
;    November, 19, 1998
;
; version with minimized virtual memory

; denoising each of 15 images by Haar wavelet thresholding and cycle averaging
; before reconstructing best image using power maximization in Fourier domain
; computing Mean,Variance,Mean Absolute Deviation,Standard Deviation
; it is saved to 'moment_h_maxpa.dat'
; "maxpowed" image is saved to:
; 'h_spimaxpa.jpg' for all  images
; 'h_spimaxpe.jpg' for even images
; 'h_spimaxpo.jpg' for odd  images
; 'h_spimaxpf.jpg' for first 7 images
; 'h_spimaxpl.jpg' for last  8 images

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

; based on algorithm proposed by Serge Koutchmy
; Written feb 25, 1990 by S. l. Keil
; Modified 27 Feb 1990 using suggestions by Phil Wiborg
; Modified  Aug 1990  by Olga Koutchmy
; Modified  Sept 1990  by T.A. Darvann
;
; name1 = '/home/gaberle2/koutchmy/olga/spicules/spi_c242'
 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

; cut boundary zeros
 n1   = 512 - (kmax+kkmax+2)
 n2   = 512 - (lmax+llmax+2)
 npix = long(n1)*n2
 npix0= long(512)*512
 k1 = kmax+1
 k2 = 511-kkmax-1
 l1 = lmax+1
 l2 = 511-llmax-1

 nwindow = 15

 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

 case choice of
 0:begin
   nb = 8  & k = 0 & kk = 2 & last = 14
   end
 1:begin
   nb = 7  & k = 1 & kk = 2 & last = 13
   end
 2:begin
   nb = 15 & k = 0 & kk = 1 & last = 14
   end
 3:begin
   nb = 7  & k = 0 & kk = 1 & last = 6
   end
 4:begin
   nb = 8  & k = 7 & kk = 1 & last = 14
   end
 else: goto,endprog
 endcase
       
 print,'si test=0,alors on ne traite que 2 images! test=?'
 read,test
 if test eq 0 then begin
 nb = 2  & k = 0 & kk = 1 & last = 1
 endif

 pp = 7.
; =====
 print,'Threshold=?'
 read,pp

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

        ii     = 0

;       for nb images' Fourier transform
	ff=complexarr(nb,npix,/noz)

 set_plot,'x'

	for i = k, last, kk do begin
;       ============================
 number = string( format='(2I)', i)
 no     = strcompress( number,/remove_all)
 s      = no
 if (i lt 10) then s = '0' + no
;binary format
 file = name1 + s + '.bin'

 if test eq 0 then begin
 print,'Read image (640,512) from file ', file
 endif

 openr, 1, file
 readu, 1, im_init
 close,1

 if test eq 0 then help

 image0  = im_init( 0:511,*)

 if test eq 0 then begin
 window,0,xsize=512,ysize=512,title='noisy image'
 tvscl, image0
 endif

; image is denoised by Haar wavelet thresholding and cycle averaging
 imager = bytscl(congrid(image0( k1:k2,l1:l2),512, 512,/cubic ))

 if test eq 0 then begin
 window,1,xsize=512,ysize=512,title='congrid noisy image'
 tvscl, imager
 endif

 denoisec, pp, imager, image0
 imager = 0

 if test eq 0 then begin
 window,2,xsize=512,ysize=512,title='denoised image'
 tvscl, image0
 endif

 ff(ii,*)= reform ( fft(image0( k1:k2, l1:l2),-1), npix)
 image0 = 0
; next file to read
        ii = ii + 1
	endfor

	fsave  = complexarr(npix,/noz)

	; print,' Now find the maximum power for npix =',npix
        temp = fltarr(ii)
        ii1  = ii-1

        for j = long(0),npix-1 do begin
	temp = abs( ff(0:ii1,j) )
	pmax = max(temp,loc)
;       loc - image's number
;       j   - pixel's number
	fsave(j) = ff(loc,j)
;       next pixel
        endfor

; free old allocation for ff
        ff = 0

        ; compute the restored image
	restored_image = fft(reform(fsave,n1,n2) ,1)

 if test eq 0 then help

; free old allocation for fsave
        fsave = 0

 array_out = bytscl( congrid( restored_image, 512, 512, /cubic ))

 window,3,xsize=512,ysize=512,title='maxpowed image (all denoised images)'
 tvscl, array_out                                       

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

 openw, 2, 'moment_h_maxpa.dat'

 case choice of
 0:printf,2,'  "maxpowed" image ( for even denoised files ):'
 1:printf,2,'  "maxpowed" image ( for odd  denoised files ):'
 2:printf,2,'  "maxpowed" image ( for all  denoised files ):'
 3:printf,2,'  "maxpowed" image ( for 0 to 6 denoised files ):'
 4:printf,2,'  "maxpowed" image ( for 7 to 14 denoised files ):'
 else: goto,endprog
 endcase
 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

; save image
 case choice of
 0:fname= 'h_spimaxpe.jpg'
 1:fname= 'h_spimaxpo.jpg'
 2:fname= 'h_spimaxpa.jpg'
 3:fname= 'h_spimaxpf.jpg'
 4:fname= 'h_spimaxpl.jpg'
 else: goto,endprog
 endcase

 write_jpeg,fname,tvrd()

 endprog:
 end