;+
; NAME:
;	RGBO_Project
; PURPOSE:
;	Merges content of z-buffer with color and opacity
;	array of volume data set
; CATEGORY:
;
; CALLING SEQUENCE:
	FUNCTION RGBO_Project, zvol, rgbo, $
		zpix		= zpix			, $
		zbuf		= zbuf			, $
		ctrgb		= ctrgb 		, $
		background	= background	, $
		bad_data	= bad_data		, $
		bad_opacity = bad_opacity	, $
		low_data	= low_data		, $
		low_opacity = low_opacity	, $
		high_data	= high_data 	, $
		high_opacity= high_opacity	, $
		matte		= matte
; INPUTS:
;	zvol			array[N,M,L]; type: byte
;						byte array for volume data
;	rgbo			array[256,4]; type; float
;						rgbo vector; rgbo[*,0]=o*r, rgbo[*,1]=o*g, rgbo[*,2]=o*b
;						rgbo[*,3] = o*255. o = opacity (0 <= o <= 1), r,g,b = red,green
;						blue color indices (0 <= r,g,b <= 255).
; OPTIONAL INPUT PARAMETERS:
;	zpix=zpix		array[n,m]; type: byte
;						Z-buffer content (the Z-buffer iis a 256-color device)
;					array[n,m,3]; type: byte
;						Z-buffer content as true color image
;	zbuf=zbuf		array[n,m]; type: integer
;						Z-buffer depth array
;	background=background
;					scalar; type: byte
;						background color index in Z-buffer
;						If not specified then !p.background is used.
;					array[3]; type: byte
;						(only if zpix is true color image)
;						background color rgb values in Z-buffer
;	ctrgb=ctrgb		array[256,3]; type: byte
;						(only if zpix is 256-color image)
;						rgb values for color table to be used to interpret zpix
;						If not specified then the current color table is retrieved
;						with the IDL tvlct function.
;	matte=matte		array[N,M,L] or array[N,M,L,4]
;						matrix of values between 0 and 1 to be applied to
;						zvol after conversion to rgbo values and before projecting onto
;						the image plane. This can be used e.g. to implement depth cueing.
;						An array[N,M,L] is applied to each channel (r,g,b,o).
; OUTPUTS:
; OPTIONAL OUTPUT PARAMETERS:
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	ArrayLocation, SuperArray
; SEE ALSO:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; MODIFICATION HISTORY:
;	JUL-2002, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-


zbuffer = n_elements(zpix) NE 0 AND n_elements(zbuf) NE 0

nvol = size(zvol, /dim)			; Size of volume data

CASE zbuffer OF
0: npix = [!d.x_size, !d.y_size]
1: npix = size(zpix, /dim)		; Size of z-buffer
ENDCASE

true_color = n_elements(npix) EQ 3

CASE true_color OF
0: BEGIN
	IF n_elements(ctrgb) EQ 0 THEN BEGIN
		tvlct, /get, red, green, blue
		ctrgb = [ [red], [green], [blue] ]
	ENDIF

	InitVar, background, !p.background
	background_rgb = reform(ctrgb[background,*])
END
1: background_rgb = background
ENDCASE

nchan = 4

; Convert the byte values in 'zvol' to rgbo quadruplets
; The result is a 2-dim float array[nvol[0],nvol[1],nvol[2],4]

bvol = reform(rgbo[zvol,*], [nvol,nchan])

; Apply matte

IF n_elements(matte) NE 0 THEN BEGIN

	; The matte is applied to r,g,b and o channels.
	; If the matte is a 3d-array then each of the r,g,b,o
	; channels is multiplied by it. If it is a 4D array
	; already then a straignt multiplication is done.

	sz = size(matte)
	CASE (where( sz[1:3] NE nvol ))[0] OF
	-1  : BEGIN

		CASE sz[0] OF
		3: bvol = bvol*SuperArray(matte,nchan,/trail)
		4: bvol = bvol*matte
		ENDCASE

	END
	ELSE: message, /info, 'matte not same size as volume data; matte ignored'
	ENDCASE

ENDIF


; bvol is the rgb volume array, still at the resolution used for the
; remote view calculations, i.e. not yet at the resolution of the display.
; Create a 'transparency' array (1 minus the opacity). The trailing
; dimension of 3 is added to make processing with the rgb values easier.

tvol = SuperArray(1.0-bvol[*,*,*,3]/255.0, 3, /trail)

