PRO TVSTACK, IMAGES, W_BORDER=W_BORDER, C_BORDER=C_BORDER, TVSCL=TVSCL, $
	POSITION=POSITION, $
	REVERSE_DIM=REVERSE_DIM, ROTATE_PAR=ROTATE_PAR, SKIP=SKIP

;+
; NAME:
;	TVSTACK
;
; PURPOSE:
;	This procedure displays a multidimensional array (e.g. 3D) as a stack of
;	images lain side by side, and by default scaled between the minimum and
;	maximum values in the array.  This is a method of visualizing 3D data
;	sets, for example reconstructions of tracer uptake in nuclear medicine.
;	For an image measuring NX by NY by NZ, the effect is somewhat the same
;	as typing
;
;		for z=0,nz-1 do tvscl, image(*,*,z), z
;
;	except that they aren't necessarily scaled to the frame maxima, and
;	some options are available through the keywords.
;
; CATEGORY:
;	IMAGE DISPLAY
;
; CALLING SEQUENCE:
;	TVSTACK, IMAGES
;
; INPUTS:
;	IMAGES: The multidimensional array to be displayed.  It may have any
;	number of dimensions but is displayed as if it is 3D.
;	
; KEYWORD PARAMETERS:
;	W_BORDER: Width of the border to be drawn around each frame displayed.
;		If framing is performed, each NX by NY image takes up
;		(NX+2*W_BORDER) pixels in the X direction and (NY+2*W_BORDER)
;		pixels in the Y direction.  Default = NO BORDER.
;	C_BORDER: Color of the border, if the W_BORDER keyword is set.  Default=
;		!d.n_colors/3.
;	TVSCL: If this keyword is set, each 2D image is scaled to its own max
;		and min.  By default, the max and min of the multidimensional
;		array are used.
;	POSITION: The screen position used by the first image frame in the
;		sequence.  Thus, images are positioned as if the command
;
;			for z=0,nz-1 do tvscl, image(*,*,z), position+z
;
;		were used.  Default = 0 (start in upper left corner of th
;		display window.)
;	REVERSE_DIM: which dimension of the image is to be reversed (using the
;		REVERSE function from the USERLIB); 1 means X, 2 means Y.
;		Default = NO REVERSAL.
;	ROTATE_PAR: what parameter to use to rotate the image frames, with the
;		IDL ROTATE function (q.v.) in which 0=no rotation, 1=90 degrees
;		CCW without transposing, and so on.  Defaul=NO ROTATION.
;	SKIP: Use this keyword if you desire to display only some of the images
;		for example every other one, which would be represented by
;		SKIP=2.  Images are positioned as if the command
;
;			for z=0,nz-1,SKIP do tvscl, image(*,*,z), z/SKIP
;
;		were used.  Default = 1, that is, do not skip.  If you specify
;		/SKIP you end up not skipping (!!!) because in this case you
;		specify steps of 1.
;
; OUTPUTS:
;	None; the procedure writes into the current display window.
;
; EXAMPLE:  A 3D gaussian can be displayed as follows:
;
;	n = 50
;	xysq = shift(dist(n)^2,n/2,n/2)	; x^2 + y^2
;	rsq = fltarr(n,n,n)		; array to hold r^2
;	for z=0,n-1 do rsq(0,0,z) = (z-n/2.0)^2 + xysq
;					; r^2 = x^2 + y^2 + z^2
;	fwhm=10.0 & sig=fwhm/(2*sqrt(2*alog(2))) & gauss_3D=exp(-rsq/(2*sig^2))
;					; 3D gaussian of 10 channels FWHM
;	tvstack, gauss_3D, /w_b		; display it
;
; MODIFICATION HISTORY:
; 	Written by:	James Hamill
;			Siemens Medical Systems
;			2501 N. Barrington Rd.
;			Hoffman Estates, IL  60195-7372
;			(708)304-7760
;			hamill@sgi.siemens.com
;			April, 1992
;-


ndim = 3		; expected # of dimensions

sizi = size(images)
if sizi(0) ne ndim then REFORM_FLAG=1 else REFORM_FLAG=0
nx = sizi(1)
ny = sizi(2)
nz = sizi(3)

if n_elements(w_border) eq 0 then w_border=0
if n_elements(c_border) eq 0 then c_border=!d.n_colors/3
if keyword_set(tvscl) eq 0 then tvscl_flag=0 else tvscl_flag=1
if n_elements(position) eq 0 then position=0
if n_elements(skip) eq 0 then skip=1

;  Max and min ...

if tvscl_flag ne 1 then begin
  mn=min(images)
  mx=max(images)
  if mx eq mn then mx=mn+1
endif

;  Set up buffer with the border color all around

xsize = nx + 2*w_border
ysize = ny + 2*w_border
buffer = replicate(c_border,fix(xsize),fix(ysize))

;  Display if the data are to be reformed first ... 

IF REFORM_FLAG THEN BEGIN

  nz = n_elements(images)/nx/ny
  imp = reform(images,nx,ny,nz)

  for z=0,nz-1,skip do begin
    if tvscl_flag eq 1 then $
	buffer(w_border,w_border) = bytscl(imp(*,*,z)) $
    else $
	buffer(w_border,w_border) = (!d.n_colors-1)*(imp(*,*,z)-mn)/(mx-mn)
    if keyword_set(reverse_dim) then buffer=reverse(buffer,reverse_dim)
    if keyword_set(rotate_par) then buffer=rotate(buffer,rotate_par)
    tv,buffer,z/skip+position
  endfor

;  ...otherwise display the data itself.

ENDIF ELSE BEGIN

  for z=0,nz-1,skip do begin
    if tvscl_flag eq 1 then $
	buffer(w_border,w_border) = bytscl(images(*,*,z)) $
    else $
	buffer(w_border,w_border) = (!d.n_colors-1)*(images(*,*,z)-mn)/(mx-mn)
    if keyword_set(reverse_dim) then buffer=reverse(buffer,reverse_dim)
    if keyword_set(rotate_par) then buffer=rotate(buffer,rotate_par)
    tv,buffer,z/skip+position
  endfor

ENDELSE


RETURN

end