pro maxpowspid, array_out
;    ==========
;    November, 19, 1998
;
; function to reconstruct best image using power maximization

; after denoising by daub4 wavelet thresholding

; best image is saved to 'spibest0.jpg' for even images
; best image is saved to 'spibest1.jpg' for odd  images

; 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
;
 im_init = bytarr(640,512)

; 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

; imcut = image( kmax+1:511-kkmax-1, lmax+1:511-llmax-1)

 n1   = 512 - (kmax+kkmax+2)
 n2   = 512 - (lmax+llmax+2)
 npix = long(n1)*n2

 nwindow = 15
 nb      = 8

	;print,' allocate the arrays'

	f =complexarr(n1,n2,/noz)
	ff=complexarr(nb,npix,/noz)
	fsave  = complexarr(npix,/noz)
	fsave2 = complexarr(n1,n2,/noz)
	p=fltarr(nb,npix,/noz)

 print,'even numbers? enter 0'
 print,'odd  numbers? enter 1'
 read,k

        jimage = 0
        ii     = 0

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

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

 image0 = im_init(0:511,*)

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

; window,2,xsize=512,ysize=512,title='Noisy image'
; tv,image

		f        = reform ( fft(image,-1), npix)
		p (ii,*) = abs(f)
		ff(ii,*) = f
        ii = ii + 1
	endfor

	; print,' Now find the maximum power for npix =',npix

        temp = fltarr(ii)
        ii1  = ii-1
        findmax:
        for j = long(0),npix-1 do begin
		temp = p(0:ii1,j)
		pmax = max(temp,loc)
		fsave(j) = ff(loc,j)
        endfor

        ; compute the restored image

	fsave2         = reform(fsave,n1,n2)
	restored_image = fft(fsave2,1)

 set_plot,'x'
 ;window,1,xsize=n1,ysize=n2,title='After MAXPOWSPID'
 ;tvscl, restored_image

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

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

 array_out = image0( kmax+1:511-kkmax-1, lmax+1:511-llmax-1)
 window,0,xsize=n1,ysize=n2,title='Denoised image'
 tv, array_out

; save the denoised image
 if k eq 0 then begin
 write_jpeg,'spibest0.jpg',tvrd()
 endif else begin
 write_jpeg,'spibest1.jpg',tvrd()
 endelse

 return
 endprog:
 end