;+ ;NAME: ; CONVOLF ;PURPOSE: ; Similar CONVOL with the /edge_truncate option set, but it ; does convolution in Fourier space. Much more efficient! ;CALLING SEQUENCE: ; result = convolf(image, kernel [, reps=reps]) ;INPUT PARAMETERS: ; IMAGE: The 2D image to be convolved with KERNEL. ; KERNEL: The KERNEL array to be convolved with IMAGE. If ; the dimensions of kernel are odd numbers, then it ; will be centered properly. If not, well, hey, you were ; warned. ;OPTIONAL INPUT KEYWORDS: ; REPS: Number of times to convolve image with kernel. Default=1. ; COMPLEX: If set, then allow complex results. By default, the ; result is returned as a float. This keyword may also be used ; to allow double precision, though the result will be complex. ;MODIFICATION HISTORY: ; 2006-Jun-08 C. Kankelborg. ;- function CONVOLF, image, kernel, reps=reps, complex=complex if n_elements(reps) ne 1 then reps=1 ;default isize = size(image) Nx = isize[1] Ny = isize[2] ksize = size(kernel) knx = ksize[1] kny = ksize[2] ;Make a periodic image by 'unfolding' image twice ;(new dimensions are 2*Nx by 2*Ny) image2 = [ image, rotate(image, 5) ] image3 = [ [image2],[rotate(image2, 7)] ] ;FFT the image imageF = fft(image3) ;Embed the kernel in Nx,Ny sized array kernel2 = fltarr(2*Nx,2*Ny) kernel2[0:knx-1, 0:kny-1] = kernel ;Shift the kernel center to 0,0 and FFT kernel3 = shift(kernel2,-knx/2,-kny/2) ;FFT the kernel kernelF = 4.0*Nx*Ny * fft(kernel3) ;Note careful normalization. ;Multiply in Fourier space and transform back. if (not keyword_set(complex)) then $ image = float( fft(kernelF^reps * imageF, /inverse) ) $ else image = fft(kernelF^reps * imageF, /inverse) ;Return lower left quadrant. return, image[0:Nx-1, 0:Ny-1] end