; Copyright (c) 1998-2013 A.P. Hitchcock  All rights reserved
;+
;NAME:
;	DISP3COL
;
;LAST CHANGED: ----------------------------------- 24-Apr-13
;
;PURPOSE:
;	This procedure generates a composite image from any 3 images using R,G,B mapping
; with user selectable absolute or relative scaling of the components.
; Interpolates to comon mesh prior to RGB generation
;
;CATEGORY:
;	AXIS: image analysis. Also executes stand-alone
;
;CALLING SEQUENCE:
; for procedures:
;	DISP3COL,im1, im2, im3, r=r, b=b, g=g, axis_on=axis_on, channel=channel
;
;CALLED FROM AXIS:
;	->Images->RGB image
;
;INPUTS:
;	All input parameters are passed as keywords.
;
;KEYWORDS:
;	IM1 - 2d array with first image
;	IM2 - 2d array with second image
;	IM3 - 2d array with third image
;	R - red color table
;	G - green color table
;	B - blue color table
;	axis_on - indicates called from AXIS
;	channel - indicates number of channels in image (AXIScomposite uses 3)
;   scalebar  - if on, plot a bar (does not work right - 23-Apr-2013)
;
;OUTPUTS:
;	The RGB-coded combination image is placed in buffer 0
; and saved as a tif file. Default format is 3-channel, byte pixels.
;
; RESTRICTIONS
;  (a) Since write as byte pixels, can be losss of intensity scale resolution
; This is ARTIFICIAL since IDL can write integer or reals (AXIS images are reals)
;  (b) Each image is INDEPENDENTLY byte scaled - relative intensities of
; the colors do NOT reflect relative intensities in the 3 component images
;
;COMMON BLOCKS:
;	@AXIS_COM	standard set of common blocks
;
;MODIFICATION HISTORY:
; ( 9-jun-99 aph) first developed from tv3d.pro (Eli's stereo projection code)
; (25-jun-99 aph) adapt to single (image, spectra) color table
; (25-feb-00 aph) AXIS standard documentation
; (16-apr-00 aph) get dimensions correct and set up to display correctly in AXIS
;               remove intermediate writing of bytscl images -> no loss of resolution
; (30-nov-00 cgz) Replaced call of wdivide with explicit !p.multi statement
; (30-nov-00 cgz) Replaced call of nwin with explicit window statement
; (06-dec-00 aph) autoscale option added
; (02-may-04 aph) replace get_num with dialog_message (default is autoscale)
; (22-feb-06 aph) interpolate to same (x,y) pixels
; (19-Mar-13 aph) add scale bar
; (24-Apr-13 aph) remove scale bar - not accurate in distance; color problems
;-

PRO disp3col,im1, im2, im3, r=r, b=b, g=g, axis_on=axis_on, channel=channel, scalebar = scalebar

@axis_com
on_error,2

if n_tags(im1) EQ 0 then im1=axb_load()
if n_tags(im2) EQ 0 then im2=axb_load()
if n_tags(im3) EQ 0 then im3=axb_load()

; -------------interpolate to image with largest dimensions
ts_red = size(im1.d)
ts_grn = size(im2.d)
ts_blu = size(im3.d)
x_vals = [ts_red(1),ts_grn(1), ts_blu(1)]
y_vals = [ts_red(2),ts_grn(2), ts_blu(2)]
xsize = max(x_vals)
if total(x_vals)/3.-xsize LT 0. then axis_log, 'unequal X-array dimensions ' + $
   string(x_vals(0)) + string(x_vals(1)) + string(x_vals(2))
ysize = max(y_vals)
if total(y_vals)/3.-ysize LT 0. then axis_log, 'unequal Y-array dimensions ' + $
	string(y_vals(0)) + string(y_vals(1)) + string(y_vals(2))
mesh = [xsize, ysize]
im1=ax_interp(im1, mesh=mesh)
im2=ax_interp(im2, mesh=mesh)
im3=ax_interp(im3, mesh=mesh)
ri=im1.d
gi=im2.d
bi=im3.d

if not keyword_set(channel) then channel = 3  ; SET default to true color

; -------- let user decide what vertical scaling to use -----------
if dialog_message('Autoscale each component?',/question, title ='RGB color composite maps') $
	  EQ 'Yes' then autoscale = 1 else autoscale = 0

current_device=!D.NAME
if keyword_set(axis_on) then WIDGET_CONTROL,/hourglass
img_size=[n_elements(im1.x), n_elements(im1.y)]
if not keyword_set(axis_on) then begin
;	nwin, 0, img_size	;generate a window if not in axis
	window, 0, xsize=img_size[0], ysize=img_size[1]
;	wdivide,0,0
	!p.multi=[0, 1, 1]
endif
if NOT keyword_set(r) then begin
	if n_elements(ax_top_color_index) EQ 0 then ax_top_color_index = 256
	r = indgen(ax_top_color_index) & g = indgen(ax_top_color_index) & b = indgen(ax_top_color_index) ; make linear color tables for R G B
endif

if channel EQ 3 then begin
	image = bytarr(3,img_size(0), img_size(1))   ; NB IDL can write integer or reals
; but Paint shop pro only deals with bytes
	if autoscale EQ 1 then begin
		image(0,*,*) = bytscl(ri)
		image(1,*,*) = bytscl(gi)
		image(2,*,*) = bytscl(bi)
	endif else begin
		rgball=[ri,gi,bi]
		all_min = min(rgball,max=all_max)
		slope = 255./(all_max-all_min)
		image(0,*,*) = byte(slope*(ri-all_min))
		image(1,*,*) = byte(slope*(gi-all_min))
		image(2,*,*) = byte(slope*(bi-all_min))
	endelse
endif else begin
;	print, 'Old code used color_quant to merge R,G,B image information'
;	print, 'removed 16-apr-00 - aph; sum of R, G, B signals only !'
	image = bytarr(img_size(0), img_size(1))
	if autoscale EQ 1 then begin
		image = bytscl(ri + gi + bi)
	endif else begin
		rgball=[ri,gi,bi]
		all_min = min(rgball,max=all_max)
		slope = 255./(all_max-all_min)
		image = byte(slope*(ri + gi + bi-3*all_min))
	endelse
endelse

if keyword_set(axis_on) then begin
	wset,MainImg
	title = 'R='+im1.dl+' G='+im2.dl+' B='+im3.dl
	x = im1.x
	y = im1.y
  	s = {t:'3d', x:x, y:y, d:image, r: R, g: G, b: B, xl:im1.xl, yl:im1.yl, dl: title}
  	CurBuf=0
	HANDLE_VALUE, Data(CurBuf), s, /SET
	PlotBuf, CurBuf
endif else begin
	tvlct, r,g,b
	if channel EQ 3 then tv, image, /true else tvscl, image
endelse

SetGraf, 'MainImg', /set

if keyword_set(scalebar) then begin
;	ax_color  ; , /show
;	device,Get_Visual_Name=thisName, Get_Visual_Depth=thisDepth
;	print, 'ax_white_color_index = ', ax_white_color_index
	length = [Xrng(1)-Xrng(0)]
	test = alog10(length/10.)-floor(alog10(length/10.))
	bar = 2*10^(alog10(length/10.)-test)
	if test GT 0.3 then bar = 5*10^(alog10(length/10.)-test)
	if test GT 0.7 then bar = 10*10^(alog10(length/10.)-test)
	xlo = bar_pos(0)*length+xrng(0)
	xhi = xlo + bar*1.16
	ybar = bar_pos(1)*(Yrng(1)-Yrng(0))+Yrng(0)
	oplot, [xlo,xhi], [ybar,ybar], thick = 7, color = ax_white_color_index
;	if bar GE 1.0 then $
;	   xyouts, xlo, ybar + 5, string(bar, format = '(I4," um")'), $
;	   color = ax_white_color_index, /device, charthick = 1.0, charsize = 1.6 $
;	  else xyouts, xlo, ybar + 5, string(1000*bar, format = '(I4," nm")'), $
;	   color = ax_white_color_index, /device, charthick = 1.0, charsize = 1.6
endif

file = pickfile2(/write, FILTER='*.tif',/lpath, DefPath=DefPath)
if file NE '' then begin
	compress = 0
	if keyword_set(channel) then begin
		if channel EQ 3 then begin
			write_tiff,file, reverse(image,3), compression=compress, Red=r, Green=g, Blue=b
		endif else write_tiff, file, reverse(image,2), compression=compress, Red=r, Green=g, Blue=b
	endif else write_tiff, file, reverse(image,2), compression=compress, Red=r, Green=g, Blue=b
	axis_log, 'RBG composite stored as ' + file
endif

end

