pro artemis,num,psf,res

; 1. computation and spline-smoothing of the MTF
; 2. construction of the Wiener filter
; 3. restoration
;
; inputs: num = seq.number of the image,
;         psf = array containing the corresponding PSF (from persefona)
; outputs: res =  array containing the restored image (fltarr(1024,1024))
;          0ptionally, res can be stored to a file.
   
on_error,2

mtf=float(fft(psf,-1))*1024     ; computation of the MTF
mtf=mtf(0:512)
 
a=findgen(513)/(0.124*1024)	; arcsec^-1 scale for plotting

; smoothing with the cubic splines (smspline.pro from KIS lib.)

no1=[0,5,10,15,20,30,40,50,75,100,150,200]  ; 1st part of nodes
no2=[250,300,400,512]                       ; 2nd part (default)
nel=4
print,'default set of nodes:'
print,[no1,no2]
goto,exesmoo	
 
newsmoo:
print,'old 2nd part of nodes:',no2
read,'no. of elements in the new 2nd part = ',nel
no2=intarr(nel)
read,'new 2nd part of nodes: ',no2

exesmoo:
tension=1.                 ; control of rigidity of splines
nodes=[no1,no2]
smtf=smspline(mtf,tension,12+nel,nodes)

; visual check of the smoothing

window,0
plot,mtf,psym=3,xra=[-10,520],/xst
oplot,smtf

newdet:
print,'do you want a detail? (y/n)'
if get_kbrd(1) eq 'n' then goto,endsmoo
read,'beginning and end positions: ',i1,i2
i2=i2<512
plot,mtf,psym=1,xra=[i1-5,i2+5],yra=[min(mtf(i1:i2)),max(mtf(i1:i2))],/xst,/yst
oplot,smtf
oplot,intarr(600)
goto,newdet

endsmoo:
print,'smoothing ok? (y/n)
if get_kbrd(1) eq 'n' then goto,newsmoo

smtf=smtf/smtf(0)      ; renormalization

; definition of the cutoff frequency...
mi=min(smtf,pmi)             ; looking for minimum
print,'minimum value of mtf = ',mi,' at position ',pmi

if mi gt 0. then begin
   print,'The minimum is positive. Shall we define it as the cutoff? (y/n)'
   if get_kbrd(1) eq 'n' then stop
   smtf(pmi+1:512)=0.
   cros0=pmi+1
endif else begin
   cros0=min(where(smtf le 0))
   smtf(cros0:512)=0.
endelse

cr0=1/(cros0/(0.124*1024))
print,'cutoff position at',cros0,' pixels, i.e.',cr0,' arcsec'

; reading the input image
; !!!! TO BE CHANGED WITH RESPECT TO WHERE THE IMAGES ARE STORED !!!!
stra=string('/scratch/sob/')
strb=string('imc.')
; !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   ima=fltarr(1315,1035)
   dcn=strtrim(num,2)
   if strlen(dcn) eq 3 then str='0'+dcn
   if strlen(dcn) eq 2 then str='00'+dcn
   if strlen(dcn) eq 1 then str='000'+dcn
   print,'Now reading image:'
   print,strb+str
   openr,1,stra+strb+str
   readu,1,ima
   close,1

; extracting 1024x1024
ima=ima(291:1314,6:1029)

; apodization...

; subtraction of the mean value
mea=mean(ima)
ima=ima-mea

; creating and applying an apodization window
perc=5.
np=rfix(1024*perc/100.)
u=fltarr(1024)+1.
u(0)=0.
for i=1,np do begin
	u(i)=(1.-cos(!pi*i/np))/2.
