pro chop_rest,filein,fileout,filedat,k,dis,iter
;
;
;
;
;NAME:
; CHOP_REST
;
;PURPOSE:
; To restore chopped and nodded images in infrared astronomy by means
; of the projected Landweber method (M.Bertero, P.Boccacci,
; F.Di Benedetto, M.Robberto, 1999, "Restoration of chopped and nodded
; images in infrared astronomy",Inverse Problems, 15, 345-372).
;
;CALLING SEQUENCE:
; CHOP_REST,filein,fileout,filedat,k,[dis,iter]
;
;INPUT:
; filein - file FITS containing the chopped and nodded image to be
; restored. The procedure will determine automatically the
; dimensions of the image: N=number of rows, Nc=number of
; columns.
; k - number of pixels corresponding to the chopping amplitude;
; the procedure assumes k<N.
; dis - minimum value of the relative discrepancy; this value must
; be of the order of the relative r.m.s error on the chopped
; and nodded image; if it is not supplied, it is set equal to
; the default value (0.08).
; iter - maximum number of iterations; if this value is not
; supplied, the number of iterations is set equal to 0 and the
; discrepancy criterion is used. If both optional parameters are
; inserted (dis,iter), then the routine CHOP_REST stops as soon
; as one of the two criteria is satisfied.
;
;OUTPUT:
; fileout - file FITS containing the restored image. The size of the
; output image ((N+2k)*Nc) is determined automatically by
; the procedure.
; filedat - file containing information about the recostruction: start
; file (filein), chopping amplitude (k), maximum number of
; iterations (iter), minimum value of the relative discrepancy
; (dis), two columns list containing the current number of
; iterations and the corresponding value of the relative
; discrepancy.
;
;
;EXAMPLES:
; CHOP_REST,'input.fits','output.fits','info.dat',k,dis,iter
;
;PROCEDURES USED:
; FITS_READ,FITS_WRITE (astronomical library
; http://hires.gsfc.nasa.gov/stis/docs/fitsdoc/fitsdoc.html)
;
;HISTORY:
; Written in C by: P.Boccacci (Department of Physics, University of
; Genova, Italy)
; Converted to IDL V5.1 by: B.Carrara and P.Zerega (Department of
; Computer and Information Science,
; University of Genova, Italy) March 1999
if N_params() LT 4 then begin ;Enough parameters supplied?
print, $
'Syntax - CHOPPING,filein,fileout,filedat,k,[dis,iter] '
return
endif
;set default's parameters
if N_params() EQ 4 then begin
iter=0
dis=0.08
endif
if N_params() EQ 5 then begin
iter=0
endif
td=1.
td1=1.
tau=1.8/16.
;read image to restore
fits_read,filein,img
s=size(img)
Nc=s[1]
N=s[2]
openw, 1,filedat
;write in file for information (filedat)
printf,1,"Start File: ",filein
printf,1,"Chopping Amplitude: ",k
IF (iter eq 0) THEN BEGIN
n1=0
td1=1
endif else begin
printf,1,"Max number of iteration: ",iter
n1=1
endelse
;write on the file for information (filedat)
if(td1) then printf,1,"Value of discrepancy: ",dis
printf,1
printf,1," n. iteration discrepancy "
M=N+2*k
A=make_array(M,N,/float,value=0.)
sol=make_array(Nc,M,/float,value=0.)
ATAF=make_array(Nc,M,/float)
for i=0,N-1 do begin
A[i,i]=-1.
A[i+k,i]=2.
A[i+2*k,i]=-1.
endfor
norm2data=sqrt(total(img^2))
AT=transpose(A)
ik=0
ATg=AT##img
;Landweber loop
while (((td GT dis) OR (td1 eq 0)) AND ((ik LT iter) OR (n1 eq 0))) do begin
;columns loop
ATAF[*,0:k-1]=sol[*,0:k-1]-2*sol[*,k:2*k-1]+sol[*,2*k:3*k-1]
ATAF[*,N+k:M-1]=sol[*,N-k:N-1]-2*sol[*,N:N+k-1]+sol[*,N+k:M-1]
if (2*k GT N) then begin
ATAF[*,k:N-1]= -2*sol[*,0:N-1-k]+5*sol[*,k:N-1] $
-4*sol[*,2*k:N-1+k]+sol[*,3*k:M-1]
ATAF[*,N:2*k-1]= -2*sol[*,N-k:k-1]+4*sol[*,N:2*k-1] $
-2*sol[*,N+k:3*k-1]
ATAF[*,2*k:k+N-1]= sol[*,0:N-1-k]-4*sol[*,k:N-1]+5*sol[*,2*k:k+N-1] $
-2*sol[*,3*k:M-1]
endif else begin
if (2*k EQ N) then begin
ATAF[*,k:2*k-1]= -2*sol[*,0:k-1]+5*sol[*,k:2*k-1] $
-4*sol[*,2*k:3*k-1]+sol[*,3*k:4*k-1]
ATAF[*,N:N-1+k]= sol[*,N-2*k:N-1-k]-4*sol[*,N-k:N-1] $
+5*sol[*,N:N-1+k]-2*sol[*,N+k:M-1]
endif else begin
ATAF[*,k:2*k-1]= -2*sol[*,0:k-1]+5*sol[*,k:2*k-1] $
-4*sol[*,2*k:3*k-1]+sol[*,3*k:4*k-1]
ATAF[*,2*k:N-1]= sol[*,0:N-2*k-1]-4*sol[*,k:N-k-1] $
+6*sol[*,2*k:N-1]-4*sol[*,3*k:N-1+k]+sol[*,4*k:M-1]
ATAF[*,N:N-1+k]= sol[*,N-2*k:N-1-k]-4*sol[*,N-k:N-1] $
+5*sol[*,N:N-1+k]-2*sol[*,N+k:M-1]
endelse
endelse
sol=sol+tau*(ATg-ATAF)
;positivity
for j=0,M-1 do for i=0,Nc-1 do if (sol[i,j] LT 0.) then sol[i,j]=0.
;total discrepancy
td=sqrt(total(((A##sol)-img)^2))/norm2data
ik=ik+1
printf,1,ik," ",td
print,ik, td
endwhile
;end of Landweber loop
close,1
;write the solution on the FITS file (fileout)
fits_write,fileout,sol
print, 'Recostruction successfully written'
end