function convolve, image, psf, FT_PSF=psf_FT, FT_IMAGE=imFT, NO_FT=noft, $
CORRELATE=correlate, AUTO_CORRELATION=auto, $
NO_PAD = no_pad
;+
; NAME:
; CONVOLVE
; PURPOSE:
; Convolution of an image with a Point Spread Function (PSF)
; EXPLANATION:
; The default is to compute the convolution using a product of
; Fourier transforms (for speed).
;
; The image is padded with zeros so that a large PSF does not
; overlap one edge of the image with the opposite edge of the image.
;
; CALLING SEQUENCE:
;
; imconv = convolve( image1, psf, FT_PSF = psf_FT )
; or:
; correl = convolve( image1, image2, /CORREL )
; or:
; correl = convolve( image, /AUTO )
;
; INPUTS:
; image = 2-D array (matrix) to be convolved with psf
; psf = the Point Spread Function, (size < or = to size of image).
;
; The PSF *must* be symmetric about the point
; FLOOR((n_elements-1)/2), where n_elements is the number of
; elements in each dimension. For Gaussian PSFs, the maximum
; of the PSF must occur in this pixel (otherwise the convolution
; will shift everything in the image).
;
; OPTIONAL INPUT KEYWORDS:
;
; FT_PSF = passes out/in the Fourier transform of the PSF,
; (so that it can be re-used the next time function is called).
; FT_IMAGE = passes out/in the Fourier transform of image.
;
; /CORRELATE uses the conjugate of the Fourier transform of PSF,
; to compute the cross-correlation of image and PSF,
; (equivalent to IDL function convol() with NO rotation of PSF)
;
; /AUTO_CORR computes the auto-correlation function of image using FFT.
;
; /NO_FT overrides the use of FFT, using IDL function convol() instead.
; (then PSF is rotated by 180 degrees to give same result)
;
; /NO_PAD - if set, then do not pad the image to avoid edge effects.
; This will improve memory and speed of the computation at the
; expense of edge effects. This was the default method prior
; to October 2009
; METHOD:
; When using FFT, PSF is centered & expanded to size of image.
; HISTORY:
; written, Frank Varosi, NASA/GSFC 1992.
; Appropriate precision type for result depending on input image
; Markus Hundertmark February 2006
; Fix the bug causing the recomputation of FFT(psf) and/or FFT(image)
; Sergey Koposov December 2006
; Fix the centering bug
; Kyle Penner October 2009
; Add /No_PAD keyword for better speed and memory usage when edge effects
; are not important. W. Landsman March 2010
;-
compile_opt idl2
sp = size( psf_FT,/str ) & sif = size( imFT, /str )
sim = size( image )
if (sim[0] NE 2) OR keyword_set( noft ) then begin
if keyword_set( auto ) then begin
message,"auto-correlation only for images with FFT",/INF
return, image
endif else if keyword_set( correlate ) then $
return, convol( image, psf ) $
else return, convol( image, rotate( psf, 2 ) )
endif
if keyword_Set(No_Pad) then begin
sc = sim/2 & npix = N_elements( image )
if (sif.N_dimensions NE 2) OR ((sif.type NE 6) AND (sif.type NE 9)) OR $
(sif.dimensions[0] NE sim[1]) OR (sif.dimensions[1] NE sim[2]) then imFT = FFT( image,-1 )
if keyword_set( auto ) then $
return, shift( npix*real_part(FFT( imFT*conj( imFT ),1 )), sc[1],sc[2] )
if (sp.N_dimensions NE 2) OR ((sp.type NE 6) AND (sp.type NE 9)) OR $
(sp.dimensions[0] NE sim[1]) OR (sp.dimensions[1] NE sim[2]) then begin
sp = size( psf )
if (sp[0] NE 2) then begin
message,"must supply PSF matrix (2nd arg.)",/INFO
return, image
endif
Loc = ( sc - sp/2 ) > 0 ;center PSF in new array,
s = (sp/2 - sc) > 0 ;handle all cases: smaller or bigger
L = (s + sim-1) < (sp-1)
psf_FT = conj(image)*0 ;initialise with correct size+type according
;to logic of conj and set values to 0 (type of psf_FT is conserved)
psf_FT[ Loc[1], Loc[2] ] = psf[ s[1]:L[1], s[2]:L[2] ]
psf_FT = FFT( psf_FT, -1, /OVERWRITE )
endif
if keyword_set( correlate ) then $
conv = npix * real_part(FFT( imFT * conj( psf_FT ), 1 )) $
else conv = npix * real_part(FFT( imFT * psf_FT, 1 ))
sc = sc + (sim MOD 2) ;shift correction for odd size images.
return, shift( conv, sc[1], sc[2] )
endif else begin
sc = floor((sim-1)/2) & npix = n_elements(image)*4.
; the spooky factor of 4 in npix is because we're going to pad the image
; with zeros
if (sif.N_dimensions NE 2) OR ((sif.type NE 6) AND (sif.type NE 9)) OR $
(sif.dimensions[0] NE sim[1]) OR (sif.dimensions[1] NE sim[2]) then begin
; here is where we make an array with twice the dimensions of image and
; pad with zeros -- thanks to Daniel Eisenstein for this fix
image_big = dblarr(sim[1]*2,sim[2]*2)
image_big[0:sim[1]-1,0:sim[2]-1] = image
imFT = FFT( image_big,-1 )
npix = n_elements(image_big)
endif
if keyword_set( auto ) then begin
intermed = shift( npix*real_part(FFT( imFT*conj( imFT ),1 )), sc[1],sc[2] )
return, intermed[0:sim[1]-1,0:sim[2]-1]
endif
if (sp.N_dimensions NE 2) OR ((sp.type NE 6) AND (sp.type NE 9)) OR $
(sp.dimensions[0] NE sim[1]) OR (sp.dimensions[1] NE sim[2]) then begin
sp = size( psf )
if (sp[0] NE 2) then begin
message,"must supply PSF matrix (2nd arg.)",/INFO
return, image
endif
; this obfuscated line determines the offset between the center of the
; image and the center of the PSF
Loc = ( sc - floor((sp-1)/2) ) > 0
psf_image = dblarr(sim[1]*2,sim[2]*2)
psf_image[Loc[1]:Loc[1]+sp[1]-1, Loc[2]:Loc[2]+sp[2]-1] = psf
psf_FT = FFT(psf_image, -1)
endif
if keyword_set( correlate ) then begin
conv = npix * real_part(FFT( imFT * conj( psf_FT ), 1 ))
conv = shift(conv, sc[1], sc[2])
endif else begin
conv = npix * real_part(FFT( imFT * psf_FT, 1 ))
conv = shift(conv, -sc[1], -sc[2])
endelse
return, conv[0:sim[1]-1,0:sim[2]-1]
endelse
end
Return to the shapelets web page or the code help menu.
| Last modified on 2nd Mar 2009 by Richard Massey. |