pro maxpowspi, array_out
;   ====================
; Modified   April, 10, 1998

; Att! version without memory optimization !

; function to reconstruct best image using power maximization
; for 15 images of spicules
; each image is denoised by denoise.pro
; (Haar filter, without cycle averaging)

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

; images non recentrees a Sunspot
name1='/home/serge/spicules/spi242'

; images recentrees
name1='spi_c242'

;if name1 eq 'spi_c242' then testfor0, kmax,kkmax,lmax,llmax

if name1 eq 'spi_c242' then begin
kmax = 27
kkmax=  3
lmax = 27
llmax=  3
endif else begin
kmax =-1
kkmax=-1
lmax =-1
llmax=-1
endelse

n1 = 512
n2 = n1
npix = long(n1)*n2

nwindow = 15

	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)

        jimage = 0

ii = 0

thresh='hard'
print,'What thresholding?(hard/soft)'
read,thresh

print,'Threshold=?'
read,pp

 print,'Number of levels=?'
 read, n

	for i = 0, nwindow-1 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'

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

print,'Read image (640,512) from file ', file

image = im_init(0:511,*)

        window,0,xsize=512,ysize=512,title='Before denoising ( Haar filter )'
	tv,image

imcut = image( kmax+1:511-kkmax-1, lmax+1:511-llmax-1)
image = congrid( imcut, 512, 512, /cubic )
;testzeros,image,kmax,kkmax,lmax,llmax

 imcut = image
 denoise,thresh,pp,n,n1,imcut,image

; image is denoised by wavelet thresholding

                window,2,xsize=512,ysize=512,title='After denoising ( Haar filter )' 
		tv,image

		f=reform ( fft(image,-1), npix)

		p(ii,*)=abs(f)

		ff(ii,*)=f
ii = ii + 1
	endfor

; now find the maximum power at each spatial frequency point for the nwindow arrays

	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)
	imager         = bytscl(restored_image)

        window,1,xsize=512,ysize=512,title='After MAXPOW'
	tv,imager
	
; return the answers

   finishup:
        array_out = imager
	return

   endprog:
end