PRO colormap, $
   x,o1,o2,o3, NAME=name, $
   MAP_DATA=map_data, $
   MAP_MODE=map_mode, $
   NORMALIZE=normalize, $
   PARAMETER=parameter

; PARAMETER
; two options:
; 1) may contain structure with fields 
;    MAP_MODE
;    MAP_DATA
;    NORMALIZE
; 2) if 2D array it is interpreted as MAP_DATA

; MAP_DATA [2...8,n]
; OPTION 1
; x , r, g, b, gamma[, map_mode]
; assume paramter is 5[+1] x N array

; OPTION 2
; x , r, g, b, gamma_r, gamma_g, gamma_b[, map_mode] 
; assume paramter is 7[+1] x N array

; OPTION 3
; x, v[, gamma[, map_mode]]
; assume paramter is 2[+1[+1]] x N array

; MAP_MODE - default: 0
; 0 RGB
; 1 HSV
; 2 HLS
; 3 YUV 
; 4 YIQ
; 5 YPbPr
; 6 YCbCr 

;; 0: RGB 
;; values are bytes in the range 0 to 255. 

;; YCbCr values are bytes with the range indicated below. 
;; Other values are in floating point. 
;; These color spaces and their ranges are: 

;; 1: HSV 
;; Hue Saturation Value (single cone). 
;; Hue is measured in degrees between 0 and 360, 
;;   with   0 corresponding to red, 
;;        120 to green, and 
;;        240 to blue. 
;; Saturation and Value are floating-point values between 0 and 1. 

;; 2: HLS 
;; Hue Lightness Saturation (double cone). 
;; Also known as HSL and HSI. 
;; Hue is measured in degrees between 0 and 360, 
;;   with   0 corresponding to red, 
;;        120 to green, and 
;;        240 to blue. 
;; Lightness and Saturation are floating-point values between 0 and 1. 

;; 3: YUV 
;; Luminance and two Chrominance. 
;; Y is a floating-point value between 0 and 1, 
;; U is a floating-point value between -0.436 and 0.436, 
;; V is a floating-point value between -0.615 and 0.615. 

;; 4: YIQ 
;; Luminance In-phase Quadrature. 
;; Y is a floating-point value between 0 and 1, 
;; I is a floating-point value between -0.596 and 0.596, 
;; Q is a floating-point value between -0.523 and 0.523. 

;; 5: YPbPr 
;; Luma ChromaBlue ChromaRed (analog). 
;; Y is a floating-point value between 0 and 1, 
;; Pb and Pr are floating-point values between -0.5 to 0.5. 

;; 6: YCbCr 
;; Luma ChromaBlue ChromaRed (digital). 
;; Y is a floating-point value between 16 and 235, 
;; Cb and Cr are floating-point values between 16 and 240. 


; interpret entry i according to color mode
; do interpolation according to colormode and gamma of entry i+1

; NORMALIZE
; 0: take values as is
; 1: valid value range mapped to 0...1
;-1: try automatic detection (default)

@compile_opt

IF N_ELEMENTS(parameter) NE 0 THEN BEGIN
   type=SIZE(parameter,/TYPE)
   IF type EQ 8 THEN BEGIN
      tags=TAG_NAMES(parameter)
      i=(WHERE(tags EQ 'MAP_MODE'))[0]
      IF i NE -1 THEN map_mode=parameter.(i)
      i=(WHERE(tags EQ 'MAP_DATA'))[0]
      IF i NE -1 THEN map_data=parameter.(i)
      i=(WHERE(tags EQ 'NORMALIZE'))[0]
      IF i NE -1 THEN normalize=parameter.(i)
   ENDIF
   dims=SIZE(parameter,/N_DIMENSIONS)
   IF dims EQ 2 THEN map_data=parameter 
ENDIF

IF N_ELEMENTS(map_mode) EQ 0 THEN map_mode=0

c=map_data
dim=SIZE(c,/DIMENSIONS)
IF N_ELEMENTS(dim) NE 2 THEN BEGIN
    PRINT,'[COLORMAP] WRONG DIMENSIONS'
    RETURN
