function minmax,array,subs,NAN=nan, DIMEN=dimen
;+
; NAME:
;      MINMAX
; PURPOSE:
;      Return a 2 element array giving the minimum and maximum of an
;      array
; EXPLANATION:
;      Using MINMAX() is faster than doing a separate MAX and MIN.
;
;      The procedure MAXMIN in
;      http://www.dfanning.com/programs/maxmin.pro
;      has a similar purpose but uses a procedure call rather than a
;      function.
; CALLING SEQUENCE:
;      value = minmax( array, [subs, /NAN, DIMEN= ] )
; INPUTS:
;      array - an IDL numeric scalar, vector or array.
;
; OUTPUTS:
;      value = a two element vector (if DIMEN is not supplied)
;            value[0] = minimum value of array
;            value[1] = maximum value of array
;
;            If the DIMEN keyword is supplied then value will be a 2 x
;            N element
;            array where N is the number of elements in the specified
;            dimension
;              
; OPTIONAL OUTPUT PARAMETER:
;       subs - two-dimensional vector; the first element gives the
;              subscript
;              of the minimum value, the second element gives the
;              subscript
;              of the maximum value.     
;
; OPTIONAL INPUT KEYWORD:
;      /NAN   - Set this keyword to cause the routine to check for
;               occurrences
;            of the IEEE floating-point value NaN in the input data.
;            Elements 
;            with the value NaN are treated as missing data.
;
;      DIMEN -  integer (either 1 or 2) specifying which dimension of
;              a 2-d 
;            array to  take the minimum and maximum.   Note that
;            (unlike the
;            DIMENSION keyword to the MIN() function) DIMEN is only
;            valid
;            for a 2-d array, larger dimensions are  not supported.  
; EXAMPLE:
;     (1)  Print the minimum and maximum of an image array, im
; 
;            IDL> print, minmax( im )
;
;     (2) Given a 2-dimension array of (echelle) wavelengths w, print
;     the
;         minimum and maximum of each order
;
;         print,minmax(w,dimen=1)
;
; PROCEDURE:
;      The MIN function is used with the MAX keyword
;
; REVISION HISTORY:
;      Written W. Landsman                January, 1990
;      Added NaN keyword.      M. Buie       June 1998
;      Added DIMEN keyword    W. Landsman  January 2002
;      Added SUBSCRIPT_MIN and SUBSCRIPT_MAX  BT Jan 2005
;      Added optional subs output parameter  W. Landsman July 2009
;-
 On_error,2
 compile_opt idl2
 if N_elements(DIMEN) GT 0 then begin
      amin = min(array, MAX = amax, NAN = nan, DIMEN = dimen,cmin,sub=cmax) 
      if arg_present(subs) then subs = transpose([[cmin], [cmax]])
      return, transpose([[amin],[amax] ])
  endif else  begin 
     amin = min( array, MAX = amax, NAN=nan, cmin, sub=cmax)
      if arg_present(subs) then subs = [cmin, cmax]
     return, [ amin, amax ]
 endelse
end



pro ms_plotimage_varbins,image2d,left_xarray,right_xarray, $
                         left_yarray, right_yarray, cbrange, logcolor = logcolor, convolve = convolve, $
                         OOB_HIGH = OOB_HIGH, OOB_LOW = OOB_LOW, _extra = _extra, had_high = had_high, $
                         had_low = had_low, ncolors = ncolors, bottom = bottom
;;;;; Purpose: Plots the 2-d histogram of scatter data points in data co-ordinates
;;;;; Inputs:  The 2-d histogram, the arrays of x & y bounds, the
;;;;; colorbar range (i.e., the density in the 2-d histogram to plot)
;;;;; and the colorbar title. cbrange should always be sent in
;;;;; "normal" units (ie without a log10 even if logcolor is set)
;;;;; Keywords: 
;;;;; logcolor -- Takes a log10 of the 2-d histogram before
;;;;;             plotting. cbrange is changed to
;;;;;             log10(cbrange). Plotting the color bar should be
;;;;;             fairly straightforward in the calling routine.
;;;;;
;;;;; (For the pedantic, technically this is not a logcolor. This is a
;;;;; linear map applied on the log of the quantity. Details..)
;;;;;
;;;;; The current color table is used. If this is used on a postscript
;;;;; device with a version of IDL that does not support decomposed=1,
;;;;; the code will break. Obviously. You are on your own there,
;;;;; upgrade your IDL license, or use the "sm" variable to set your
;;;;; plotting color directly. However, the plot will likely look
;;;;; horrible (depending on the color table loaded - since
;;;;; !p.color/background will be set to the color table; good luck).
;;;;;
;;;;;
;;;;; This is how I envision that data for image2d, (left/right)_xarray and (left/right)_yarray
;;;;; was created:
;;;;; 
;;;;; image2d = hist_nd(transpose([[x],[y]]),[xbinsize,ybinsize],min=[xmin,ymin],max=[xmax,ymax])
;;;;; nx = (size(image2d))[1]
;;;;; ny = (size(image2d))[2]
;;;;; left_xarray  = findgen(nx)*xbinsize + xmin
;;;;; right_xarray = left_xarry + xbinsize
;;;;; left_yarray  = findgen(ny)*ybinsize + ymin
;;;;; right_yarray = left_yarray + ybinsize
;;;;; cbrange = minmax(image2d) 
;;;;; logcolor = 0
;;;;; Ncolors = 100
;;;;; OOB_HI  = 'red'
;;;;; OOB_LOW = 'white'
;;;;  cgloadct, 4, Ncolors=Ncolors
;;;;; ms_plotimage_varbins, image2d, left_xarry, right_xarry, $
;;;;; left_yarray, right_yarray, cbrange, logcolor = [0/1], convole = [0/1],$
;;;;; position=position, OOB_HI=OOB_HI,OOB_LOW=OOB_LOW,had_hi = had_hi, $
;;;;; had_low = had_low, Ncolors=Ncolors
;;;;; if had_hi eq 0 then delvarx, OOB_HI
;;;;; if had_low eq 0 then delvarx, OOB_LOW
;;;;; cgcolorbar, range = logcolor eq 1 ? alog10(cbrange):cbrange,$
;;;;; /fit, OOB_HI=OOB_HI, OOB_LOW=OOB_LOW, Ncolors=Ncolors
;;;;; 
;;;;; 
;;;;; Rev 1.00 : MS 12/11/2014.

