;+
; NAME:
;   PLOTIMAGE
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;
; PURPOSE:
;   Displays an image via a "PLOT"-like interface.
;
; CALLING SEQUENCE:
;   PLOTIMAGE, img, [xrange=xrange,] [yrange=yrange,] ...
;
; DESCRIPTION: 
;
;   PLOTIMAGE displays an image (or slice of an image) on the current
;   graphics device.  The syntax is very similar to the PLOT command,
;   in the sense that an XRANGE and YRANGE for the plot can be
;   specified.  
;
;   More importantly, coordinate ranges can be specified for the
;   *image* itself, and then PLOTIMAGE will intelligently crop the
;   image so that only the part that falls within the XRANGE and
;   YRANGE is displayed.  Images often have physical units (other than
;   pixels) assigned to their X and Y dimensions.  PLOTIMAGE allows
;   you to slice the image appropriately and display appropriate
;   coordinate axes.
;
;   Image coordinates are specified via the IMGXRANGE and IMGYRANGE
;   keywords.  IMGXRANGE specifies the x-values for the first and last
;   pixels in each row.  IMGYRANGE gives the y-values for the first
;   and last pixels in each column.  Default image coordinates are in
;   units of pixels.
;
;   Plot coordinates are specified in the usual XRANGE and YRANGE
;   keywords.  The [XY]RANGE may specify a range smaller than the
;   image size, so that the image is cropped; or a range larger than
;   the image size, in which case the appropriate portion of the image
;   is displayed.
;
; INPUTS:
;
;   IMG - A byte array to be displayed.  An image declared as
;         ARRAY(M,N) will be M pixels in the x-direction and N pixels
;         in the y-direction.  The image is resampled via
;         interpolation to fill the desired display region.
; 
; OPTIONAL INPUTS:
;   NONE
;
; INPUT KEYWORD PARAMETERS:
;
;   XRANGE, YRANGE - Each is a two component vector that specifies the
;                    X and Y plot ranges, respectively.  These values
;                    are not required to coincide with IMG[XY]RANGE.
;                    Default: the image ranges
;
;   IMGXRANGE, IMGYRANGE - Each is a two component vector that
;                          describes the X and Y position of the first
;                          and last pixels.
;                          Default: the size of the image in pixels
;
;   POSITION - Position of the inner plot window in the standard
;              graphics keyword format.  Overrides PANEL and SUBPANEL.
;
;   PANEL, SUBPANEL - An alternate way to more precisely specify the
;                     plot and annotation positions.  See SUBCELL.
;
;   NOERASE - If set, the display is not erased before graphics
;             operations.
;
;   NODATA - If set, the image is not actually displayed, but
;            coordinate axes may be drawn.
;
;   NOAXES - An attempt is made to render the image without coordinate
;            axes.  This is possible if POSITION, PANEL, or SUBPANEL
;            are  given.
;
;   PLOTIMAGE will pass other keywords directly to the PLOT command
;   used for generating the plot axes.  XSTYLE=1 and YSTYLE=1 are
;   enforced.
;
; OUTPUTS:
;   NONE
;
; PROCEDURE:
;
; EXAMPLE:
;
;   This example constructs an image whose values are found by
;       z(x,y) = cos(x) * sin(y)
;   and x and y are in the range [-2,2] and [4,8], respectively.
;   The image is then plotted, with appropriate axes.
;   
;   x = findgen(20)/5. - 2.
;   y = findgen(20)/5. + 4.
;   zz = cos(x) # sin(y)
;   imgxrange = [min(x), max(x)]
;   imgyrange = [min(y), max(y)]
;   plotimage, bytscl(zz), imgxrange=imgxrange, imgyrange=imgyrange
;
;   This second example plots the same image, but with a plot range
;   much larger than the image's.
;
;   xr=[-10.,10]
;   yr=[-10.,10]
;   plotimage, bytscl(zz), imgxrange=imgxrange, imgyrange=imgyrange, $
;      xrange=xr, yrange=yr
;
; SEE ALSO:
;
;   OPLOTIMAGE, BYTSCL
;
; EXTERNAL SUBROUTINES:
;
;   SUBCELL, DEFSUBCELL, TVIMAGE
;
; MODIFICATION HISTORY:
;   Written, CM, 1997
;
;-

