; XFPR : XWINDOW oriented Fractional Pixel Restoration Program ; Routine to REMOVE Cosmic Ray Hits pro fixcrh,d,dx,dy,wdebug ; Remove Cosmic Ray hits ; Set a pixel that exceed 8 standard deviations above ; the local average to the average value. widget_control,wdebug, set_value= "Removing Cosmic Ray Hits" for i=0,dx-1 do begin for j=0,dy-1 do begin thres = d(*,*,0,0)*0 replace = thres*0 for k=0,dx-1 do begin for l=0,dy-1 do begin if((i ne k) or (j ne l)) then begin thres = thres + d(*,*,k,l) replace = replace + d(*,*,k,l) endif endfor endfor thres = thres/float(dx*dx-1.0) replace = replace/float(dx*dx-1.0) var = thres*0 for k=0,dx-1 do begin for l=0,dy-1 do begin if((i ne k) or (j ne l)) then begin var = var + (thres-d(*,*,k,l))^2 endif endfor endfor var = thres + 8.0*sqrt(var) hits = where(var lt d(*,*,i,j),count) tmp=thres*0 tmp(hits) = 1 tmp = d(*,*,i,j) tmp(hits)= replace(hits) d(*,*,i,j)=tmp endfor endfor return end ; Routine for Removing bad or dead pixels pro fixpix,a ; Set a bad pixel to ones that exceed 2.5 standard deviations below ; local average. a=float(a) s=smooth(a,9) v=smooth(a^2,9) v=sqrt(v-s^2) ind = where(s-2.5*v gt a) ; check and make sure that there are some bad pixels s=size(ind) if(s(0) eq 0) then return s=smooth(a,3) a(ind) = (9.0*s(ind)-a(ind))/8.0 a=fix(a) return end ;Richardson-Lucy Algorithm pro rl,d,ah,p,niter,wdebug eta = 1e-6 ah=d p = fft(p,1) for i = 1,niter do begin widget_control,wdebug, set_value= $ "Preforming RL Iteration #" + string(i,format='(i3)') tvscl,abs(ah) ah = fft(ah,1) dh = conj(p)*ah dh = fft(dh,-1) dh = d/(dh+eta) dh = fft(dh,1) dh = p*dh dh = fft(dh,-1) ah = fft(ah,-1)*dh endfor p = fft(p,-1) end ; FPR algorithm pro fpr,data,output,psf,sizex,sizey,dx,dy,niter,avg,crh,rmpsf,wdebug FALSE = 0 TRUE = 1 window,2,xs=dx*sizex,ys=dy*sizey if(crh) then fixcrh,data,dx,dy,wdebug ; co-add data based on the parameter avg output = output*0.0 for i=0,dy-1 do begin for j=0,dx-1 do begin if(avg) then begin output = output + $ shift(rebin(data(*,*,j,i),sizex*dx,sizey*dy,sample=1),j,i) endif else begin for k=0,sizey-1 do begin for l=0,sizex-1 do output(dx*l+j,dy*k+i)=data(l,k,j,i) endfor endelse endfor endfor tvscl,output ; created distoration filter indx = (indgen(dx*sizex)-dx*sizex/2)/float(dx*sizex/2) indy = (indgen(dy*sizex)-dy*sizey/2)/float(dy*sizey/2) maskx = 1.0+complex(cos(!pi*indx),sin(!pi*indx)) for i=2,dx-1 do begin maskx = maskx+complex(cos(float(i)*!pi*indx),sin(float(i)*!pi*indx)) endfor masky = 1.0+complex(cos(!pi*indy),sin(!pi*indy)) for i=2,dx-1 do begin masky = masky+complex(cos(float(i)*!pi*indy),sin(float(i)*!pi*indy)) endfor mask = maskx#masky if(avg) then begin mask = mask*conj(mask) endif ; add psf to distoration filter if known mask = shift(mask,dx*sizex/2,dy*sizey/2) if(rmpsf) then begin psf=fft(shift(psf,dx*sizex/2,dy*sizey/2),1) mask = mask*psf endif ; normalize and uncenter distortion filter mask = mask/total(mask) mask = fft(mask,-1) tmp = output ; apply the RL algorithm if(niter ne 0) then begin rl,tmp,output,mask,niter,wdebug output = abs(output) endif return end ; XFPR EVENT HANDLER pro xfpr_event, ev common parameters, xnum,ynum,ftype,avg,crh,sx,sy,rli,oname,rmpsf,$ files,mkout,pname,lfn common data, output common zoom, sizex,sizey,dx,dy,zscl common debug, wdebug bell = bytarr(2) & bell(0)=7 rmbadpix = 0 widget_control, ev.id, get_uvalue=value if (n_elements(value) eq 0) then value = '' name = strmid(tag_names(ev, /structure_name), 7, 4) case (name) of ; Action buttons control case statement "BUTT": begin case (value) of ; EXIT button Procedure "DONE": begin WIDGET_CONTROL, /destroy, ev.top return endcase ; Zoom out procedure rebins input without interpolation ; decrease magnification by an integer factor until the ; desplay size equals that of the data "ZOUT" : begin if((zscl gt 1) and mkout) then begin zscl = zscl-1 window,2,xs=dx*sizex*zscl,ys=dy*sizey*zscl tvscl,rebin(output,dx*sizex*zscl,dy*sizey*zscl,sample=1) endif return endcase ; Zoom in procedure rebins and displays the input without ; interpolation, increases magnification by an integer amount ; in both the x and y direction until the display size is ; > 512x512 "ZIN" : begin if(mkout) then begin if((zscl*dy*sizey lt 512) and zscl*dy*sizey lt 512) then begin zscl = zscl+1 window,2,xs=dx*sizex*zscl,ys=dy*sizey*zscl tvscl,rebin(output,dx*sizex*zscl,dy*sizey*zscl,sample=1) endif endif return endcase ; Load File Table button procedure "LOAD" : begin ; Write out starting message WIDGET_CONTROL, wdebug , $ set_value="Loading Optional File Table" ; read input file ; set jump tp ioerr (unable to open file) on_ioerror,ioerr ; get file name widget_control,lfn,get_value=value openr,1,value(0) ;read file names on_ioerror,ioerr2 name = string(80) for i=0,1+ynum*2 do begin for j=0,1+xnum*2 do begin readf,1,name widget_control,files(j,i),set_value=name endfor endfor close,1 WIDGET_CONTROL, wdebug , set_value='' return endcase ; Co-add button procedure "CA" : begin WIDGET_CONTROL, wdebug , $ set_value="Preforming Direct Co-Addition" niter = 0 goto,rla endcase ; Richardson - Lucy button procedure "RL" : begin ; Write out starting message widget_control,rli,get_value=value niter = fix(value(0)) ; read input files depending on file type rla: if(ftype eq 0) then begin ; Read FITS input files ; reads first fits to determine sizex and sizey ; open file to see if it exist and check if its ; simple FITS file ; set jump to ioerr if one can not open file on_ioerror,ioerr ; get file name widget_control,files(0,0),get_value=value ; open file and check if it is a FITS file ; portion of code was lifted from readfits openr,1,value(0) hdr = bytarr(80) readu,1,hdr header = string( hdr > 32b ) if strmid( header(0), 0, 8) NE 'SIMPLE ' then begin goto, ioerr1 endif close,1 ; Read first FITS file d=readfits(value(0),silent) if(rmbadpix) then begin widget_control,wdebug, $ set_value="Fixing bad pixels in " + string(value(0)) fixpix,d endif tmp=size(d) dx = (xnum+1)*2 & dy = (ynum+1)*2 sizex = tmp(1) & sizey = tmp(2) ; set and check input file x and y sizes if( sizex ne sizey) then begin print,string(bell) widget_control,wdebug, $ set_value= "Unable to deal with non square data sets" return endif input = fltarr(sizex,sizey,dx,dy) input(*,*,0,0) = float(d) widget_control, sx, set_value=string(sizex,format='(i4)') widget_control, sy, set_value=string(sizey,format='(i4)') ; read remaining FITS files for i=0,1+ynum*2 do begin for j=0,1+xnum*2 do if (i ne 0 or j ne 0) then begin ; set jump to ioerr and check fits files on_ioerror,ioerr widget_control,files(j,i),get_value=value openr,1,value(0) hdr = bytarr(80) readu,1,hdr header = string( hdr > 32b ) if strmid( header(0), 0, 8) NE 'SIMPLE ' then begin goto, ioerr1 endif close,1 ; read FITS files d=readfits(value(0)) if(rmbadpix) then begin widget_control,wdebug, $ set_value="Fixing bad pixels in " $ + string(value(0)) fixpix,d endif input(*,*,j,i) = float(d) endif endfor endif else begin ; READ raw input files ; set jump tp ioerr (unable to open file) on_ioerror,ioerr dx = (xnum+1)*2 & dy = (ynum+1)*2 ; get x and y file sizes set by users widget_control,sx,get_value=value sizex = fix(value(0)) widget_control,sy,get_value=value sizey = fix(value(0)) ; check to make sure the data is a square array if( sizex ne sizey) then begin print,string(bell) widget_control,wdebug, set_value= $ "Unable to deal with non square data sets" return endif ; check file size if( sizex lt 3 or sizey lt 3) then begin print,string(bell) widget_control,wdebug, set_value = $ "X or Y input size is to SMALL" return endif ; allocate data space input = fltarr(sizex,sizey,dx,dy) d = intarr(sizex,sizey) ; read raw input images for i=0,1+ynum*2 do begin for j=0,1+xnum*2 do begin on_ioerror,ioerr widget_control,files(j,i),get_value=value openr,1,value(0) on_ioerror,ioerr2 readu,1,d close,1 if(rmbadpix) then begin widget_control,wdebug, $ set_value="Fixing bad pixels in " + $ string(value(0)) fixpix,d endif input(*,*,j,i) = float(d) endfor endfor endelse ; open display window window,2,xs=sizex*dx,ys=sizey*dy zscl = 1 ; allocate space for restoration output = fltarr(sizex*dx,sizey*dy) output = output*0 mkout = 1 ; TRUE ; READ PSF : XFPR assumes that the PSF is the same size as ; the RESTORATION i.e. size = (sizex*dx,sizey*dy) ; The routine also assumes that the psf is centered in the ; middle of the array if(rmpsf) then begin if(ftype eq 0) then begin ; READ FITS PSF : This routine will adjust PSF array size on_ioerror,ioerr widget_control,pname,get_value=value openr,1,value(0) hdr = bytarr(80) readu,1,hdr header = string( hdr > 32b ) if strmid( header(0), 0, 8) NE 'SIMPLE ' then begin goto, ioerr1 endif close,1 d=readfits(value(0),silent) tmp=size(d) if(tmp(1) gt sizex*dx) or (tmp(2) gt sizey*dx) then begin print,string(bell) widget_control,wdebug, set_value = $ "ERROR : PSF is Larger than data" return endif ; Convert PSF to reals and center it psf=fltarr(sizex*dx,sizey*dy) psf(0:tmp(1)-1,0:tmp(2)-1) = float(d) psf = shift(psf,(tmp(1)+sizex*dx)/2,(tmp(2)+sizey*dy)/2) psf = psf/total(psf) endif else begin ; READ raw PSF file : This routine will only read ; one size = (sizex*dx,sizey*dy) on_ioerror,ioerr d = intarr(sizex*dx,sizey*dy) on_ioerror,ioerr widget_control,pname,get_value=value openr,1,value(0) on_ioerror,ioerr2 readu,1,d close,1 psf = float(d) endelse endif widget_control,wdebug, $ set_value="Executing the FPR Algorithm " fpr, input, output, psf, sizex, sizey, dx, dy, niter, $ avg, crh, rmpsf, wdebug WIDGET_CONTROL, wdebug , set_value='' return ioerr: ; Error procedure if file can not be opened close,1 print,string(bell) WIDGET_CONTROL, wdebug , set_value= $ 'ERROR : Unable to open ' + string(value(0)) return ioerr1: ; Error procedure if file is not a SIMPLE FITS file close,1 print,string(bell) WIDGET_CONTROL, wdebug , set_value= $ 'ERROR : ' + string(value(0)) + ' is NOT a SIMPLE FITS File' return ioerr2: ; Error procedure if the input file can not be read close,1 print,string(bell) WIDGET_CONTROL, wdebug , set_value = $ 'ERROR : Unable to read ' + string(value(0)) return endcase "WRITE": begin ; WRITE button procedure if(mkout eq 0) then begin print,string(bell) WIDGET_CONTROL, wdebug , set_value = $ 'ERROR: No output to Write' return endif WIDGET_CONTROL, wdebug , $ set_value='Writing Display to a File' WIDGET_CONTROL, oname, get_value=name ; scale output into integer format if(max(output) gt (2^15-1)) then begin output = (2^15-1)*output/max(output) endif output = fix(output) if(ftype eq 0) then begin writefits,name(0),output endif else begin openw,1,name(0) writeu,1,output close,1 endelse endcase "CT": begin ; Modify Color Table button procedure if (mkout) then xloadct return endcase else: endcase WIDGET_CONTROL, wdebug , set_value='' endcase "TEXT": begin ; General text widget procedure handles all text input boxes ; Does Nothing at the moment return endcase else: begin ; Procedure to handle all toggle button widgets and their ; associated parameters case (string(value)) of ; Average toggle buttons "AVG": begin if(avg eq 0) then avg = 1 else avg = 0 endcase ; Cosmic Ray Hit Removal toggle buttons "CRH": begin if(crh eq 0) then crh = 1 else crh = 0 endcase ; Remove PSF toggle buttons "PSF": begin if(rmpsf eq 0) then rmpsf = 1 else rmpsf = 0 endcase ; Number of samples in the x direction toggle buttons "XSIZE": begin if(xnum eq 0) then xnum = 1 else xnum = 0 endcase ; Number of samples in the y direction toggle buttons "YSIZE": begin if(ynum eq 0) then ynum = 1 else ynum = 0 endcase ; File type toggle buttons "FTYPE": begin WIDGET_CONTROL, oname, get_value=name if(ftype eq 0) then begin ftype = 1 ; Set default name to xfpr.fits if no name entered if(name(0) eq "xfpr.fits") then begin WIDGET_CONTROL, oname, set_value="xfpr.raw" endif endif else begin ftype = 0 ; Set default name to xfpr.raw if no name entered if(name(0) eq "xfpr.raw") then begin WIDGET_CONTROL, oname, set_value="xfpr.fits" endif endelse endcase else: endcase endcase endcase end ; main procedure xfpr pro xfpr common parameters common data common zoom common debug ; define space for widget pointers and initialize states of ; toggle buttons files = lonarr(4,4) xnum=0 & ynum=0 & ftype=0 & avg = 0 & crh = 0 & mkout = 0 & rmpsf = 0 zscl = 0 ; create base widget base = WIDGET_BASE(title='XFPR', /column) b1 = WIDGET_BASE(base, /frame, /column) t1 = WIDGET_BASE(b1,/row) ; label the data file widgets t2 = widget_label(t1,value = 'Data File Names',/align_center,xsize=720) t1 = WIDGET_BASE(b1,/row) t2 = widget_label(t1,value = ' ',xsize=40,/align_center) for i=0,3 do begin str = "x = " + string(i,format='(i1)') t2 = widget_label(t1, xsize = 161,value = str,/align_center) endfor ; create the data file widgets for i = 0,3 do begin t1 = WIDGET_BASE(b1,/row) str = "y = " + string(i,format='(i1)') t2 = widget_label(t1, value = str,xsize=40,/align_center) for j = 0,3 do begin files(j,i) = widget_text(t1, /editable, xsize=24, ysize=1) endfor endfor ; create and label the toggle widget for the number of samples in x t1 = WIDGET_BASE(b1,/row) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='# of Observations in X:', $ /EXCLUSIVE, /NO_REL, ['Two', 'Four'], UVALUE='XSIZE', SET_VALUE=xnum) t2 = widget_label(t1, $ value ='(Number of observations acquire in the x direction)') ; create and label the toggle widget for the number of samples in y t1 = WIDGET_BASE(b1,/row) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='# of Observations in Y:', $ /EXCLUSIVE, /NO_REL, ['Two', 'Four'], UVALUE='YSIZE', SET_VALUE=ynum) t2 = widget_label(t1, $ value ='(Number of observations acquire in the y direction)') ; create and label the file type toggle widget t1 = WIDGET_BASE(b1,/row) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='File Types:', /EXCLUSIVE, /NO_REL, $ ['Fits', 'Raw'], UVALUE='FTYPE', SET_VALUE=ftype) t2 = widget_label(t1, $ value ='(The Format Used by the Files Listed Above and Output)') ; create and label the input x and y size text widget t1 = WIDGET_BASE(b1,/row) t2 = widget_label(t1, value = 'Size of Files in Pixel : X') sx = widget_text(t1, /editable, xsize=4, ysize=1) t2 = widget_label(t1, value = 'Y') sy = widget_text(t1, /editable, xsize=4, ysize=1) t2 = widget_label(t1, value = ' Name of Optional File List :') lfn = widget_text(t1, /editable, xsize=36, ysize=1) ; create and label the Pixel averaging toggle widget t1 = WIDGET_BASE(b1,/row) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='Use Pixel Averaging:', $ /EXCLUSIVE, /NO_REL, ['No', 'Yes'], UVALUE='AVG', SET_VALUE=avg) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='Remove Cosmic Ray Hits:', $ /EXCLUSIVE, /NO_REL, ['No', 'Yes'], UVALUE='CRH', SET_VALUE=crh) ; create and label the remove PSF toggle widget and PSF name text widget t1 = WIDGET_BASE(b1,/row) t2 = CW_BGROUP(t1, /ROW, LABEL_LEFT='Remove PSF:', $ /EXCLUSIVE, /NO_REL, ['No', 'Yes'], UVALUE='PSF', SET_VALUE=rmpsf) t2 = widget_label(t1, value = '',xsize = 46) t2 = widget_label(t1, value = 'PSF File Name :') pname = widget_text(t1, /editable, xsize=60, ysize=1) ; create and label the number of RL iteration text widget t1 = WIDGET_BASE(b1,/row) t2 = widget_label(t1, value = 'Number of RL Iterations :') rli = widget_text(t1, /editable, xsize=4, ysize=1,value='25') t2 = widget_label(t1, value = 'Output File Name :') if(ftype eq 0) then begin oname = widget_text(t1, /editable, xsize=60, ysize=1,value="xfpr.fts") endif else begin oname = widget_text(t1, /editable, xsize=60, ysize=1,value="xfpr.raw") endelse ; create and label the error message box text widget t1 = WIDGET_BASE(b1,/row) t2 = widget_label(t1, value = 'Message Box:') wdebug = WIDGET_TEXT(t1, xsize=100, ysize=1,value='Starting XFPR') ; create the action button widgets b2 = WIDGET_BASE(base, /frame, /row,/align_center) t1 = widget_button(b2, value = 'Load Image List', uvalue = 'LOAD') t1 = widget_button(b2, value = 'Direct Co-add', uvalue = 'CA') t1 = widget_button(b2, value = 'FPR', uvalue = 'RL') t1 = widget_button(b2, value="Zoom Out", uvalue = "ZOUT") t1 = widget_button(b2, value="Zoom In", uvalue = "ZIN") t1 = widget_button(b2, value="Modify Color Table", uvalue = "CT") t1 = widget_button(b2, value="Write Output", uvalue = "WRITE") t1 = widget_button(b2, value="Exit", uvalue = "DONE") ; manage the base widget widget_control, base, /realize XMANAGER, 'xfpr', base end