endfor
u(1024-np:1023)=reverse(u(1:np))
ima=ima*(u#u)

;adding the mean
ima=ima+mea

; 2-D Fourier transform of the original image
fima=fft(ima,-1)

; --------------- creation of the smoothed Moon limb ----------------

psfs=shift(psf,256)

mol=fltarr(512)
for i=0,511 do mol(i)=total(psfs(0:i))   ; integration of psf = limb
mol=[mol,reverse(mol)]

fmol=fft(mol,-1)

heavis=fltarr(512)			; Heaviside step function
heavis(256:511)=1.

; ---------------------- filtering ----------------------------------

; filter selection

window,4,xsize=550,ysize=700
plot,a,smtf,linestyle=5,xra=[-0.2,4.2],yra=[-0.2,6],/xst,/yst

newc:
wdelete,0,2,3

read,'Wiener filter parameter C= ',c
wfilt=(c+1)*smtf/(c*smtf^2+1)

wset,4
plot,a,wfilt,xra=[-0.2,4.2],yra=[-0.2,6],/xst,/yst,/noerase

fmols=fmol*[wfilt,reverse(wfilt(1:511))]
rmol=float(fft(fmols,1))

window,6,xsize=600,ysize=700
plot,heavis,linestyle=1,xra=[206,306],yra=[-0.1,1.1],/xst,/yst
oplot,mol,linestyle=1
oplot,rmol

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
	res=0
	stop
endif

; --------- transition to a 2-dim Wiener filter ---------------

wfi2=fltarr(1024,1024)
f1=fltarr(513,513)
print,'creating a 2D Wiener filter'

cros0=min(where(wfilt le 0))		; zero crossing point of wfilt,
					; 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)=wfilt(r1)+(wfilt(r1+1)-wfilt(r1))*(r-r1)
		endelse
	endfor
endfor

f2=reverse(f1(1:511,*),1)
f4=reverse(f1(*,1:511),2)
f3=reverse(f2(*,1:511),2)

wfi2(0:512,0:512)=f1
wfi2(513:1023,0:512)=f2
wfi2(513:1023,513:1023)=f3
wfi2(0:512,513:1023)=f4
f1=0
f2=0
f3=0
f4=0

print,'filter ready, now doing the Fourier restoration'
;tvscl,rebin(wfi2,512,512)

; ----------- filtering and inverse Fourier transform -----------

res=float(fft(fima*wfi2,1))
wfi2=0

; visualization of the restored image

rre=rebin(res,512,512)
arcs=findgen(400)*0.062

selbox:
print,'select (with a cursor) center of the zoom (left), right to skip'
window,0,xsize=512,ysize=512
tvscl,rre
cursor,x,y,/device,/down
if !err eq 4 then goto,fin

  x=99>x*2<923
  y=99>y*2<923

rbox=rebin(res(x-99:x+100,y-99:y+100),400,400)
ibox=rebin(ima(x-99:x+100,y-99:y+100),400,400)
print,'rms contrast of unrestored field =',stdev(ibox,me)/me
print,'rms contrast of restored field   =',stdev(rbox,me)/me

window,2,xsize=800,ysize=400
tvscl,ibox,0   ; original image
tvscl,rbox,1   ; restored image
window,3,xsize=800,ysize=300,ypos=200

; plot of a selected row
selrow:
wset,2
print,'select (with a cursor) the row to be plotted (left), right to skip'
cursor,x,y,/device,/down
if !err eq 4 then goto,endsel
wset,3
plot,arcs,rbox(*,y),xtitle='arcsec',yra=[min(rbox(*,y)),max(rbox(*,y))],/yst
oplot,arcs,ibox(*,y),color=160
goto,selrow
endsel: goto,selbox

fin:
print,'Restoration finished.'
print,'will you try another filter? (y/n)'
if get_kbrd(1) eq 'y' then goto,newc

print,'do you want to save the restored image in a file? (y/n)'
if get_kbrd(1) eq 'y' then begin
   strb=string('imr.')
   dcn=strtrim(num,2)
   if strlen(dcn) eq 3 then str='0'+dcn
   if strlen(dcn) eq 2 then str='00'+dcn
   if strlen(dcn) eq 1 then str='000'+dcn
   print,'saving image:'
   print,strb+str
   openw,1,stra+strb+str
   writeu,1,res
   close,1	
endif

end