function fft_clean, image, scale, order = order

;+
; FFT_CLEAN
;	Apply an FFT-based low-pass filter to an image.
;
; Usage:
;	cimg=fft_filter(image[, scale])
;
; Return Value:
; 	cimg	The filtered image (float unless the input is double)
;
; Arguments:
;	image	numeric	The image to be filtered.
;	scale	float	The scale of the filter (larger is a stronger filter).
;
; Keyword:
;	order	float	The power to raise the filter.
;
; Notes:
;	The filter is created as exp(-(dist(d0,d1)/scale)^order),
;	where the d's are the dimensions of the mapping image (the
;	next power of 2 greater than twice the dimensions of the input
;	image).
;
; History:
;	Original: 4/8/10; SJT
;-

if ~keyword_set(order) then order = 2.0

sz = size(image)
if sz[0] ne 2 then begin
    print, "WARNING input image is not a 2-D array"
    return, image
endif

d = 2l^ceil(alog(sz[1:2])/alog(2.)+1)
if sz[3] eq 5 then ei = dcomplexarr(d) $
else ei = complexarr(d)

llc = d/4
ei[llc[0], llc[1]] = image


tei = fft(ei)
if n_params() eq 1 then scale = float(min(llc))

md = max(d)
filt = exp(-(dist(md)*2.d0*scale/md)^order)
filt /= max(filt)
filt = rebin(filt, d)

fei = fft(tei*filt, 1)
if sz[3] eq 5 then $
  rv = double(fei[llc[0]:llc[0]+sz[1]-1, llc[1]:llc[1]+sz[2]-1]) $
else rv = float(fei[llc[0]:llc[0]+sz[1]-1, llc[1]:llc[1]+sz[2]-1])

return, rv

end