;+
;
; NAME: reduce.pro
;
; USAGE: sparse_stuff = reduce(data, xfactor, sx=sx, sy=sy)
;
; PURPOSE: Truncates and/or reduces resolution of 2-d "data" array, 
;      for easier visualization by VELOVECT, etc.
;
;      Accomplished by finding number of times "[xy]factor" divides 
;      into "data" and calling REBIN, which uses blinear interp- 
;      olation (better than CONGRID, which uses sampling.)
;
; DEFAULT: Reduce resolution by factor of 2, without truncation.
;
; REQUIRED ARGUMENTS:
;     data := (nx,ny) array of data, of any type (but complex?)
;
; KEYWORD PARAMETERS:
;     xfactor := if set, keeps every "xfactor" point in 1st dim.
;     yfactor := if set, keeps every "yfactor" point in 2nd dim.
;                DEFAULT: yfactor = xfactor  
;     sub := truncates dimensions of "data" to central "xsub" fraction
;                DEFAULT: sub = 1. 
;     xsub := truncates 1st dim of "data" to central "xsub" fraction
;                DEFAULT: xsub = sub 
;     ysub := truncates 2nd dim of "data" to central "ysub" fraction
;                DEFAULT: ysub = sub 
;     tx := floating array of indices for 1st dim. of truncated data
;     ty := floating array of indices for 2nd dim. of truncated data
;     sx := floating array of indices for 1st dim. of sparsified data
;     sy := floating array of indices for 2nd dim. of sparsified data
;
; HISTORY: scraped together from other routines by BTW on 29 April '03
; 
;

function reduce, data, xfactor, yfactor, tdata=tdata, sub=sub, $
            xsub=xsub, ysub=ysub, tx=tx, ty=ty, sx=sx, sy=sy

nparms=N_Params()
if (nparms eq 0) then begin
   Message, '', /info   
   Message, 'NAME: reduce.pro', /info
   Message, '', /info   
   Message, 'USAGE: result = reduce(data, xfactor, sx=sx, sy=sy)', /info
   Message, '', /info   
   Message, 'PURPOSE: Truncates &/or reduces resolution of 2-d array', /info 
   Message, '    "data", by "xfactor", for display w/VELOVECT, etc.', /info
   Message, '', /info   
   Message, '    If "[xy]sub" is set, truncates "data" outside input', /info   
   Message, '    fraction(s) of original array.', /info   
   Message, '', /info   
   Message, '    Then finds number of times "[xy]factor" divides', /info
   Message, '    into "data", truncates remainder, and calls REBIN. ', /info   
   Message, '', /info   
   Message, '    Default: no truncation, halves resolution.', /info   
   Message, '', /info   
   Message, 'lo_res = reduce(hi_res) halves resolution in x & y', /info
   Message, 'lo_res = reduce(hi_res,4) drops res. by 4 in x & y',/info
   Message, 'lo_res = reduce(hi_res,4,2) drops res.by 4 in x, 2 in y',/info
   Message, 'lo_res = reduce(hi_res,4,2,sx=sx,sy=sy) drops res. by 4 & 2',/info
   Message, '         (in x & y), and returns coord arrays for lo_res',/info 
   Message, '', /info   
   Message, 'WRITTEN: B Welsch UCB/SSL 29 April, 2002', /info
   return, -1
endif

if (not(keyword_set(sub))) then sub = 1.
if (not(keyword_set(xsub))) then xsub = sub
if (not(keyword_set(ysub))) then ysub = sub

if (not(keyword_set(xfactor))) then xfactor = 2.
if (not(keyword_set(yfactor))) then yfactor = xfactor

; get relevant sizes
dims = size(data, /dim)
nx = dims(0) 
ny = dims(1)
x = findgen(nx)
y = findgen(ny)

; if called for, truncate data to "sub" fraction around image center
if (xsub eq 1.) then xmin = 0 else $
  xmin = max([    0, round((.5-.5*xsub) *nx)])    
if (xsub eq 1.) then xmax = nx-1 else $
  xmax = min([ nx-1, round((.5+.5*xsub) *nx ) - 1])       
if (ysub eq 1.) then ymin = 0 else $
        ymin = max([    0, round((.5-.5*ysub) *ny )])    
if (ysub eq 1.) then ymax = ny-1 else $
        ymax = min([ ny-1, round((.5+.5*ysub) *ny ) - 1])     


; width of truncated array, and new x & y arrays
dx =  xmax-xmin+1
dy =  ymax-ymin+1

; ordinate & abscissa for plotting truncated results
tx =  findgen(dx) + xmin
ty =  findgen(dy) + ymin


; make arrays sparse
case 1 of
    ((xfactor ne 1) or (yfactor ne 1)): begin
        ; if (__factor) doesn't divide evenly, compute remainder
        xremain = float((dx)mod(xfactor))
        yremain = float((dy)mod(yfactor))

        ; displace start from left/bottom edges by half remainder
        dxmin = round(xremain/2.) ; errs to right by 1 if split pix.
        dymin = round(yremain/2.) ; errs to   top by 1 if split pix.

        ; how many elements in each dim'n of sparse arrays
        n_sx = (dx - xremain)/xfactor
        n_sy = (dy - yremain)/yfactor
        
        ; ordinate & abscissa for plotting sparsified results
        sxmin = xmin + dxmin
        symin = ymin + dymin

        sxmax = sxmin + n_sx*xfactor - 1
        symax = symin + n_sy*yfactor - 1

        sx = sxmin + xfactor*findgen(n_sx)
        sy = symin + yfactor*findgen(n_sy)

        tdata =  data( sxmin:sxmax, symin:symax)
        sdata = rebin( tdata, n_sx, n_sy, /sample) 

    END
    else: begin  ; no sparsification!  
        sx = tx
        sy = ty
        tdata =  data( xmin:xmax, ymin:ymax)
        sdata = tdata ; sparse and truncated data are identical!
    end
endcase 

return, sdata

end