forward_function defsubcell, subcell

function defsubcell, default

  if n_elements(default) EQ 0 then default = [-1.,-1,-1,-1]
  mysubcell = default
  defaultsubpos = [ 0.08, 0.08, 0.95, 0.95 ]

  iwh = where(mysubcell LT 0, ict)
  if ict GT 0 then $
    mysubcell(iwh) = defaultsubpos(iwh)

  return, mysubcell
end
function subcell, subpos, position, margin=margin

  ;; Default value for subposition
  if n_elements(subpos) EQ 0 then mysubpos = [-1.,-1,-1,-1] $
  else mysubpos = subpos

  ;; Default value for position - full screen
  if n_elements(position) EQ 0 then position = [0.,0.,1.,1.]

  ;; Get margins if necessary
  if keyword_set(margin) EQ 1 OR n_elements(subpos) EQ 0 then $
    mysubpos = defsubcell(mysubpos)

  ;; Compute new window position
  x0 = position(0)
  y0 = position(1)
  dx = position(2)-position(0)
  dy = position(3)-position(1)

  newsubpos = reform(mysubpos * 0, 4)
  newsubpos([0,2]) = x0 + dx * mysubpos([0,2])
  newsubpos([1,3]) = y0 + dy * mysubpos([1,3])

  return, newsubpos
end
  