ENDIF

IF (dim[0] LT 2) OR (dim[0] GT 8) THEN BEGIN
    PRINT,'[COLORMAP] WRONG NUMBER OF ENTRIES'
    RETURN
ENDIF


; IF SIZE(x,/n_dimensions) GT 1 then STOP

one=x
one[*]=1

o1=x
o1[*]=0
o2=o1
o3=o1

; deal with template overwrite
n=N_ELEMENTS(c[0,*])
w=WHERE(c[0,*] GT -1.d99)
c=c[*,w]
IF N_ELEMENTS(map_mode) EQ n THEN map_mode=map_mode[w]
IF N_ELEMENTS(normalize) EQ n THEN normalize=normalize[w]

col_min=MIN(c[0,*],MAX=col_max)
y=(x>col_min)<col_max

c=DOUBLE(c[*,SORT(c[0,*])])
m=dim[0]
n=dim[1]

IF (m EQ 4) OR (m EQ 6) OR (m EQ 8) THEN BEGIN 
   map_mode=c[m-1,*]
   c=c[0:m-2,*]
   m-=1
ENDIF

IF m EQ 2 THEN BEGIN
   c=c[[0,1,1],*]
   c[2,*]=1
   m+=1
ENDIF

IF N_ELEMENTS(map_mode) EQ n-1 THEN map_mode=[map_mode[0],map_mode]
IF N_ELEMENTS(map_mode) LT n   THEN map_mode=[map_mode,REPLICATE(map_mode[N_ELEMENTS(map_mode)-1],n-N_ELEMENTS(map_mode))]

IF N_ELEMENTS(normalize) EQ 0 THEN BEGIN 
   normalize=0
   autonorm=1
ENDIF ELSE BEGIN
   autonorm=0
ENDELSE
IF N_ELEMENTS(normalize) EQ n-1 THEN normalize=[normalize[0],normalize]
IF N_ELEMENTS(normalize) LT n   THEN normalize=[normalize,REPLICATE(normalize[N_ELEMENTS(normalize)-1],n-N_ELEMENTS(normalize))]

IF N_ELEMENTS(autonorm) EQ n-1 THEN autonorm=[autonorm[0],autonorm]
IF N_ELEMENTS(autonorm) LT n   THEN autonorm=[autonorm,REPLICATE(autonorm[N_ELEMENTS(autonorm)-1],n-N_ELEMENTS(autonorm))]

w=WHERE(normalize EQ -1)
IF w[0] NE -1 THEN BEGIN
   autonorm[w]=1
   normalize[w]=0
ENDIF

IF m EQ 3 THEN BEGIN
   c=c[[0,1,1,1,2],*]
   m+=2
   FOR i=0,n-1 DO BEGIN
      CASE map_mode[i] OF
         1: BEGIN 
            c[2,i]=1
            c[3,i]=1
         END
         2: BEGIN 
            c[2,i]=0.5
            c[3,i]=1
         END
         3: BEGIN 
            c[2,i]=0
            c[3,i]=0
         END
         4: BEGIN 
            c[2,i]=0
            c[3,i]=0
         END
         5: BEGIN 
            c[2,i]=0
            c[3,i]=0
         END
         6: BEGIN   
; this may conflict with auto-normalization for this color model         
            IF c[1,i] LT 16 THEN BEGIN 
               c[2,i]=0.5
               c[3,i]=0.5
            ENDIF ELSE BEGIN
               c[2,i]=128
               c[3,i]=128
            ENDELSE
         END
         ELSE:
      ENDCASE
   ENDFOR
ENDIF

