PRO restele,in,out

; RESTORATION of an array IN (2-D frame or 3-D movie) for
; the 2-D theoretical MTF of a diffraction-limited telescope.
; Restored array = OUT.
; The IN array must have xdim = ydim = even number.
; calling:   RESTELE,in,out
; Note: program asks for Wiener filter parameter C. C=20 seems to be adequate.

;on_error,2

; PARAMETERS OF RESTORATION ---------------------------------------------

	stp=0.1658  ;0.06218   ; SCALE (PIXEL SIZE) OF THE IMAGE  (arcsec/pix)
	lambda=1.5542d-4  ;5.250d-5   ; WORKING WAVELENGTH (cm)
	d=44.  ;50.		  ; DIAMETER OF THE TELESCOPE (cm)

; array analysis --------------------------------------------------------

siz=size(in)
if siz(0) lt 2 then goto,wrong
if siz(1) ne siz(2) then goto,wrong
if (siz(1) mod 2) ne 0 then goto,wrong     ; checking even dimensions
if siz(0) eq 3 then nim=siz(3) else nim=1  ; number of frames in array
n=siz(1)
a=findgen(n/2+1)/(n*stp)                   ; scale array for plotting
print,'Input array - size and no. of frames:',n,nim

; 1-D calculation of MTF ------------------------------------------------

mtf=dblarr(n)
c=lambda*2.062648d5/(d*n*stp)   ; 206264.8 arcsec in radian
ss=0

   for i=1,n do begin
	      
      omega=c*(i-1)
      if omega lt 1.d0 then begin
        tt=acos(omega)-omega*sqrt(1.-omega*omega)
      endif else begin
        ss=ss+1
	if ss eq 1 then print,'MTF cutoff: ',i-1
	tt=0.d0
      endelse
      mtf(i-1)=2.*tt/3.141592654d0

   endfor

window,2
plot,a,mtf,linestyle=5,xra=[-0.2,5],yra=[-0.2,5],/xst,/yst

; creating the 1-D Wiener filter ----------------------------------------

newc:
read,'Wiener filter parameter C =',c
t1=(c+1)*mtf/(c*mtf*mtf+1)

wset,2
plot,a,t1,xra=[-0.2,5],yra=[-0.2,5],/xst,/yst,/noerase
print,'Change filter? (y/n)'
if get_kbrd(1) eq 'y' then goto,newc

print,'Do you want to stop here? (y/n)'
if get_kbrd(1) eq 'y' then begin
	out=0
	goto,fin
endif

t1=float(t1)

; transition to 2-D Wiener filter ---------------------------------------

t=fltarr(n,n)
f1=fltarr(n/2+1,n/2+1)

cros0=min(where(t1 le 0))		; zero crossing point of t1,
					; i.e. we can compute only
for j=0,cros0 do begin			; the non-zero part. The rest
	for i=0,cros0 do begin		; are zeros by definition.
		ii=float(i)
		jj=float(j)
		r=sqrt(ii*ii+jj*jj)
		r1=fix(r)
		if r1 ge cros0 then begin
		  f1(i,j)=0.
		endif else begin
		  f1(i,j)=t1(r1)+(t1(r1+1)-t1(r1))*(r-r1)
		endelse
	endfor
endfor

f2=reverse(f1(1:n/2-1,*),1)

t(0:n/2,0:n/2)=f1			; t - 2-D Wiener filter
t(n/2+1:n-1,0:n/2)=f2
t(n/2+1:n-1,n/2+1:n-1)=reverse(f2(*,1:n/2-1),2)
t(0:n/2,n/2+1:n-1)=reverse(f1(*,1:n/2-1),2)
f1=0
f2=0

; RESTORATION ------------------------------------------------------------

out=in

;apodization window

   perc=5.		  ; percentage of frame size to be damped
   np=nint(n*perc/100.)
   u=fltarr(n)+1.
   u(0)=0.
   for i=1,np do u(i)=(1.-cos(!pi*i/np))/2.
   u(n-np:n-1)=reverse(u(1:np))
   u=u#u		  	; 2-D apodization window

for k=0,nim-1 do begin    ; loop over frames

   if nim eq 1 then ima=in else ima=in(*,*,k)

; apodization

;   mea=mean(ima)
;   ima=ima-mea  
;   ima=ima*u
;   ima=ima+mea		  ; apodized frame

; 2-D Fourier filtering
   
   f1=fft(ima,-1)
   f1=float(fft(f1*t,1))

   if nim eq 1 then out=f1 else out(0,0,k)=f1

endfor

goto,fin

wrong:
print,'*** ERROR: WRONG INPUT ARRAY *** STOPPED'
out=0
goto,fin

fin:
print,'End of program. Filter parameter C =',c
END