pro maxpowspis, array_out
;    ==========
;    November, 19, 1998
;
; version with minimized virtual memory
; function to reconstruct best image using power maximization

; denoising "maxpowed" image by daub4 wavelet thresholding and cycle averaging
; computing Mean,Variance,Mean Absolute Deviation,Standard Deviation
; it is on 'moment_maxpa_d.dat'
; denoised best image is saved to:
; 'spimaxpa_d.jpeg' for all  images
; 'spimaxpe_d.jpeg' for even images
; 'spimaxpo_d.jpeg' for odd  images
; 'spimaxpf_d.jpeg' for first 7 images
; 'spimaxpl_d.jpeg' 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

 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
       
        ii     = 0

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

	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'

 im_init = bytarr(640,512)

; print,'Read image (640,512) from file ', file
 openr, 1, file
 readu, 1, im_init
 close,1

 ff(ii,*)= reform ( fft(im_init( kmax+1:511-kkmax-1, lmax+1:511-llmax-1),-1), npix)

; free old allocation for im_init
 im_init=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)

; free old allocation for fsave
        fsave = 0

 set_plot,'x'

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

; image is denoised by daub4 wavelet thresholding and cycle averaging
 imager = bytscl( congrid( restored_image, 512, 512, /cubic ))
 denoises, pp, imager, image0
 imager = 0

 array_out = bytscl(congrid(image0( kmax+1:511-kkmax-1,lmax+1:511-llmax-1),512 512,/cubic ))
 image0 = 0
 window,0,xsize=512,ysize=512,title='Denoised "maxpowed" image'
 tv, array_out

 mom = moment(image0(kmax+1:511-kkmax-1,lmax+1:511-llmax-1),MDEV=adevi,SDEV=sdevi)

 openw, 2, 'moment_maxpa_d.dat'

 case choice of
 0:printf,2,'Denoised "maxpowed" image ( for even files ):'
 1:printf,2,'Denoised "maxpowed" image ( for odd  files ):'
 2:printf,2,'Denoised "maxpowed" image ( for all  files ):'
 3:printf,2,'Denoised "maxpowed" image ( for 0 to 6 files ):'
 4:printf,2,'Denoised "maxpowed" image ( for 7 to 14 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= 'spimaxpe_d.jpeg'
 1:fname= 'spimaxpo_d.jpeg'
 2:fname= 'spimaxpa_d.jpeg'
 3:fname= 'spimaxpf_d.jpeg'
 4:fname= 'spimaxpl_d.jpeg'
 else: goto,endprog
 endcase

 write_jpeg,fname,tvrd()

 endprog:
 end