pro plotimage, img, xrange=xrange, yrange=yrange, $
               imgxrange=imgxrange, imgyrange=imgyrange, $
               position=position, panel=panel, subpanel=subpanel, $
               noerase=noerase, nodata=nodata, noaxes=noaxes, $
               pixtolerance=pixtolerance, $
               _EXTRA=extra

  ;; Usage message
  if n_params() EQ 0 then begin
      message, 'PLOTIMAGE, image, xrange=, yrange=, imgxrange=, imgyrange=,..'$
        , /info
      return
  endif

  ;; Must have a byte-scaled image already
  imgsize  = size(img)
  if imgsize(imgsize(0)+1) NE 1 then begin
      message, 'ERROR: image must be of type BYTE', /cont
      return
  endif
  
  ;; Parameter checking
  if keyword_set(nodata) then mynodata = 1 else mynodata = 0
  if n_elements(pixtolerance) EQ 0 then pixtolerance = 1.e-2
  imgpanel = [0., 0., 1., 1.]

  ;; Image size
  imgdim   = imgsize(0)
  imgxsize = imgsize(1)
  if imgdim NE 2 then $
    message, 'ERROR: image must have 2 dimensions'
  imgysize = imgsize(2)

  ;; By default, we have no info about the image, and display the
  ;; whole thing
  if n_elements(imgxrange) EQ 0 then $
    imgxrange = [ 0L, imgxsize ]
  if n_elements(xrange) EQ 0 then $
    xrange = imgxrange

  if xrange(1) LT xrange(0) then $
    message, 'ERROR: Do not reverse XRANGE in PLOTIMAGE'

  srcxpix  = [ 0L, imgxsize-1 ]
  ;; Size of one x pix
  dx = (imgxrange(1) - imgxrange(0)) / imgxsize 

  if xrange(0) GE imgxrange(1) OR xrange(1) LE imgxrange(0) then begin
      message, 'WARNING: No image data in specified XRANGE.', /info
      mynodata = 1
      goto, PLOTDATA
  endif

  ;; Case where xrange cuts off image at left
  if xrange(0) GT imgxrange(0) then begin
      offset = double(xrange(0)-imgxrange(0))/dx
      srcxpix(0) = long( offset )
      froffset = offset - long(offset)
      if abs(froffset) GT pixtolerance then $
        xrange(0) = xrange(0) - froffset * dx
  endif

  ;; Case where xrange cuts off image at right
  if xrange(1) LT imgxrange(1) then begin
      offset = double(imgxrange(1)-xrange(1))/dx
      srcxpix(1) = imgxsize - long(offset) - 1
      froffset = offset - long(offset)
      if abs(froffset) GT pixtolerance then begin
          xrange(1) = xrange(1) + froffset * dx
          srcxpix(1) = srcxpix(1) + 1
      endif
  endif 

  if xrange(0) LT imgxrange(0) then $
    imgpanel(0) = (imgxrange(0) - xrange(0))/(xrange(1)-xrange(0))
  if xrange(1) GT imgxrange(1) then $
    imgpanel(2) = (imgxrange(1) - xrange(0))/(xrange(1)-xrange(0))

  ;; By default, we have no info about the image, and display the
  ;; whole thing
  if n_elements(imgyrange) EQ 0 then $
    imgyrange = [ 0L, imgysize ]
  if n_elements(yrange) EQ 0 then $
    yrange = imgyrange

  if yrange(1) LT yrange(0) then $
    message, 'ERROR: Do not reverse YRANGE in PLOTIMAGE'

  srcypix  = [ 0L, imgysize-1 ]
  ;; Size of one y pix
  dy = (imgyrange(1) - imgyrange(0)) / imgysize 

  if yrange(0) GE imgyrange(1) OR yrange(1) LE imgyrange(0) then begin
      message, 'WARNING: No image data in specified YRANGE.', /info
      mynodata = 1
      goto, PLOTDATA
  endif

  ;; Case where yrange cuts off image at bottom
  if yrange(0) GT imgyrange(0) then begin
      offset = double(yrange(0)-imgyrange(0))/dy
      srcypix(0) = long( offset )
      froffset = offset - long(offset)
      if abs(froffset) GT pixtolerance then $
        yrange(0) = yrange(0) - froffset * dy
  endif


  ;; Case where yrange cuts off image at top
  if yrange(1) LT imgyrange(1) then begin
      offset = double(imgyrange(1)-yrange(1))/dy
      srcypix(1) = imgysize - long(offset) - 1
      froffset = offset - long(offset)
      if abs(froffset) GT pixtolerance then begin
          yrange(1) = yrange(1) + froffset * dy
          srcypix(1) = srcypix(1) + 1
      endif
  endif 

  if yrange(0) LT imgyrange(0) then $
    imgpanel(1) = (imgyrange(0) - yrange(0))/(yrange(1)-yrange(0))
  if yrange(1) GT imgyrange(1) then $
    imgpanel(3) = (imgyrange(1) - yrange(0))/(yrange(1)-yrange(0))

  PLOTDATA:

  if n_elements(position) EQ 0 AND n_elements(panel) EQ 0 $
    AND n_elements(subpanel) EQ 0 then begin
      ;; If neither POSITION nor PANEL/SUBPANEL are given, then plot
      ;; once to set up axes, despite NOAXES
      plot, xrange, yrange, noerase=noerase, /nodata, xstyle=1, ystyle=1, $
        _EXTRA=extra
      mypanel = fltarr(4)
      mypanel([0,2]) = !x.window
      mypanel([1,3]) = !y.window
      imgposition = subcell(imgpanel, mypanel)
  endif else begin
      ;; Construct the plot size from panel info.  Default is full-screen
      if NOT keyword_set(noerase) then erase
      if n_elements(position) GT 0 then begin
          imgposition = subcell(imgpanel, position)
      endif else begin
          if n_elements(panel) EQ 0 then panel=[0.0,0.0,1.0,1.0]
          if n_elements(subpanel) EQ 0 then subpanel = [-1., -1, -1, -1]
          subpanel = defsubcell(subpanel)
          imgpanel = [0., 0., 1, 1]

          imgposition = subcell(subcell(imgpanel, subpanel), panel)
          position = subcell(subpanel, panel)
      endelse
  endelse

  if NOT keyword_set(mynodata) then begin
      myimg = img(srcxpix(0):srcxpix(1), srcypix(0):srcypix(1))
      myimg = reform(myimg, srcxpix(1)-srcxpix(0)+1, srcypix(1)-srcypix(0)+1, $
                     /overwrite)
      tvimage, myimg, position=imgposition, /half_half
  endif

  if NOT keyword_set(noaxes) then begin
      if n_elements(mypanel) GT 0 then begin
          axis, xaxis=0, xtickformat='(A1)', xstyle=1
          axis, xaxis=1, xtickformat='(A1)', xstyle=1
          axis, yaxis=0, ytickformat='(A1)', ystyle=1
          axis, yaxis=1, ytickformat='(A1)', ystyle=1
      endif else begin
          plot, xrange, yrange, /noerase, /nodata, /normal, xstyle=1, $
            ystyle=1, position=position, _EXTRA=extra
      endelse
  endif

  return
end