; Do the projection to the image plane. Combine planes in the z-direction
; processing from the back of the volume (i=0) to the front (i=nvol[2]-1).
; The transparency at plane i is used to add in the effects of all planes
; behind it (i' < i).

; Create true color image array (we create a float array for now).
; Initialize the array to white, rgb = [255,255,255]


zout = make_array(type=IsType(bvol), dimension=[nvol[0:1],3])
zout[*] = SuperArray(background_rgb, nvol[0]*nvol[1], /lead)

FOR i=0,nvol[2]-1 DO zout = bvol[*,*,i,0:2]+tvol[*,*,i,*]*zout

zout = reform(zout, /overwrite)					; Remove dummy 3rd dimension of 1

; Rebin to zbuffer dimension 'npix'

IF npix[0] NE nvol[0] OR npix[1] NE nvol[1] then	$
	zout = congrid(zout, npix[0], npix[1], 3, /minus_one)


; zout at this point contains the true-color image for the volume data alone.
; We still need to merge the content of the z-buffer. To do this we only need
; to update those pixels which are part of an object in the z-buffer (i.e.
; where the depth array zbuf > zmin. We also need to include pixels with
; zero depth that do not have the background color.
; The z-buffer depth array is a 16-bit (short) signed integer array
; According to IDL manual the depth has a range from -32765 to +32765
; (corresponding to normalized coordinates 0.0 to 1.0) and is initialized to
; -32765 (normalized coordinate 0.0).

zmin = -32765.0										; Min and max possible values in z-buffer
zmax = +32765.0


; Get a list of pixels in the npix[0] x npix[1] image that need to be updated.

CASE zbuffer OF
0: nbuf = 0
1: BEGIN
	CASE true_color OF
	0: pbuf = where(zbuf GT zmin OR zpix NE background, nbuf)
	1: pbuf = where(zbuf GT zmin OR zpix[*,*,0] NE background_rgb[0] OR		$
									zpix[*,*,1] NE background_rgb[1] OR		$
									zpix[*,*,2] NE background_rgb[2], nbuf)
	ENDCASE
END
ENDCASE


IF nbuf NE 0 THEN BEGIN								; Something in the zbuffer

	; Convert linear index to index pairs into npix[0] x npix[1]

	pvol = ArrayLocation(pbuf, dim=npix[0:1])

	; Determine the index pairs into nvol[0] x nvol[1] voxels of the input volume data

	pvol = round( pvol* ( ((nvol[0:1]-1.0)/(npix[0:1]-1.0))#replicate(1,nbuf) ) )

	; Convert to linear index into nvol[0] x nvol[1] voxels

	pvol = ArrayLocation(pvol, dim=nvol[0:1], /onedim)

	; pvol is a list of nbuf indices 0 <= pvol <= nvol[0]*nvol[1].
	; We need indices for the corresponding pixels in all z-layers (nvol[2])
	; and all 4 color/opacity channels.

	pvol = pvol#replicate(1.0,nvol[2]*nchan)+replicate(1.0,nbuf)#(nvol[0]*nvol[1]*lindgen(nvol[2]*nchan))
	pvol = reform(pvol, nbuf*nvol[2]*nchan, /overwrite)

	; Convert all rgbo quadruplets from 'bvol'. We now collected nbuf columns of
	; nvol[2] voxels for all 4 rgbo channels.

	bvol = reform( bvol[pvol], nbuf, nvol[2], nchan )


	; Scale the depth for the nbuf pixels to the range [0,nvol[2]-1]
	; Then replicate the depths over the entire column of nvol[2].

	zdepth = (nvol[2]-1)/(zmax-zmin)*(zbuf[pbuf]-zmin)	; array[nbuf]
	zdepth = SuperArray(zdepth, nvol[2], /trail)		; array[nbuf,nvol[2]]

	; Set up the positions in the z-direction for all nbuf pixels

	pvol = SuperArray(lindgen(nvol[2]), nbuf, /lead)	; array[nbuf,nvol[2]]

	; Find the voxels that are behind the depth value in the z-buffer.
	; Set opacities of voxels behind something in the z-buffer to zero.

	pvol = where(zdepth GT pvol, n)

	; pvol is a list of n indices into nbuf*nvol[2] voxels. We need the
	; corresponding elements for all 4 rgbo channels.

	IF n GT 0 THEN BEGIN

		pvol = pvol#replicate(1.0,nchan)+replicate(1.0,n)#(nbuf*nvol[2]*lindgen(nchan))
		pvol = reform(pvol, n*nchan, /overwrite)

		; Set the opacities behind the z-buffer to zero (i.e. make them fully
		; transparent).

		bvol[pvol] = 0.0

	ENDIF

	; Do the projection to the image plane. The background is set to the
	; color from the z-buffer. Since the opacities are zero (and hence the
	; transparency is one) up to the point where the object in the z-buffer
	; is located, this color is not diluted up to that point.

	tvol = SuperArray(1.0-bvol[*,*,3]/255.0, 3, /trail)


	zout = reform(zout, npix[0]*npix[1],3, /overwrite)

	; If zpix is input as a 2D-array of color indices we need to check the color
	; table to find the corresponding r,g,b values.

	CASE true_color OF
	0: zout[pbuf,*] = ctrgb[zpix[pbuf],*]						; Use color table to get rgb values
	1: zout[pbuf,*] = (reform(zpix,npix[0]*npix[1],3))[pbuf,*]
	ENDCASE

	FOR i=0,nvol[2]-1 DO zout[pbuf,*] = bvol[*,i,0:2]+tvol[*,i,*]*zout[pbuf,*]

	zout = reform(zout, npix[0],npix[1],3, /overwrite)

ENDIF


RETURN, byte(zout)  &  END