; try some automatic fixing of wrong ranges.
; (0..1-->0..255 for RBG)
; (0..1-->0..360 for H)
; etc.
;Auto-detect wrong normalization
ii=WHERE(autonorm EQ 1)
IF ii[0] NE -1 THEN BEGIN
   w=WHERE(map_mode[ii] EQ 0)
   IF w[0] NE -1 THEN BEGIN
      IF MAX(c[1:3,ii[w]]) LE 1 THEN normalize[ii[w]]=1
   ENDIF
   w=WHERE((map_mode[ii] EQ 1) OR (map_mode[ii] EQ 2))
   IF w[0] NE -1 THEN BEGIN
      IF MAX(c[1,ii[w]]) LE 1 THEN normalize[ii[w]]=1
   ENDIF
   w=WHERE(map_mode[ii] EQ 3)
   IF w[0] NE -1 THEN BEGIN
      IF (MAX(c[2,ii[w]]) GT 0.436) OR $
         (MAX(c[3,ii[w]]) GT 0.615) $
      THEN normalize[ii[w]]=1
   ENDIF
   w=WHERE(map_mode[ii] EQ 4)
   IF w[0] NE -1 THEN BEGIN
      IF (MAX(c[2,ii[w]]) GT 0.596) OR $
         (MAX(c[3,ii[w]]) GT 0.523) $
      THEN normalize[ii[w]]=1
   ENDIF
   w=WHERE(map_mode[ii] EQ 5)
   IF w[0] NE -1 THEN BEGIN
      IF (MAX(c[2,ii[w]]) GT 0.5) OR $
         (MAX(c[3,ii[w]]) GT 0.5) $
      THEN normalize[ii[w]]=1
   ENDIF
   w=WHERE(map_mode[ii] EQ 6)
   IF w[0] NE -1 THEN BEGIN
      IF MAX(c[1:3,ii[w]]) LT 16 THEN normalize[ii[w]]=1
   ENDIF
ENDIF
FOR i=0,n-1 DO BEGIN
   IF normalize[i] EQ 1 THEN BEGIN
      CASE map_mode[i] OF
         0: c[1:3,i]*=255.D0 
         1: c[1,i]*=360.D0 
         2: c[1,i]*=360.D0 
         3: BEGIN
            c[2:3,i]-=0.5
            c[2,i]*=0.436*2
            c[3,i]*=0.615*2
         END
         4: BEGIN
            c[2:3,i]-=0.5
            c[2,i]*=0.596*2
            c[3,i]*=0.523*2
         END
         5: c[2:3,i]-=0.5
         6: BEGIN 
            c[1,i]*=(235-16)
            c[2:3,i]*=(240-16)
            c[1:3,i]+=16
         END
      ENDCASE
   ENDIF
ENDFOR

c0=c[*,0]
FOR i=1,n-1 DO BEGIN
   IF map_mode[i-1] NE map_mode[i] THEN BEGIN
      cc1=c0[1]
      cc2=c0[2]
      cc3=c0[3]
      color_translate,cc1,cc2,cc3,map_mode[i-1],map_mode[i]
      c0[1:3]=[cc1,cc2,cc3]
   ENDIF

   c1=c[*,i]
   w=WHERE((y GE c0[0]) AND (y LE c1[0]))
   IF w[0] NE -1 THEN BEGIN
      dc=c1-c0
      cc=((y[w]-c0[0])/dc[0])
      IF m EQ 5 THEN BEGIN
         cc^=c1[4]
         o1[w]=c0[1]+cc*dc[1]
         o2[w]=c0[2]+cc*dc[2]
         o3[w]=c0[3]+cc*dc[3]
      ENDIF ELSE BEGIN
         o1[w]=c0[1]+cc^c1[4]*dc[1]
         o2[w]=c0[2]+cc^c1[5]*dc[2]
         o3[w]=c0[3]+cc^c1[6]*dc[3]
      ENDELSE
      IF map_mode[i] NE 0 THEN BEGIN
         cc0=o1[w]
         cc1=o2[w]
         cc2=o3[w]
         color_translate,cc0,cc1,cc2,map_mode[i],0
         o1[w]=cc0
         o2[w]=cc1
         o3[w]=cc2
      ENDIF

   ENDIF
   c0=c1
ENDFOR

o1=BYTE(ROUND(o1)<255>0)
o2=BYTE(ROUND(o2)<255>0)
o3=BYTE(ROUND(o3)<255>0)

END