; Time-stamp: <97/12/26 15:41:38 dph> 
;             MIT Directory: ~dph/h1/ASC/TG/Anal/3color_spec/
;                            /wiwaxia/d4/ASC/lib/IDL
; File: mk_color_image
; Author: D. Huenemoerder
;         Original version: 971226 
;
; (this header is ~dph/libidl/time-stamp-template.el)
; to auto-update the stamp in emacs, put this in your .emacs file:
;      (add-hook 'write-file-hooks 'time-stamp)
;====================================================================

PRO mk_color_image, x, y, xy_bin, color_vector, image, x_axis, y_axis, rgb_table
;+
; NAME:
; mk_color_image
;
;
; PURPOSE:
;   create an image and color-code it according to the mean of a
;   specified parameter in each x,y bin.
;
;
; CATEGORY:
;   image idsplay
;
;
; CALLING SEQUENCE:
;   mk_color_image, x, y, xy_bin, color_vector, image, x_axis, y_axis, rgb_table
;
; 
; INPUTS:
;   x = [np] =  input vectors to be binned into an image.
;   x = [np] =  input vectors to be binned into an image.
;   xy_bin = [dx,dy] = bin-sizes on x,y for image
;   color_vector = [np] = vector commensurate w/ x,y for calculation
;     of color for x,y bins.
;
;
; OPTIONAL INPUTS:
;    NONE
;
;       
; KEYWORD PARAMETERS:
;    NONE
;
;
; OUTPUTS:
;
;   image = [np/xbin,ny/ybin] = x,y binned intensity image
;   x_axis = [np] = x coordinate array for image
;   y_axis = [np] = y coordinate array for image
;   rgb_table = [np,3] = color table, in r,g,b system, encoding
;          <color_vector> in each x,y bin.
;
; OPTIONAL OUTPUTS:
;
;   NONE
;
; COMMON BLOCKS:
;
;  NONE
;
; SIDE EFFECTS:
;
;   NONE
;
; RESTRICTIONS:
;
;   x,y,color_vector must be same lengths and commensurate.
;
;   Saturation is pegged at 1.0.
;   Hue table is set to 0-359.
;
; PROCEDURE:
;   use make_image to form image in x,y, returning axes and
;   reverse_indices.  
;   Use reverse_indices to bin the color_vector into an image.
;   Use x,y image for intensity, use color_vector_image for hue, set
;   saturation to 1, then convert HSV planes to RGB table.
;   Color_quantize to RGB system.
;
;
; EXAMPLE:
;   marx_dir = './marx'
;   yp = read_marx_file(marx_dir+'/ypos.dat')
;   zp = read_marx_file(marx_dir+'/zpos.dat')
;   detector = read_marx_file(marx_dir+'/detector.dat')
;   cvec = read_marx_file(marx_dir+'/pha.dat')
;   mk_color_image, yp, zp,[8,2]*0.024, cvec, $
;    im_out, x_out, y_out, rgb_out
;   tvlct,rgb
;   tv,im_out
;
; MODIFICATION HISTORY:
;
;   971226 (dph) original procedure, based on 971212 prototype.
;-
   
;;;;;;;;;;;;;;;;;; USAGE string array: ;;;;;;;;;;;;;;;;

u_sarray = [$
     ' PURPOSE:', $
     'mk_color_image, x, y, xy_bin, color_vector, image, x_axis, y_axis, rgb_table', $
     '', $
     'color-code an image according to the mean of a specified parameter in each x,y bin.', $
     'INPUTS:  ', $
     '  x = [np] =  input vectors to be binned into an image.', $
     '  y = [np] =  input vectors to be binned into an image.',$
     '  xy_bin = [dx,dy] = bin-sizes on x,y for image.', $
     '  color_vector = [np] = vector commensurate w/ x,y for calculation of color for x,y bins.', $
     'OUTPUTS:  ', $
     '  image = [(xmax-xmin)/xbin,(ymax-ymin)/ybin] = x,y binned intensity image', $
     '  x_axis = [np] = x coordinate array for image', $
     '  y_axis = [np] = y coordinate array for image', $
     '  rgb_table = [np,3] = color table, in r,g,b system, encoding', $
     '       <color_vector> in each x,y bin.'$
    ]

   np = n_params()
   IF np NE 8 THEN BEGIN
      print, ''
      message, /inf, 'Incorrect number of parameters.'
      print, 'USAGE:  '+ u_sarray, format = '(a)'
      print, ''
      return
   ENDIF
   
   xra = [min(x, max = tmp), tmp]
   yra = [min(y, max = tmp), tmp]
   
   im = make_image(x, y, xra=xra, yra=yra, xbin=xy_bin[0], ybin=xy_bin[1], $
                   xax=x_axis, yax=y_axis, $
                   index_list=idx, reverse_indices=rev_idx)
   
   im_colors = im*0
   l_colors = color_vector(idx)              ; color_vector's for ypos,zpos indices used in image.

   lnz = where(im GT 0)   ;; only need to colorize non-zero intensity image bins
   nlnz = n_elements(lnz)

;   message, /inf, 'Begin making <color_vector> image' ;;;; debug

   FOR ii = 0, nlnz-1 DO BEGIN  ;;; loop over non-zero image bins

      pixnum = lnz(ii)  ;;; serial bin number

      IF rev_idx(pixnum) NE rev_idx(pixnum+1) THEN BEGIN   
         ;; if there is a contribution (this should always be true
         ;; since we already selected non-zero bins)

         plist = rev_idx( rev_idx(pixnum):rev_idx(pixnum+1)-1 ) ; who contributed?
         
                                ;  IF n_elements(plist) GT 0 THEN
                                ;  BEGIN ; must be true to get here 
         im_colors(pixnum) = avg(l_colors(plist))  

      ENDIF ELSE BEGIN
         message,'FATAL error: no contribution to this non-zero bin!'
      ENDELSE
   
   ENDFOR

;   message, /inf, 'Converting color image' ;;;; debug
                                ; make H,S,V planes, then
                                ; color_convert, and color_quan

   im_H = float(im_colors) / max(im_colors) * 359. ; all Hues (0-359).
   im_S = float(im NE 0)        ; "saturation" - make it 1.0 or 0
                                ; Caller should impose a threshold on the
                                ; brightness, since the dynamic range
                                ; could be large.  
   im_V = float(im)/max(im)     ; "value" - the intensity; 0.0-1.0
   thresh = max(im_V)           ;
   im_V = float(im < thresh) / thresh ; "value" - the intensity; 0.0-1.0
                                ;this clips out a single large spike
                                ;which dominates the intensity scale.

   color_convert, im_H, im_S, im_V, im_r, im_g, im_b, /hsv ; returns r,g,b images.
   imtrue=[ [[im_r]],[[im_g]],[[im_b]]]  

   image = color_quan(imtrue, 3, r, g, b, /dither, colors = 256) ; convert to 8-bit.
   
   rgb_table = [ [r], [g], [b]]


END