/
convolve.pro
84 lines (77 loc) · 2.78 KB
/
convolve.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
function convolve, image, psf, FT_PSF=psf_FT, FT_IMAGE=imFT, NO_FT=noft, $
CORRELATE=correlate, AUTO_CORRELATION=auto
;+
; 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).
;
; 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).
;
; 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)
; METHOD:
; When using FFT, PSF is centered & expanded to size of image.
; HISTORY:
; written, Frank Varosi, NASA/GSFC 1992.
; Converted to IDL V5.0 W. Landsman September 1997
;-
sp = size( psf_FT ) & sif = size( imFT )
sim = size( image ) & sc = sim/2 & npix = N_elements( 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 (sif[0] NE 2) OR (sif[sif[0]+1] NE 6) OR $
(sif[1] NE sim[1]) OR (sif[2] NE sim[2]) then imFT = FFT( image,-1 )
if keyword_set( auto ) then $
return, shift( npix*float( FFT( imFT*conj( imFT ),1 ) ), sc[1],sc[2] )
if (sp[0] NE 2) OR (sp[sp[0]+1] NE 6) OR $
(sp[1] NE sim[1]) OR (sp[2] 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 = complexarr( sim[1], sim[2] )
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 * float( FFT( imFT * conj( psf_FT ), 1 ) ) $
else conv = npix * float( FFT( imFT * psf_FT, 1 ) )
sc = sc + (sim MOD 2) ;shift correction for odd size images.
return, shift( conv, sc[1], sc[2] )
end