function orbits_mfilter, image, number, threshold, kill = kill, $
                         fill_spatial = fill_spatial, no_point = no_point

;+
; ORBITS_MFILTER
;	Make an image that is the result of median filtering in the time-domain
;
; Usage:
;	im=orbits_mfilter(image, number)
;
; Return value:
;	A floating-point array with the image.
;
; Arguments:
;	image	obj	The smei_image object with the central orbit.
;	number	int	The number of orbits before and after to
;			compare
;	threshold float	The threshold at which to replace the value.
;
; Keyword:
;	/kill		If set, then replace the value with NaN
;			instead of the median.
;	fill_spatial int If this is set to a non-zero value then
;			replace flagged pixels with the nxn spatial
;			median rather than a temportal median.
;	/no_point	Do not stop when frames have different pointing.
;
; Example:
;	In a sop file to do a 5-orbit filter with a threshold of 3ADU use:
;	(orbits_mfilter ss{-1} #2 #3.0)
;
; History:
;	Original (after orbits average): 13/11/09; SJT
;-

level = image -> get_pipe_level()
compact = image -> get_pflags(/compact)

sz = image -> get_size()
img = fltarr([sz, 2*number+1])

; The current image

im0 = image -> get_image()
im1 = im0

if keyword_set(fill_spatial) then sim = median(im0, fill_spatial)

; Don't include pixels with no content, the level 1 compact is a "best
; guess"

if (~compact) then begin
    imw = image -> get_raw_image(/weight)
    locs = where(imw eq 0)
endif else if level ge 2 then begin
    imf = image -> get_flag_data()
    locs = where(imf eq 'ffff'xu)
endif else locs = where(im0 eq 0)
im1[locs] = !values.f_nan 

img[*, *, 0] = im1   ; Note that the order we store them doesn't matter

ih = image -> get_header(/fits)
solon = sxpar(ih, "SUN_ELOP", count = ns)
if ns eq 0 then begin
    solon = round(sxpar(ih, "SUN_ELON"))
    rflag = 0b
endif else rflag = 1b

imn = image -> get_next()
imp = image -> get_prev()

for j = 0, number-1 do begin

; Images after the current are stored in planes 1,3,...
    if obj_valid(imn) then begin
        ih = imn -> get_header(/fits)
        if rflag then sn = sxpar(ih, "SUN_ELOP") $
        else sn = round(sxpar(ih, "SUN_ELON"))
        if sn eq solon or keyword_set(no_point) then begin
            im = imn -> get_image()

            if (~compact) then begin
                imw = image -> get_raw_image(/weight)
                locs = where(imw eq 0)
            endif else if level ge 2 then begin
                imf = image -> get_flag_data()
                locs = where(imf eq 'ffff'xu)
            endif else locs = where(im eq 0)
            im[locs] = !values.f_nan 
            img[*, *, 2*j+1] = im
        endif else img[*, *, 2*j+1] = !values.f_nan
        imn = imn -> get_next()
    endif else img[*, *, 2*j+1] = !values.f_nan

; Images before the current are stored in planes 2,4,...
    if obj_valid(imp) then begin
        ih = imp -> get_header(/fits)
        if rflag then sn = sxpar(ih, "SUN_ELOP") $
        else sn = round(sxpar(ih, "SUN_ELON"))
        if sn eq solon  or keyword_set(no_point) then begin
            im = imp -> get_image()

            if (~compact) then begin
                imw = image -> get_raw_image(/weight)
                locs = where(imw eq 0)
            endif else if level ge 2 then begin
                imf = image -> get_flag_data()
                locs = where(imf eq 'ffff'xu)
            endif else locs = where(im eq 0)
            im[locs] = !values.f_nan 
            img[*, *, 2*j+2] = im
        endif else img[*, *, 2*j+2] = !values.f_nan
        imn = imp -> get_prev()
    endif else img[*, *, 2*j+2] = !values.f_nan
endfor

mim = median(img, dimension = 3)
locs = where(abs(im0-mim) gt threshold and finite(im1), nfl)

if nfl ne 0 then begin
    if keyword_set(kill) then im0[locs] = !values.f_nan $
    else if keyword_set(fill_spatial) then im0[locs] = sim[locs] $
    else im0[locs] = mim[locs]
endif

return, im0

end