Viewing contents of file '../idllib/contrib/buie/clean.pro'
;+
; NAME:
; clean
; PURPOSE:
; Remove a PSF from an image via the ``clean'' algorithm.
; DESCRIPTION:
; CATEGORY:
; CCD data processing
; CALLING SEQUENCE:
; clean,image,psf,xloc,yloc,maxdist,iters,new_image,resid $
; DISPLAY=display,VERBOSE=verbose
; INPUTS:
; image - Original source image to be cleaned.
; psf - PSF image at same sampling resolution as image.
; xloc - X location of "object"
; yloc - Y location of "object"
; maxdist - Maximum distance from xyloc to look for local max.
; iters - Number of cleaning iterations to perform.
; OPTIONAL INPUT PARAMETERS:
; KEYWORD INPUT PARAMETERS:
; CROSS - If set to a number greater than or equal to 0, the peak at
; each iteration is found by using a cross-correlation
; of the psf against the image rather than just using raw DN.
; If 0, then the full size of the psf is used. If 1 or greater,
; a sub-region from the center that is 2*CROSS+1 pixels square
; is used as the convolution kernal.
; DISPLAY - Display intermediate results, if 0, this is suppressed,
; if non-zero, this is the interval for the display, that is,
; DISPLAY=10 would cause a display every 10th iteration.
; GAIN - "Gain" of the clean process, the default value is 0.05 and
; is the scaled amount of the psf removed at each step.
; VERBOSE - Verbose printout of intermediate steps to the screen. Just
; like display, VERBOSE=0 suppresses output, VERBOSE=n will
; print information every nth iteration.
; OUTPUTS:
; new_image - Clean-ed image result.
; resid - Remains of the original image after clean-ed image is removed.
; KEYWORD OUTPUT PARAMETERS:
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
; 1/26/94, written by Marc W. Buie, Lowell Observatory, algorithmic insight
; graciously provided by Tod Lauer (NOAO, Tucson).
; 2/25/94, MWB, added GAIN and CROSS keywords.
; 5/19/94, MWB, changed psf normalization from peak to total. Should
; now be strictly flux conserving.
; 5/21/94, MWB, Changed CROSS to allow for only a portion of psf for
; convolution.
;-
pro clean,image,psf,xloc,yloc,maxdist,iters,new_image,resid, $
CROSS=cross,DISPLAY=display,GAIN=gain,VERBOSE=verbose
self='CLEAN: '
if badpar(image, [1,2,3,4,5],2,caller=self+'(image) ') then return
if badpar(psf, [1,2,3,4,5],2,caller=self+'(psf) ') then return
if badpar(xloc, [1,2,3], 0,caller=self+'(xloc) ') then return
if badpar(yloc, [1,2,3], 0,caller=self+'(yloc) ') then return
if badpar(maxdist,[1,2,3], 0,caller=self+'(maxdist) ') then return
if badpar(iters, [1,2,3], 0,caller=self+'(iters) ') then return
if badpar(display,[0,1,2,3], 0,caller=self+'(display) ',default=0) then return
if badpar(verbose,[0,1,2,3], 0,caller=self+'(verbose) ',default=0) then return
if badpar(gain, [0,4,5], 0,caller=self+'(GAIN) ',default=0.05) then return
if badpar(cross, [0,1,2,3], 0,caller=self+'(CROSS) ',default=-1) then return
; Set up information on the input image and make a working copy.
work = image
s=size(work)
imw = s[1]
imh = s[2]
; Set up the new image
new_image = fltarr(imw,imh)
; Set up information on the psf array.
npsf = psf / max(psf)
psfnorm = total(npsf)
s=size(psf)
psfw = s[1]
psfh = s[2]
boxm,npsf,psfw/2,psfh/2,psfw/2,psfh/2,psfx,psfy
if cross eq 0 then begin
kernal = psf
if (verbose gt 0) or (psfx ne psfw-psfx) or (psfy ne psfh-psfy) then begin
print,'Complete psf used for cross correlation.'
if psfx ne psfw-psfx then $
print,'Warning: Cross correlation not symmetric in X'
if psfy ne psfh-psfy then $
print,'Warning: Cross correlation not symmetric in Y'
endif
endif else if cross gt 0 then begin
dw = min([cross,psfx,psfy,psfw-psfx,psfh-psfy])
x1 = psfx - dw
x2 = psfx + dw
y1 = psfy - dw
y2 = psfy + dw
kernal = psf[x1:x2,y1:y2]
if (verbose gt 0) or (dw ne cross) then begin
print,'PSF core used for cross correlation.'
if dw ne cross then $
print,'Warning, cross too big, reduced from ',cross,' to ',dw $
else $
print,'Half-width of core used is ',dw
endif
endif
for i=0L,iters do begin
if cross gt 0 then begin
workc = convol(work,kernal)
boxm,workc,xloc,yloc,maxdist,maxdist,xpos,ypos,/ABSMAX
endif else begin
boxm,work,xloc,yloc,maxdist,maxdist,xpos,ypos,/ABSMAX
endelse
peak = work[xpos,ypos]
scale_fac = peak * gain
; Increment the output image.
new_image[xpos,ypos] = new_image[xpos,ypos] + scale_fac*psfnorm
; Determine the sub-array for subtraction
il = max([0,xpos-psfx])
pl = max([0,psfx-xpos])
ir = min([imw-1,xpos+(psfw-1-psfx)])
pr = min([psfw-1,psfx+(imw-1-xpos)])
ib = max([0,ypos-psfy])
pb = max([0,psfy-ypos])
it = min([imh-1,ypos+(psfh-1-psfy)])
pt = min([psfh-1,psfy+(imh-1-ypos)])
; Subtract 5% at the max location
work[il:ir,ib:it] = work[il:ir,ib:it] - npsf[pl:pr,pb:pt]*scale_fac
if keyword_set(verbose) then begin
if i mod verbose eq 0 then begin
print,'#',i,' x,y ',xpos,ypos,' sf ',scale_fac,' (',il,':',ir,',', $
ib,':',it,') (',pl,':',pr,',',pb,',',pt,')', $
peak,work[xpos,ypos], $
format='(a,i6.6,a,i3,1x,i3,a,f6.1,a,i3.3,a,i3.3,a,' + $
'i3.3,a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,i3.3,a,1x,f6.1,1x,f6.1)'
endif
endif
; Optional display
if display ne 0 then begin
if i mod display eq 0 then begin
mn = min(work[il:ir,ib:it])
mx = max(work[il:ir,ib:it])
sz=20 ; 35 ; 20
zf= 6 ; 4 ; 6
new_wid = (2*sz+1)*zf
tmp1 = work[xloc-sz:xloc+sz,yloc-sz:yloc+sz]
tmp2 = new_image[xloc-sz:xloc+sz,yloc-sz:yloc+sz]
if !d.x_size lt new_wid or !d.y_size lt new_wid then $
setwin,0,xsize=new_wid*2,ysize=new_wid
ypos = (!d.y_size - new_wid)/2
xpos1 = (!d.x_size)/2 - new_wid
xpos2 = (!d.x_size)/2
tvscl,rebin(tmp1,new_wid,new_wid,/sample),xpos1,ypos
tvscl,rebin(tmp2,new_wid,new_wid,/sample),xpos2,ypos
endif
endif
endfor
resid = work
end