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