compile_opt idl2,strictarrsubs
on_error,0 ;;; stop inside this routine on error


if n_elements(image2d) eq 0 or n_elements(left_xarray) eq 0 or n_elements(right_xarray) eq 0 $
   or n_elements(left_yarray) eq 0 or n_elements(right_yarray) eq 0 or n_elements(cbrange) eq 0 then begin
   print, 'usage: ms_plotimage image2d [2-d image data] left_xarray right_xarray [x min/max values] left_yarray right_yarray [y min/max values] cbrange [range of data in image2d to plot]'
   print, 'returning..'
    return
endif

sz = size(image2d)
if size(image2d, /n_dim) ne 2 then $
  message, 'image2d must be a 2D array'

if size(left_xarray,/n_dim) ne 1 or size(right_xarray, /n_dim) ne 1 then $
  message, 'xarrays must be an 1D'

if n_elements(left_xarray) ne n_elements(right_xarray) then $
   message, 'xarrays must have the same number of elements'

if size(left_yarray,/n_dim) ne 1 or size(right_yarray, /n_dim) ne 1 then $
  message, 'yarrays must be an 1D'

if n_elements(left_yarray) ne n_elements(right_yarray) then $
   message, 'yarrays must have the same number of elements'
   

if size(cbrange,/n_dim) ne 1 then $
  message,'cbrange must be an 1D array'

if n_elements(cbrange) ne 2 then $
  message,'cbrange must contain only two elements: min and max of the data to plot'

device,get_decomposed=olddecomposed
device,decomposed=1


if n_elements(bottom) eq 0 then $
   BOTTOM = 0

if n_elements(ncolors) eq 0 then $
   NCOLORS = 256

TOP = BOTTOM+NCOLORS-1 < 255


image = double(image2d)
fixed_left_x  = double(left_xarray)
fixed_right_x = double(right_xarray)

fixed_left_y = double(left_yarray)
fixed_right_y = double(right_yarray)

cb_range = double(cbrange)

nx = n_elements(fixed_left_x)
ny = n_elements(fixed_left_y)
if nx ne sz[1] or ny ne sz[2] then begin
   message, 'Number of elements in xarray/yarray does not match that in image2d'
endif

if n_elements(xrange) eq 0 then begin
   xrange = minmax([fixed_left_x, fixed_right_x], /nan)
endif

if n_elements(yrange) eq 0 then begin
   yrange = minmax([fixed_left_y, fixed_right_y], /nan)
endif
   
if keyword_set(logcolor) then log_color = 1 else log_color = 0
if min(cb_range) le 0.0  and log_color eq 1 then begin
    print,'WARNING: log color is enabled but the min. value of the image is less than 0.0'
    stop
endif
if keyword_set(convolve) then begin
    kernel =double( [$
                      [ 1,  8, 15,  8, 1], $
                      [ 8, 63,127, 63, 8], $
                      [15,127,255,127,15], $
                      [ 8, 63,127, 63, 8], $
                      [ 1,  8, 15,  8, 1]] ) ;;; gaussian kernel 

    convol_image = convol(image,kernel,invalid=!VALUES.F_NAN, missing=!VALUES.F_NAN,/edge_truncate,/normalize,/center)
    image = convol_image
endif

actual_cb_range = cb_range
plot,[0],/nodata,xrange=xrange,yrange=yrange,_extra=_extra

