pro Gr_display_img, zin, xin, yin, range=range, colour_range=colour_range, $
                   pixel_size=pixel_size

;+
; GR_DISPLAY_IMG
;	Colour/greyscale image display for GRAFFER
;
; Usage:
;	gr_display_img, zin, xin, yin, range=range, $
;	colour_range=colour_range pixel_size=pixel_size
;
; Arguments:
;	zin	float	input	The data to be displayed
;	xin	float	input	The X coordinates of the data.
;	yin	float	input	The Y coordinates of the data.
;
; Keywords:
;	range	float	input	The range from "black" to "white"
;	colour. int	input	The range of colour indices to use.
;	pixel.  float	input	For devices with scalable pixels, the
;				size of a displayed pixel.
;
; History:
;	Original: 10/12/96; SJT
;	Add code to clip to the viewport: 12/12/96; SJT
;-

if (!D.flags and 1) then begin  ; PS or similar with scalable pixels
    if (not keyword_set(pixel_size)) then pixel_size = 0.5 ; default
                                ; 0.5 mm pixels
    scfac = 10./([!D.x_px_cm, !D.y_px_cm] * pixel_size)
endif else scfac = [1, 1]

;	Select out those parts which are within the viewport

if (!X.type eq 1) then  $
  locsx = where(xin ge 10^!X.crange(0) and xin le 10^!X.crange(1), nx) $
else locsx = where(xin ge !X.crange(0) and xin le !X.crange(1), nx)
if (!Y.type eq 1) then  $
  locsy = where(yin ge 10^!Y.crange(0) and yin le 10^!Y.crange(1), ny) $
else locsy = where(yin ge !Y.crange(0) and yin le !Y.crange(1), ny)

if (nx eq 0 or ny eq 0) then return ; Image is wholly outside the VP

x = xin(locsx)
y = yin(locsy)
z = zin(locsx, *)
z = z(*, locsy)

xrange = [(mnx = min(x, max = mxx)), mxx]
yrange = [(mny = min(y, max = mxy)), mxy]

corners = convert_coord(xrange, yrange, /data, /to_device)
corners = round(corners(0:1, *)*scfac(*, intarr(2)))
dvsize = (corners(*, 1)-corners(*, 0))
cmsize = (corners(*, 1)-corners(*, 0))/([!D.x_px_cm, !D.y_px_cm]*scfac)
cmll = corners(*, 0)/([!D.x_px_cm, !D.y_px_cm]*scfac)

xd = convert_coord((findgen(dvsize(0))+corners(0, 0)) / scfac(0), $
                   fltarr(dvsize(0)), /device, /to_data)
xd = reform(xd(0, *))
yd = convert_coord(fltarr(dvsize(1)),  $
                   (findgen(dvsize(1))+corners(1, 0)) / scfac(1), /device, $
                   /to_data)
yd = reform(yd(1, *))

sz = size(z)
xx = interpol(findgen(sz(1)), x, xd)
yy = interpol(findgen(sz(2)), y, yd)

zz = bilinear(z, xx, yy)

if (range(0) eq range(1)) then zrange = [min(z, max = mxz), mxz] $
else zrange = range

if (not keyword_set(colour_range)) then colour_range = [16, $
                                                        !D.n_colors-1]

locs = where(finite(zz) eq 0, nnan)

if (nnan gt 0) then zz(locs) = range(0)-1

if (!d.flags and 1) then $
  tv, bytscl(zz, min = zrange(0), max = zrange(1), top = $
           colour_range(1)-colour_range(0)) + colour_range(0), $
  cmll(0), cmll(1), xsize = cmsize(0), ysize = cmsize(1), /centi $
else $
  tv, bytscl(zz, min = zrange(0), max = zrange(1), top = $
             colour_range(1)-colour_range(0)) + colour_range(0), $
  corners(0, 0), corners(1, 0)

end