function median_filtered_mean,vec
;
;   computes the median filtered mean value of the data in the
;   vec vector
;
;   returns the mean and standard deviation
;
;   DPC     20040101  legacy code
;           20051012  trapped stdev to no smaller than 0.001 to halt sigma collapse
;           20051023  tuned to improve speed
;
this_vec = vec
nvec = size(vec)
nvec = nvec[1]
;
value = fltarr(2)
;
case nvec of
    0: begin
         value[0] = -1.0e6
         value[1] = -1.0e6
       end
    1: begin
           value[0] = vec[0]
           value[1] = vec[0]
       end
    else: begin
         ;
         i = 0L
         ind = INDGEN(nvec)   ; indices of all elements
         nv1 = nvec-1L
         old_count = 0
         for iloop = 0, 4 do begin
             scale = 2.2 + 0.2 * iloop
             the_vector = this_vec[ind]
             med = d_MEDIAN(the_vector)     ; compute the median for this pixel for all allowed elements
             mx = d_MOMENT(the_vector)
             stddev = 0.001 > SQRT(mx[1])         ; standard dev, with trap to prevent zero
             ;
             ;   keep or reject each original data value
              ;
             ind = INDGEN(nvec)
             for i = 0L, nv1 do begin       ; for each element
                 if(abs(this_vec[i]-med) gt scale*stddev) then ind[i] = -1
             endfor
             ind = WHERE(ind ne -1,count)         ; turn on only values within scale*sigma of median
             ;
             ;   test for count -- if no change with loop, quit
             ;
             if(count eq old_count) then goto, done else old_count = count
         endfor
         ;
         done:
         mx = d_MOMENT(this_vec[ind])
         value[0] = mx[0]     ; the final mean value
         value[1] = sqrt(mx[1])
         ;
       end
endcase
;
return, value
end