if log_color  eq 1 then begin
    xx = where(finite(image) eq 1 and image gt 0.0,cnt)
    yy  = where(finite(image) eq 1 and image le 0.0,cnt1)
    tiny = 1d-5
    if cnt1 gt 0 then begin
        minimg = min(image[xx],/NAN)
        image[yy] = tiny*minimg
        if minimg lt cb_range[0] then minimg = cb_range[0]
    endif else begin
        minimg = cb_range[0]
    endelse
    cb_range = alog10(cb_range)
    log_image = alog10(image)
endif else begin
    minimg = min(image,/NAN)
    if minimg lt cb_range[0] then minimg = cb_range[0]
endelse

tvlct,red,green,blue,/get
had_high=0
had_low=0
cg_oob_hi = cgcolor(OOB_HI)
cg_oob_lo = cgcolor(OOB_LOW)

for iy = 0l, ny-1 do begin
   if fixed_left_y[iy] ge yrange[1] then break
   if fixed_right_y[iy] lt yrange[0] then continue
   ind = where(image[*, iy] ge actual_cb_range[0] and image[*, iy] le actual_cb_range[1] and finite(image[*, iy]) eq 1, cnt)
   if cnt gt 0 then begin
      if log_color eq 0 then $
         sm = reform(image[ind, iy]) $
      else $
         sm = reform(log_image[ind, iy])
      
      if min(sm) eq max(sm) then begin
         sm = (sm - cb_range[0]) /(cb_range[1] - cb_range[0])*Ncolors + BOTTOM
      endif else begin
         min_sm = min(sm, max = max_sm)
         scaled_range = (max_sm-min_sm)/(cb_range[1]-cb_range[0])*Ncolors
         scaled_min = ((min_sm-cb_range[0])/(cb_range[1]-cb_range[0])*Ncolors + BOTTOM) > BOTTOM
         scaled_max = scaled_min + scaled_range
         if scaled_min eq scaled_max then stop
         sm = scale_vector(sm, scaled_min, scaled_max)
      endelse
      sm = (fix(sm) < TOP ) > BOTTOM
      color = long(red[sm] + 256L*green[sm] + 256L^2*blue[sm])
      for ix = 0l, n_elements(ind)-1 do begin
         if fixed_left_x[ind[ix]] ge xrange[1] then break
         if fixed_right_x[ind[ix]] lt xrange[0] then continue
         xpoly = [fixed_left_x[ind[ix]], fixed_right_x[ind[ix]], fixed_right_x[ind[ix]], $
                  fixed_left_x[ind[ix]], fixed_left_x[ind[ix]]]
         ypoly = [fixed_left_y[iy], fixed_left_y[iy],                             fixed_right_y[iy], $
                  fixed_right_y[iy], fixed_left_y[iy]]
         polyfill, xpoly, ypoly, color = color[ix], /data, noclip = 0
      endfor
   endif

    ;;; Now plot the out of bounds
    if n_elements(OOB_LOW) gt 0 then begin
        ind = where(image[*,iy] lt actual_cb_range[0] and finite(image[*,iy]) eq 1,cnt)
        if cnt gt 0 then begin
            for ix=0l,n_elements(ind)-1 do begin
               if fixed_left_x[ind[ix]] ge xrange[1] then break
               if fixed_right_x[ind[ix]] lt xrange[0] then continue
               xpoly = [fixed_left_x[ind[ix]], fixed_right_x[ind[ix]], fixed_right_x[ind[ix]], $
                        fixed_left_x[ind[ix]], fixed_left_x[ind[ix]]]
               ypoly = [fixed_left_y[iy], fixed_left_y[iy],                             fixed_right_y[iy], $
                        fixed_right_y[iy], fixed_left_y[iy]]
               polyfill, xpoly, ypoly, color = cg_oob_lo, /data, noclip = 0
            endfor
            had_low=1
        endif
    endif

    if n_elements(OOB_HIGH) gt 0 then begin
        ind = where(image[*,iy] gt actual_cb_range[1] and finite(image[*,iy]) eq 1,cnt)
        if cnt gt 0 then begin
            for ix=0l,n_elements(ind)-1 do begin
               if fixed_left_x[ind[ix]] ge xrange[1] then break
               if fixed_right_x[ind[ix]] lt xrange[0] then continue
               xpoly = [fixed_left_x[ind[ix]], fixed_right_x[ind[ix]], fixed_right_x[ind[ix]], $
                        fixed_left_x[ind[ix]], fixed_left_x[ind[ix]]]
               ypoly = [fixed_left_y[iy], fixed_left_y[iy],                             fixed_right_y[iy], $
                        fixed_right_y[iy], fixed_left_y[iy]]
                polyfill,xpoly,ypoly,color=cg_oob_hi,/data,noclip=0
            endfor
            had_high=1
        endif
    endif
endfor

plot,[0],/nodata,/noerase,xrange=xrange,yrange=yrange,_extra=_extra
device,decomposed=olddecomposed

end