R/colorscience.R

Defines functions saturationCIECAM02 saturationCIELABEvaLubbe saturationCIELAB saturationCIELUV xyChromaticitiesVos1978 PhotoYCC2RGB RGB2PhotoYCC XYZ2RxRyRz RxRyRz2XYZ spectra2ISObrightness sprague XYZMoonSpencer1945 XYZTannenbaum1974 emittanceblackbodyPlanck wlnm2xyz wlnm2XYZ kelvin2xy XYZ2CCT.Robertson xy2CCT.HernandezAndres xy2CCT.McCamy StearnsStearnscorrection zFit_1931 yFit_1931 xFit_1931 heuristic.wlnm2RGB Adjust huedegree CIEluminanceY2NCSblackness MunsellV2relativeLuminanceY deltaECMC deltaE2000 deltaE1994 deltaE1976 compuphaseDifferenceRGB InvCompand Compand Luv2XYZ Lab2XYZ XYZ2xyY xyY2XYZ XYZ2Luv XYZ2Lab Y2MunsellV MunsellV2Y Y2MunsellVtable1D1535 DIN99toCIELab CIELabtoDIN99 spectra2CRIGAIFSCI spectra2CCT createIsoTempLinesTable spectra2XYZ spectra2lux LMS2RGB RGB2LMS LMS2XYZ XYZ2LMS LSLM2RGB RGB2LSLM WestlandBlacknessIndex CIE1976uvSaturation CIE1976hueangle CIE1976chroma RGB2hue PreucilPercentHueError PreucilPercentGreyness PreucilAngle RGB2XYZ XYZ2RGB DominantWavelength CIETint GanzGrieser.Tint GanzGrieser.Whiteness Hunter60.WhitenessIndex Taube60.Whiteness Stensby68.Whiteness Berger59.Whiteness CIE.Whiteness ASTM.E313.Whiteness ASTM.E313.YellownessIndex ASTM.D1925.YellownessIndex DIN6167.YellownessIndex LCHuv2Luv Luv2LCHuv LCHab2Lab Lab2LCHab CCT2XYZ NickersonColorDifference HunterLab2XYZ XYZ2HunterLab CMYK2CMY CMY2CMYK CMY2RGB RGB2CMY HSV2RGB RGB2HSV Hue.2.RGB HSL2RGB RGB2HSL YUV2RGB RGB2YUV YIQ2RGB RGB2YIQ YCbCr2RGB RGB2YCbCr YPbPr2RGB RGB2YPbPr LEF2RGB RGB2LEF BVR2XYZ XYZ2BVR CIE1960UCS2CIE1964 CIE1960UCS2xy Yxy2CIE1960UCS CIE1931xy2CIE1976uv CIE1931XYZ2CIE1976uv CIE1931xy2CIE1960uv CIE1931XYZ2CIE1960uv CIE1931XYZ2CIE1931xyz CIE1976uv2CIE1960uv CIE1976uv2CIE1931xy Yuv2xy Yuv2XYZ Yxy2Yuv XYZ2Yuv Yuv2Luv Luv2Yuv dklCart2dkl dkl2dklCart rgb2dklV rgb2dklCart dklCart2rgb dkl2rgb footcandle2candela.steradian.sqmeter footcandle2watt.sqcentimeter footcandle2lux makeChromaticAdaptationMatrix XYZ2xyY LAB2LUV LUV2LAB chromaticity.diagram.color.fill Maxwell.triangle.color.fill Maxwell.triangle.color chromaticity.diagram.color Maxwell.triangle chromaticity.diagram XYZtoRGB LMS2DKL CheckColorLookup ColorBlockFromMunsell MunsellSpecToHVC

Documented in Adjust ASTM.D1925.YellownessIndex ASTM.E313.Whiteness ASTM.E313.YellownessIndex Berger59.Whiteness BVR2XYZ CCT2XYZ CheckColorLookup chromaticity.diagram chromaticity.diagram.color chromaticity.diagram.color.fill CIE1931xy2CIE1960uv CIE1931xy2CIE1976uv CIE1931XYZ2CIE1931xyz CIE1931XYZ2CIE1960uv CIE1931XYZ2CIE1976uv CIE1960UCS2CIE1964 CIE1960UCS2xy CIE1976chroma CIE1976hueangle CIE1976uv2CIE1931xy CIE1976uv2CIE1960uv CIE1976uvSaturation CIELabtoDIN99 CIEluminanceY2NCSblackness CIETint CIE.Whiteness CMY2CMYK CMY2RGB CMYK2CMY ColorBlockFromMunsell Compand compuphaseDifferenceRGB createIsoTempLinesTable deltaE1976 deltaE1994 deltaE2000 deltaECMC DIN6167.YellownessIndex DIN99toCIELab dkl2dklCart dkl2rgb dklCart2dkl dklCart2rgb DominantWavelength emittanceblackbodyPlanck footcandle2candela.steradian.sqmeter footcandle2lux footcandle2watt.sqcentimeter GanzGrieser.Tint GanzGrieser.Whiteness heuristic.wlnm2RGB HSL2RGB HSV2RGB Hue.2.RGB huedegree Hunter60.WhitenessIndex HunterLab2XYZ InvCompand kelvin2xy Lab2LCHab LAB2LUV Lab2XYZ LCHab2Lab LCHuv2Luv LEF2RGB LMS2DKL LMS2RGB LMS2XYZ LSLM2RGB LUV2LAB Luv2LCHuv Luv2XYZ Luv2Yuv makeChromaticAdaptationMatrix Maxwell.triangle Maxwell.triangle.color Maxwell.triangle.color.fill MunsellSpecToHVC MunsellV2relativeLuminanceY MunsellV2Y NickersonColorDifference PhotoYCC2RGB PreucilAngle PreucilPercentGreyness PreucilPercentHueError RGB2CMY rgb2dklCart rgb2dklV RGB2HSL RGB2HSV RGB2hue RGB2LEF RGB2LMS RGB2LSLM RGB2PhotoYCC RGB2XYZ RGB2YCbCr RGB2YIQ RGB2YPbPr RGB2YUV RxRyRz2XYZ saturationCIECAM02 saturationCIELAB saturationCIELABEvaLubbe saturationCIELUV spectra2CCT spectra2CRIGAIFSCI spectra2ISObrightness spectra2lux spectra2XYZ sprague StearnsStearnscorrection Stensby68.Whiteness Taube60.Whiteness WestlandBlacknessIndex wlnm2xyz wlnm2XYZ xFit_1931 xy2CCT.HernandezAndres xy2CCT.McCamy xyChromaticitiesVos1978 xyY2XYZ XYZ2BVR XYZ2CCT.Robertson XYZ2HunterLab XYZ2Lab XYZ2LMS XYZ2Luv XYZ2RGB XYZ2RxRyRz XYZ2xyY XYZ2Yuv XYZMoonSpencer1945 XYZTannenbaum1974 XYZtoRGB Y2MunsellV Y2MunsellVtable1D1535 YCbCr2RGB yFit_1931 YIQ2RGB YPbPr2RGB Yuv2Luv YUV2RGB Yuv2xy Yuv2XYZ Yxy2CIE1960UCS Yxy2Yuv zFit_1931

    
    
    
#   MunsellSpecToHVC()
#
#   convert Munsell notation to numeric Hue,Value,Chroma
#       Hue     is ASTM Hue in the interval (0,100]  (except Hue=0 for neutrals)
#       Value   is in the interval [0,10]
#       Chroma  is positive  (except Chroma=0 for neutrals)
#
#   arguments:
#       MunsellSpecString   a character vector of length n
#
#   return value: an nx3 matrix with the computed HVCs placed in the rows
#
#   author:  Glenn Davis
MunsellSpecToHVC <- function( MunsellSpecString )
    {
    if( ! is.character(MunsellSpecString) )
        {
        #   mess = "MunsellSpecToHVC().  MunsellSpecString is not a string.\n"
        #   cat( mess )
        return(NULL)
        }

    n = length(MunsellSpecString)
    
    out = matrix( as.numeric(NA), n, 3 )
    colnames(out)   = c( "Hue", "Value", "Chroma" )
    rownames(out)   = MunsellSpecString
        
    # Remove whitespace from input Munsell string
    MunsellSpecString <- gsub( ' |\t', '', MunsellSpecString )
    
    # Make all letters in Munsell string upper case
    MunsellSpecString <- toupper(MunsellSpecString)  

    #-----   chromatic colors   ----------#
    pattern = "^([0-9.]+)(R|YR|Y|GY|G|BG|B|PB|P|RP)([0-9.]+)/([0-9.]+)$" 

    #   time_start  = as.double( Sys.time() )
        
    sub1234 = sub( pattern, "\\1 \\2 \\3 \\4", MunsellSpecString )
    
    mask    = (nchar(MunsellSpecString) < nchar(sub1234) )     # a matched string gets longer          #mask = grepl( pattern, MunsellSpecString )   #   ; print( mask )
    
    if( any(mask) )
        {
        #   make a little LUT
        hue = seq( 0, 90, by=10 )
        names(hue)  = c("R","YR","Y","GY","G","BG","B","PB","P","RP")
        
        sub1234    = sub1234[mask]        
        
        dat     = unlist( strsplit( sub1234, ' ', fixed=T ) ) #;  print( dat )
        dat     = matrix( dat, length(sub1234), 4, byrow=T ) #;    print( dat )
        
        out[ mask, 1]   = hue[ dat[ ,2] ]  +  as.numeric(dat[ ,1])
        out[ mask, 2]   = as.numeric( dat[ ,3] )
        out[ mask, 3]   = as.numeric( dat[ ,4] )
        }
        
    #   cat( "Elapsed: ", as.double( Sys.time() ) - time_start, '\n' )
        
    #------   achromatic colors     ----#
    pattern = "^N([0-9.]+)(/|/0)?$"

    sub1    = sub( pattern, "\\1", MunsellSpecString )
        
    mask    = (nchar(sub1) < nchar(MunsellSpecString) ) # a matched string gets shorter.    #mask = grepl( pattern, MunsellSpecString )   #   ; print( mask )
 
    if( any(mask) )
        {
        out[ mask, 1]    = 0
        out[ mask, 2]    = as.numeric( sub1[mask] )
        out[ mask, 3]    = 0
        }
        
    return( out )
    }



#   ColorBlockFromMunsell()
#
#   HVC        a numeric 3-vector, or an nx3 matrix
#
#   if HVC[] is a 3-vector
#       HVC[1]     Munsell hue, on the ASTM D1535 100 point circular scale. All values are valid.
#       HVC[2]     Munsell value, must be between 0 and 10
#       HVC[3]     Munsell chroma, must be non-negative
#
#   if HVC[,] is an nx3 matrix
#       each row is defined as above
#
#   return value:
#       for a 3-vector:     list with: HVC, ISCC-NBS Number, ISCC-NBS Name
#       for an nx3 matrix:  data.frame with columns: HVC, ISCC-NBS Number, ISCC-NBS Name
#
#   the function requires these global data.frames:
#       get("SystemISCCNBS", envir = environment())
#       get("CentralsISCCNBS", envir = environment())
#
#   author:  Glenn Davis
ColorBlockFromMunsell  <-  function( HVC )
    {    
    if( ! is.numeric(HVC) )
        {
        #   mess = "ColorBlockFromMunsell().  HVC is not numeric.\n"
        #   cat( mess )
        return(NULL)
        }

    if( is.matrix(HVC) )
        {
        if( ncol(HVC) != 3 )
            {
            #   mess = "ColorBlockFromMunsell().  HVC does not have 3 columns.\n"
            #   cat( mess )
            return(NULL)
            }        
        
        colnames(HVC)  = c( "H", "V", "C" )        
        class(HVC)     = "model.matrix" 
        
        out = data.frame( HVC=HVC, Number=as.integer(NA), Name=as.character(NA), stringsAsFactors=FALSE )
        
        for( i in 1:nrow(out) )
            {
            result  = ColorBlockFromMunsell( HVC[i, ] )
        
            out$Number[i]   = result$Number
            out$Name[i]     = result$Name
            }
        
        return(out)
        }
    
    if( length(HVC) != 3 )
        {
        #   mess = "ColorBlockFromMunsell().  length of HVC is not 3.\n"
        #   cat( mess )
        return(NULL)
        }        
        

    if( is.null( names(HVC) ) )
        names(HVC) = c( "Hue", "Value", "Chroma" )
    
    out = list( HVC=HVC, Number=as.integer(NA), Name=as.character(NA) )
            
    vmax    = 10            
    valid   = all( is.finite(HVC) )  &&  (0 <= HVC[2])  &&  (HVC[2] <= vmax)  &&  (0 <= HVC[3])
    
    if( ! valid )   return( out )
    
    #   translate hue to the interval [1,101)   (101 is not included)
    HVC[1] = ((HVC[1] - 1 ) %% 100) + 1
        
    if( HVC[2] == vmax )   HVC[2] = vmax - 1.e-6     # because upper comparison below is strict
    
    #   do 6-way comparison.  
    #   Note upper comparisons are strict, and lower comparisons are not strict.
    #   So a point on a boundary is in only 1 block.
    mask.H  = get("SystemISCCNBS", envir = environment())$Hmin <= HVC[1]  &  HVC[1] < get("SystemISCCNBS", envir = environment())$Hmax
    mask.V  = get("SystemISCCNBS", envir = environment())$Vmin <= HVC[2]  &  HVC[2] < get("SystemISCCNBS", envir = environment())$Vmax
    mask.C  = get("SystemISCCNBS", envir = environment())$Cmin <= HVC[3]  &  HVC[3] < get("SystemISCCNBS", envir = environment())$Cmax
        
    theRow  = get("SystemISCCNBS", envir = environment())[ mask.H & mask.V & mask.C, ]
    
    if( nrow(theRow) != 1 )
        {
        #   mess = sprintf( "Expected exactly 1 match, but found %d\n", nrow(theRow) )
        #   cat( mess )
        return( out )
        }
        
    out$Number  = theRow$Number
    
    out$Name    = get("CentralsISCCNBS", envir = environment())$Name[ out$Number ] #
    
    return( out )
    }

#
#   CheckColorLookup()
#   performs lookup with colorBlockFromMunsell() and verifies that the color block number is correct
#
#   DataISCCNBS   data.frame with columns
#       MunsellSpec     a string in Munsell notation
#       Number          the correct ISCC-NBS Number
#
#   return value:  TRUE or FALSE
#
#   author:  Glenn Davis
CheckColorLookup <- function( DataISCCNBS=get("CentralsISCCNBS", envir = environment()) )#
    {
    hvc = MunsellSpecToHVC( DataISCCNBS$MunsellSpec )
    
    block   = ColorBlockFromMunsell( hvc )  #;      print( v )
    
    #   compute errors
    idx = which( block$Number != DataISCCNBS$Number )
    
    out = (length(idx) == 0)
    
    if( out )
        cat( "There are no errors !\n" )
    else
        {
        mess    = sprintf( "There are %d errors\n", length(idx) )
        cat(mess)
        df  = cbind( DataISCCNBS[ idx, ], NumberLookup=block$Number[idx] )
        print( df )
        }

    return( out )
    }

# David H. Brainard
# Cone Contrast and Opponent Modulation Color Spaces
# pp. 563
# PART IV: CONE CONTRAST AND
# OPPONENT MODULATION COLOR SPACES
# DKL Example
# 7/6/95 dhb Wrote it.
# September 2014 Jose Gama converted the code to R
LMS2DKL<-function(bg,diffcone.coords,DKL2LMS=FALSE){
# bg= the background vector for the conversion
# diffcone.coords = the vector we wish to convert
# STEP 1, 2: GET bg and diffcone.coords
# STEP 3: Set M.raw as in equation A.4.9.
# This is found by inserting the background
# values into equation A.4.8. Different
# backgrounds produce different matrices.
# The MATLAB notation below just
# fills the desired 3-by-3 matrix.
M.raw <- matrix(c(1,1,0,1,-bg[1]/bg[2],0,-1,-1,(bg[1]+bg[2])/bg[3]),3,3,byrow=TRUE)
# STEP 4: Compute the inverse of M for
# equation A.4.10. The MATLAB inv() function
# computes the matrix inverse of its argument.
M.raw.inv <- solve(M.raw)
# STEP 5: Find the three isolating stimuli as
# the columns of M.inv.raw. The MATLAB
# notation X[,i] extracts the i-th column
# of the matrix X.
isochrom.raw <- M.raw.inv[,1]
rgisolum.raw <- M.raw.inv[,2]
sisolum.raw <- M.raw.inv[,3]
# STEP 6: Find the pooled cone contrast of each
# of these. The MATLAB norm() function returns
# the vector length of its argument. The MATLAB
# / operation represents entry-by-entry division.
isochrom.raw.pooled <- norm(isochrom.raw / bg,'f')
rgisolum.raw.pooled <- norm(rgisolum.raw / bg,'f')
sisolum.raw.pooled <- norm(sisolum.raw / bg,'f')
# STEP 7: Scale each mechanism isolating
# modulation by its pooled contrast to obtain
# mechanism isolating modulations that have
# unit length.
isochrom.unit <-  matrix(isochrom.raw/ isochrom.raw.pooled)
rgisolum.unit <-  matrix(rgisolum.raw/ rgisolum.raw.pooled)
sisolum.unit <- matrix(sisolum.raw / sisolum.raw.pooled)
# STEP 8: Compute the values of the normalizing
# constants by plugging t~e unit isolating stimuli
# into A.4.9 and seeing what we get. Each vector
# should have only one non-zero entry. The size
# of the entry is the response of the unscaled
# mechanism to the stimulus that should give unit
# response.
lum.resp.raw <- M.raw %*% isochrom.unit
l.minus.m.resp.raw <- M.raw %*% rgisolum.unit
s.minus.lum.resp.raw <- M.raw %*% sisolum.unit
# STEP 9: We need to rescale the rows of M.raw
# so that we get unit response. This means
# mUltiplying each row of M.raw by a constant.
# The easiest way to accomplish the multiplication
# is to form a diagonal matrix with the desired
# scalars on the diagonal. These scalars are just
# the multiplicative inverses of the non-zero
# entries of the vectors obtained in the previous
# step. The resulting matrix M provides the
# entries of A.4.11. The three .resp vectors
# computed should be the three unit vectors
# (and they are).
D.rescale <- matrix(c(1/lum.resp.raw[1],0,0,0,1/l.minus.m.resp.raw[2],0,0,0,1/s.minus.lum.resp.raw[3]),3,3,byrow=TRUE)
M <- D.rescale %*% M.raw
lum.resp <- M %*% isochrom.unit
l.minus.m.resp <- M %*% rgisolum.unit
s.minus.lum.resp <- M %*% sisolum.unit
# STEP 10: Compute the inverse of M to obtain
# the matrix in equation A.4.12.
if (!DKL2LMS) {
M.inv <- solve(M)
M<-M.inv
}
# STEP 11: Multiply the vector we wish to
# convert byM to obtain its DKL coordinates.
DKL.coords <- M %*% diffcone.coords
# STEP 12: convert to spherical coordinates.
# According to the conventions in the original DKL
# paper, azimuth of 0 is along our rgisolum axis,
# azimuth of 90 is along our negative sisolum
# axis. The isochromatic axis has an elevation
# of 90 degrees. To do the conversion. we flip the
# sign of the sisolum coordinate and then do a
# standard conversion to polar coordinates .
RADS.TO.DEGS <- 360/(2*pi)
azimuth.rads <- atan(mrdivide(-DKL.coords[3], DKL.coords[2]))#pracma::
isolum.len <- sqrt(DKL.coords[2]^2 + DKL.coords[3]^2)
elevation.rads <- atan(mrdivide(DKL.coords[1], isolum.len))#pracma::
azimuth <- RADS.TO.DEGS %*% azimuth.rads
elevation <- RADS.TO.DEGS %*% elevation.rads
list(azimuth.rads=azimuth.rads,isolum.len=isolum.len,elevation.rads=elevation.rads,azimuth=azimuth,elevation=elevation)
}

XYZtoRGB<-function(xc,yc,zc, ColorSystem=c(0.67,0.33,0.21,0.71,0.14,0.08,0.310,0.316)){
# simplified XYZ to RGB conversion, mostly for plots
xr <- ColorSystem[1]
yr <- ColorSystem[2]
zr <- 1.0 - xr - yr
xg <- ColorSystem[3]
yg <- ColorSystem[4]
zg <- 1.0 - xg - yg
xb <- ColorSystem[5]
yb <- ColorSystem[6]
zb <- 1.0 - xb - yb
d <- xr*yg*zb - xg*yr*zb - xr*yb*zg + xb*yr*zg + xg*yb*zr - xb*yg*zr
R <- (-xg*yc*zb + xc*yg*zb + xg*yb*zc - xb*yg*zc - xc*yb*zg + xb*yc*zg) / d
G <- ( xr*yc*zb - xc*yr*zb - xr*yb*zc + xb*yr*zc + xc*yb*zr - xb*yc*zr) / d
B <- ( xr*yg*zc - xg*yr*zc - xr*yc*zg + xc*yr*zg + xg*yc*zr - xc*yg*zr) / d
cbind(R,G,B)
}

chromaticity.diagram<-function(chromaticityCoordinates=get("cccie31", envir = environment()), conversionFunction=NULL,...){#
# plot the chromaticity diagram AKA "horse shoe"
# conversionFunction CIE1931XYZ2CIE1976uv
pLen<-length(chromaticityCoordinates[["wlnm"]])
n=seq(1,pLen,by=1)
x<-chromaticityCoordinates$x[n]
y<-chromaticityCoordinates$y[n]
if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
z<-conversionFunction(cbind(chromaticityCoordinates$x[n],chromaticityCoordinates$y[n],chromaticityCoordinates$z[n]))
x<-z[,1]
y<-z[,2]
}
dots <-  list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
if (!lNameDot) dots <- modifyList(dots, list(x=x,y=y,type='l')) else dots <- modifyList(dots, list(x=x,y=y ))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='l'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
do.call(plot, dots ) # horseshoe
segments(x[1],y[1],x[pLen],y[pLen]) # line
}

Maxwell.triangle<-function(primariesRGB=get("whitepointsRGB", envir = environment()), conversionFunction=NULL,...){
# plot the Maxwell triangle
# conversionFunction CIE1931XYZ2CIE1976uv
x<-as.numeric(primariesRGB[1,c('xRed','xGreen','xBlue')])
y<-as.numeric(primariesRGB[1,c('yRed','yGreen','yBlue')])
if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
z<-1-x-y
z<-conversionFunction(cbind(x,y,z))
x<-z[,1]
y<-z[,2]
}
dots <- list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
if (!lNameDot) dots <- modifyList(dots, list(x=c(x,x[1]),y=c(y,y[1]),type='l')) else dots <- modifyList(dots, list(x=c(x,x[1]),y=c(y,y[1])))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='l'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
do.call(plot, dots )
}

chromaticity.diagram.color<-function(chromaticityCoordinates=get("cccie31", envir = environment()), conversionFunction=NULL,granularity=10,...){
# plot the chromaticity diagram AKA "horse shoe"
# conversionFunction CIE1931XYZ2CIE1976uv
pLen<-length(chromaticityCoordinates[["wlnm"]])
n=seq(1,pLen,by=1)
wlnmVector<-seq(chromaticityCoordinates[1,"wlnm"],chromaticityCoordinates[pLen,"wlnm",],by=1/granularity)
pLen2<-length(wlnmVector)
XYZ<-matrix(unlist(wlnm2xyz(wlnmVector)),ncol=3,byrow=FALSE)
x<-XYZ[,1]
y<-XYZ[,2]
z<-XYZ[,3]
X<-x
Y<-y
Z<-z
if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
az<-conversionFunction(cbind(x,y,z))
x<-az[,1]
y<-az[,2]
z<-1-x-y
}
# temp<-XYZtoRGB(X,Y,Z)/255
# t1<-temp/max(temp)*255
temp1<-xyz2srgb(cbind(X,Y,Z))
temp<-temp1[["sRGB"]]/255
t1<-temp/max(temp)*255
t1<-round(t1)
t1[which(t1<0)]<-0
t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
dots <- list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
if (!lNameDot) dots <- modifyList(dots, list(x=x,y=y,type='p')) else dots <- modifyList(dots, list(x=x,y=y))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
if (!(any(nameDots=='col'))) dots <- modifyList(dots, list(col=t2))
if (!(any(nameDots=='pch'))) dots <- modifyList(dots, list(pch=20))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='p'))
do.call(plot, dots );par(new=TRUE)
#line
x1<-seq(x[1],x[pLen2],length.out=pLen2/7)#
y1<-seq(y[1],y[pLen2],length.out=pLen2/7)#
X1<-seq(X[1],X[pLen2],length.out=pLen2/7)#
Y1<-seq(Y[1],Y[pLen2],length.out=pLen2/7)#
# temp<-XYZtoRGB(X1,Y1,1-X1-Y1)/255
# t1<-temp/max(temp)*255
temp1<-xyz2srgb(cbind(X1,Y1,1-X1-Y1))
temp<-temp1[["sRGB"]]/255
t1<-temp/max(temp)*255
t1<-round(t1)
t1[which(t1<0)]<-0
t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
dots <- modifyList(dots, list(x=x1))
dots <- modifyList(dots, list(y=y1))
dots <- modifyList(dots, list(col=t2))
do.call(plot, dots )
}

Maxwell.triangle.color<-function(primariesRGB=get("whitepointsRGB", envir = environment()), conversionFunction=NULL,granularity=10,...){
# Maxwell triangle
# conversionFunction CIE1931XYZ2CIE1976uv
pLen<-100
n=seq(1,pLen,by=1)
x<-as.numeric(primariesRGB[1,c('xRed','xGreen','xBlue')])
y<-as.numeric(primariesRGB[1,c('yRed','yGreen','yBlue')])
x1<-seq(x[1],x[2],length.out=pLen)
y1<-seq(y[1],y[2],length.out=pLen)
x2<-seq(x[1],x[3],length.out=pLen)
y2<-seq(y[1],y[3],length.out=pLen)
x3<-seq(x[2],x[3],length.out=pLen)
y3<-seq(y[2],y[3],length.out=pLen)
x<-c(x1,x2,x3)
y<-c(y1,y2,y3)
z<-1-x-y
X<-x
Y<-y
Z<-z
if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
z<-conversionFunction(cbind(x,y,z))
x<-z[,1]
y<-z[,2]
}
# temp<-XYZtoRGB(X,Y,Z)/255
# t1<-temp/max(temp)*255
temp1<-xyz2srgb(cbind(X,Y,Z))
temp<-temp1[["sRGB"]]/255
t1<-temp/max(temp)*255
t1<-round(t1)
t1[which(t1<0)]<-0
t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
dots <- list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
if (!lNameDot) dots <- modifyList(dots, list(x=x,y=y,type='p')) else dots <- modifyList(dots, list(x=x,y=y))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
if (!(any(nameDots=='col'))) dots <- modifyList(dots, list(col=t2))
if (!(any(nameDots=='pch'))) dots <- modifyList(dots, list(pch=20))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='p'))
do.call(plot, dots)
}

Maxwell.triangle.color.fill<-function(chromaticityCoordinates=get("cccie31", envir = environment()), conversionFunction=NULL,granularity=10,...){
# plot the Maxwell triangle
# conversionFunction CIE1931XYZ2CIE1976uv
pLen<-length(chromaticityCoordinates[["wlnm"]])
n=seq(1,pLen,by=1)
wlnmVector<-seq(chromaticityCoordinates[1,"wlnm"],chromaticityCoordinates[pLen,"wlnm",],by=1/granularity)
pLen2<-length(wlnmVector)
XYZ<-matrix(unlist(wlnm2xyz(wlnmVector)),ncol=3,byrow=FALSE)
x<-XYZ[,1]
y<-XYZ[,2]
z<-XYZ[,3]
X<-x
Y<-y
Z<-z
if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
az<-conversionFunction(cbind(x,y,z))
x<-az[,1]
y<-az[,2]
}
maxx<-max(x)
maxy<-max(y)
minx<-min(x)
miny<-min(y)
pLen2<-10*granularity;k=1
mx<-seq(minx,maxx, length.out=pLen2)
my<-seq(miny,maxy, length.out=pLen2)
p2<-matrix(0,pLen2^2,2)
p3<-vector(mode='character',pLen2^2)
for (n in 1:pLen2)# triangle
for (m in 1:pLen2){
   XC=mx[m]
   YC=my[n]
   ZC=1-XC-YC
   if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
az<-conversionFunction(cbind(XC,YC,ZC))
XC<-az[,1]
YC<-az[,2]
ZC=1-XC-YC
}
# temp<-XYZtoRGB(mx[m],my[n],1-mx[m]-my[n])/255
# t1<-temp/max(temp)*255
temp1<-xyz2srgb(cbind(mx[m],my[n],1-mx[m]-my[n]))
temp<-temp1[["sRGB"]]/255
t1<-temp/max(temp)*255
t1<-round(t1)
if (any(t1<0)) t2<-NA else t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
   p2[k,1]<- XC
   p2[k,2]<- YC
   p3[k]<- t2
   k<-k+1
}
dots <- list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
if (!lNameDot) dots <- modifyList(dots, list(x=p2[,1],y=p2[,2],col=p3,type='p')) else dots <- modifyList(dots, list(x=p2[,1],y=p2[,2],col=p3))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
if (!(any(nameDots=='col'))) dots <- modifyList(dots, list(col=p3))
if (!(any(nameDots=='pch'))) dots <- modifyList(dots, list(pch=15))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='p'))
do.call(plot, dots );par(new=TRUE)
}


chromaticity.diagram.color.fill<-function(chromaticityCoordinates=get("cccie31", envir = environment()), conversionFunction=NULL,granularity=10, conversionFunctionInv=NULL,...){
# plot the chromaticity diagram AKA "horse shoe"
# conversionFunction CIE1931XYZ2CIE1976uv
pLen<-length(chromaticityCoordinates[["wlnm"]])
#n=seq(1,pLen,by=1)
wlnmVector<-seq(chromaticityCoordinates[1,"wlnm"],chromaticityCoordinates[pLen,"wlnm",],by=1/granularity)
pLen2<-length(wlnmVector)
XYZ<-matrix(unlist(wlnm2xyz(wlnmVector)),ncol=3,byrow=FALSE)
x<-XYZ[,1]
y<-XYZ[,2]
z<-XYZ[,3]
X<-x
Y<-y
Z<-z
if (is.function(conversionFunction)) {
if (!is.function(conversionFunctionInv)){
mystrfun<-as.character(substitute(conversionFunction))
if (mystrfun=='CIE1931xy2CIE1976uv') conversionFunctionInv<-CIE1976uv2CIE1931xy
if (mystrfun=='CIE1931xy2CIE1960uv') conversionFunctionInv<-CIE1960UCS2xy # CIE1960uv2CIE1931xy
if (mystrfun=='CIE1976uv2CIE1931xy') conversionFunctionInv<-CIE1931xy2CIE1976uv
if (mystrfun=='CIE1960uv2CIE1931xy') conversionFunctionInv<-CIE1931xy2CIE1960uv
}
az<-conversionFunction(cbind(x,y,z))
x<-az[,1]
y<-az[,2]
}
boundariesXY<-cbind(x,y)
#line
#x1<-seq(x[1],x[length(x)],length.out=pLen2)
#y1<-seq(y[1],y[length(y)],length.out=pLen2)
#limits
maxx<-max(X)
maxy<-max(Y)
minx<-min(X)
miny<-min(Y)
pLen2<-10*granularity;k=1
dots <- list(...)
nameDots<-names(dots)
lNameDot<-length(dots)>0
boolDots<-FALSE
if ((any(nameDots=='xlim'))) { minx<-dots$xlim[1]
maxx<-dots$xlim[2]
boolDots<-TRUE
}
if ((any(nameDots=='ylim'))) { miny<-dots$ylim[1]
maxy<-dots$ylim[2]
boolDots<-TRUE
}
containedInTheHorseShoe<-FALSE
if (boolDots){
if (is.function(conversionFunctionInv)){
minx<-conversionFunctionInv(c(minx, miny))[1]
maxx<-conversionFunctionInv(c(maxx, maxy))[1]
miny<-conversionFunctionInv(c(minx, miny))[2]
maxy<-conversionFunctionInv(c(maxx, maxy))[2]
}
# are xlim, ylim contained in the horse shoe?
if (all(point.in.polygon(c(minx, maxx, maxx, minx),c(miny, maxy, miny, maxy),boundariesXY[,1],boundariesXY[,2]))) containedInTheHorseShoe<-TRUE
else { # 
maxx<-max(X)
maxy<-max(Y)
minx<-min(X)
miny<-min(Y)
}
}
mx<-seq(minx,maxx, length.out=pLen2)
pLen3<-round(pLen2 /( (maxx-minx)/(maxy-miny) ))
my<-seq(miny,maxy, length.out=pLen3)
# p2<-matrix(0,pLen2*pLen3,2)
# p3<-vector(mode='character',pLen2*pLen3)
# for (n in 1:pLen3)# Horse Shoe
# for (m in 1:pLen2){
#    XC=mx[m]
#    YC=my[n]
#    ZC=1-XC-YC
#    if (!is.null(conversionFunction)) if (is.function(conversionFunction)) {
# az<-conversionFunction(cbind(XC,YC,ZC))
# XC<-az[,1]
# YC<-az[,2]
# ZC=1-XC-YC
# }
# temp<-XYZtoRGB(mx[m],my[n],1-mx[m]-my[n])/255
# t1<-temp/max(temp)*255
# t1<-round(t1)
# 
# #boundariesXY<-rbind(boundariesXY,cbind(x1,y1))
# #boundariesXY<-rbind(boundariesXY,boundariesXY[1,])
# if (any(t1<0)) t1[which(t1<0)]<-1
# #if (!in.out(boundariesXY,cbind(XC,YC))) t2<-NA else t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
# #if (!pnt.in.poly(cbind(XC,YC),boundariesXY)) t2<-NA else t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
# 
# if (containedInTheHorseShoe) {
# t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
# } else {
# if (!point.in.polygon(XC,YC,boundariesXY[,1],boundariesXY[,2])) t2<-NA else t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
# }
# 
#    p2[k,1]<- XC
#    p2[k,2]<- YC
#    p3[k]<- t2
#    k<-k+1
# }
#vectorized version
n  = rep(1:pLen3,each=pLen2)
m  = rep(1:pLen2,pLen3)
   XC=mx[m]
   YC=my[n]
   ZC=1-XC-YC
if (is.function(conversionFunction)) {
az<-conversionFunction(cbind(XC,YC,ZC))
XC<-az[,1]
YC<-az[,2]
ZC=1-XC-YC
}
# temp<-XYZtoRGB(mx[m],my[n],1-mx[m]-my[n])/255
# t1<-temp/apply(temp,1,max)*255
temp1<-xyz2srgb(cbind(mx[m],my[n],1-mx[m]-my[n]))

temp<-temp1[["sRGB"]]/255
t1<-temp/apply(temp,1,max)*255
t1<-round(t1)

# t1<-temp1[["sRGB"]]

if (any(t1<0)) t1[which(t1<0)]<-1
if (containedInTheHorseShoe) {
t2<-apply(t1,1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
} else {
boolInPolygon<-point.in.polygon(XC,YC,boundariesXY[,1],boundariesXY[,2])
t2<-rep(NA,pLen2*pLen3)
t2[which(boolInPolygon==1)]<-apply(t1[which(boolInPolygon==1),],1,function(x) sprintf("#%02X%02X%02X",x[1],x[2],x[3]))
}
p2<-cbind(XC,YC)
p3<-t2
if (!lNameDot) dots <- modifyList(dots, list(x=p2[,1],y=p2[,2],col=p3,type='p')) else dots <- modifyList(dots, list(x=p2[,1],y=p2[,2],col=p3))
if (!(any(nameDots=='xlab'))) dots <- modifyList(dots, list(xlab='x'))
if (!(any(nameDots=='ylab'))) dots <- modifyList(dots, list(ylab='y'))
if (!(any(nameDots=='xlim'))) dots <- modifyList(dots, list(xlim=c(0,0.85)))
if (!(any(nameDots=='ylim'))) dots <- modifyList(dots, list(ylim=c(0,0.85)))
#if (!(any(nameDots=='col'))) dots <- modifyList(dots, list(col=p3))
if (!(any(nameDots=='pch'))) dots <- modifyList(dots, list(pch=15))
if (!(any(nameDots=='type'))) dots <- modifyList(dots, list(type='p'))
do.call(plot, dots );par(new=TRUE)
}


LUV2LAB<-function(Luvmatrix) XYZ2Lab(Luv2XYZ(Luvmatrix))
LAB2LUV<-function(Labmatrix) XYZ2Luv(Lab2XYZ(Labmatrix))
XYZ2xyY <- function(XYZmatrix) {
rSum<-apply(XYZmatrix,1,sum)
cbind(XYZmatrix[,1]/ rSum ,XYZmatrix[,2]/ rSum,XYZmatrix[,3])
}

makeChromaticAdaptationMatrix<-function(ChromaticAdaptationAlgorithm='VonKries', illuminantSource='C', illuminantDestination='D65', 
observer=2, ChromaticAdaptationArray=get("ChromaticAdaptation", envir = environment()), referenceWhiteArray=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Chromatic Adaptation
which(referenceWhiteArray[,"Illuminant"]==illuminantDestination)
if (observer==2) observerPos<-2:4 else observerPos<-5:7
coneResponseDomainSource<-ChromaticAdaptationArray[,,ChromaticAdaptationAlgorithm,'direct'] %*% 
matrix(unlist(referenceWhiteArray[which(referenceWhiteArray[,"Illuminant"]==illuminantSource),observerPos]),3,1,byrow=TRUE)
coneResponseDomainDestination<-ChromaticAdaptationArray[,,ChromaticAdaptationAlgorithm,'direct'] %*% 
matrix(unlist(referenceWhiteArray[which(referenceWhiteArray[,"Illuminant"]==illuminantDestination),observerPos]),3,1,byrow=TRUE)
d <- diag(3)
diag(d) <- coneResponseDomainDestination/coneResponseDomainSource
ChromaticAdaptationArray[,,ChromaticAdaptationAlgorithm,'inverse'] %*%  d   %*%  ChromaticAdaptationArray[,,ChromaticAdaptationAlgorithm,'direct'] #  adaptation matrix
}

footcandle2lux<-function(ftcl) ftcl * 10.7639104167
# converts foot candle to Lumens/lux
# http://www.translatorscafe.com/cafe/EN/units-converter/illumination
footcandle2watt.sqcentimeter<-function(ftcl) ftcl*0.00000157597517082
# converts foot candle to watts / square centimeter [w/cm^2] (at 555 nm) 
# http://www.translatorscafe.com/cafe/EN/units-converter/illumination
footcandle2candela.steradian.sqmeter<-function(ftcl) ftcl * 10763.9104167 
# converts foot candle to candela steradian / square meter [cd*sr/m^2]
# http://www.translatorscafe.com/cafe/EN/units-converter/illumination

dkl2rgb <- function(dklMatrix, conversionMatrix=NA)
{ # Convert from DKL color space (Derrington, Krauskopf & Lennie) to RGB.
# Package psychopy for Python
# http://www.psychopy.org/epydoc/psychopy.misc-pysrc.html#dkl2rgb
if (!is.matrix(dklMatrix)) dklMatrix <- matrix(dklMatrix, ncol=3,byrow=TRUE)
elev = dklMatrix[,1]
azim = dklMatrix[,2]
radius = dklMatrix[,3]
RGBYLUM <- cbind(radius * sin(pi/180*(elev)), radius * cos(pi/180*(elev))*cos(pi/180*(azim)), radius * cos(pi/180*(elev))*sin(pi/180*(azim)) )

if (is.na(conversionMatrix)) conversionMatrix <- matrix(c(
  #LUMIN      %L-M    %L+M-S  (note that dkl has to be in cartesian coords first!) 
  1.0000, 1.0000, -0.1462,#R 
  1.0000, -0.3900, 0.2094,#G 
  1.0000, 0.0180, -1.0000#B 
),3,3,byrow=TRUE)

rgb <- conversionMatrix %*% t(RGBYLUM)
return (t(rgb))#return in the shape we received it
}

dklCart2rgb <- function(dklMatrix, conversionMatrix=NA)
{ # Like dkl2rgb except that it uses cartesian coords (LM,S,LUM) rather than
# spherical coords for DKL (elev, azim, contr)
if (!is.matrix(dklMatrix)) dklMatrix <- matrix(dklMatrix, ncol=3,byrow=TRUE)
if (is.na(conversionMatrix)) conversionMatrix <- matrix(c(
  #LUMIN      %L-M    %L+M-S  (note that dkl has to be in cartesian coords first!) 
  1.0000, 1.0000, -0.1462,#R 
  1.0000, -0.3900, 0.2094,#G 
  1.0000, 0.0180, -1.0000#B 
),3,3,byrow=TRUE)

rgb <- conversionMatrix %*% t(dklMatrix)
return (t(rgb))#return in the shape we received it
}

rgb2dklCart <- function(rgbMatrix, conversionMatrix=NA)
{ # Converts RGB into Cartesian DKL space
if (!is.matrix(rgbMatrix)) rgbMatrix <- matrix(rgbMatrix, ncol=3,byrow=TRUE)
if (is.na(conversionMatrix)) conversionMatrix <- matrix(c(
0.25145542,  0.64933633,  0.09920825,
0.78737943, -0.55586618, -0.23151325,
0.26562825,  0.63933074, -0.90495899
),3,3,byrow=TRUE)

dkl <- conversionMatrix %*% t(rgbMatrix)
return (t(dkl))#return in the shape we received it
}



rgb2dklV <- function(RGB){
# convert RGB to DKL
# Graph-Based Visual Saliency (MATLAB source code)
# Jonathan Harel
# California Institute of Technology
# http://www.klab.caltech.edu/~harel/share/gbvs.php
# constants used for RGB -> DKL conversion:
lutRGB <- matrix(c(
    0.024935, 0.0076954, 0.042291,
    0.024974, 0.0077395, 0.042346,
    0.025013, 0.0077836, 0.042401,
    0.025052, 0.0078277, 0.042456,
    0.025091, 0.0078717, 0.042511,
    0.02513, 0.0079158, 0.042566,
    0.025234, 0.007992, 0.042621,
    0.025338, 0.0080681, 0.042676,
    0.025442, 0.0081443, 0.042731,
    0.025545, 0.0082204, 0.042786,
    0.025649, 0.0082966, 0.042841,
    0.025747, 0.0084168, 0.042952,
    0.025844, 0.0085371, 0.043062,
    0.025942, 0.0086573, 0.043172,
    0.026039, 0.0087776, 0.043282,
    0.026136, 0.0088978, 0.043392,
    0.026234, 0.0090581, 0.043502,
    0.026331, 0.0092184, 0.043612,
    0.026429, 0.0093788, 0.043722,
    0.026526, 0.0095391, 0.043833,
    0.026623, 0.0096994, 0.043943,
    0.026818, 0.0099198, 0.044141,
    0.027013, 0.01014, 0.044339,
    0.027208, 0.010361, 0.044537,
    0.027403, 0.010581, 0.044736,
    0.027597, 0.010802, 0.044934,
    0.027857, 0.010994, 0.04522,
    0.028117, 0.011186, 0.045507,
    0.028377, 0.011379, 0.045793,
    0.028636, 0.011571, 0.046079,
    0.028896, 0.011764, 0.046366,
    0.029104, 0.012068, 0.046652,
    0.029312, 0.012373, 0.046938,
    0.029519, 0.012677, 0.047225,
    0.029727, 0.012982, 0.047511,
    0.029935, 0.013287, 0.047797,
    0.030273, 0.013663, 0.048326,
    0.03061, 0.01404, 0.048855,
    0.030948, 0.014417, 0.049383,
    0.031286, 0.014794, 0.049912,
    0.031623, 0.01517, 0.050441,
    0.032156, 0.015707, 0.051035,
    0.032688, 0.016244, 0.05163,
    0.033221, 0.016782, 0.052225,
    0.033753, 0.017319, 0.052819,
    0.034286, 0.017856, 0.053414,
    0.034961, 0.018693, 0.054493,
    0.035636, 0.019531, 0.055573,
    0.036312, 0.020369, 0.056652,
    0.036987, 0.021206, 0.057731,
    0.037662, 0.022044, 0.058811,
    0.038623, 0.023246, 0.060044,
    0.039584, 0.024449, 0.061278,
    0.040545, 0.025651, 0.062511,
    0.041506, 0.026854, 0.063744,
    0.042468, 0.028056, 0.064978,
    0.043857, 0.029659, 0.066806,
    0.045247, 0.031263, 0.068634,
    0.046636, 0.032866, 0.070463,
    0.048026, 0.034469, 0.072291,
    0.049416, 0.036072, 0.074119,
    0.051221, 0.038156, 0.076476,
    0.053026, 0.04024, 0.078833,
    0.054831, 0.042325, 0.081189,
    0.056636, 0.044409, 0.083546,
    0.058442, 0.046493, 0.085903,
    0.06039, 0.048737, 0.087996,
    0.062338, 0.050982, 0.090088,
    0.064286, 0.053226, 0.092181,
    0.066234, 0.055471, 0.094273,
    0.068182, 0.057715, 0.096366,
    0.070519, 0.06012, 0.098921,
    0.072857, 0.062525, 0.10148,
    0.075195, 0.06493, 0.10403,
    0.077532, 0.067335, 0.10659,
    0.07987, 0.069739, 0.10914,
    0.082208, 0.072345, 0.11176,
    0.084545, 0.07495, 0.11438,
    0.086883, 0.077555, 0.117,
    0.089221, 0.08016, 0.11963,
    0.091558, 0.082766, 0.12225,
    0.094026, 0.085611, 0.12533,
    0.096494, 0.088457, 0.12841,
    0.098961, 0.091303, 0.1315,
    0.10143, 0.094148, 0.13458,
    0.1039, 0.096994, 0.13767,
    0.10688, 0.10028, 0.14119,
    0.10987, 0.10357, 0.14471,
    0.11286, 0.10685, 0.14824,
    0.11584, 0.11014, 0.15176,
    0.11883, 0.11343, 0.15529,
    0.12208, 0.11695, 0.15903,
    0.12532, 0.12048, 0.16278,
    0.12857, 0.12401, 0.16652,
    0.13182, 0.12754, 0.17026,
    0.13506, 0.13106, 0.17401,
    0.1387, 0.13499, 0.17819,
    0.14234, 0.13892, 0.18238,
    0.14597, 0.14285, 0.18656,
    0.14961, 0.14677, 0.19075,
    0.15325, 0.1507, 0.19493,
    0.15727, 0.15519, 0.19956,
    0.1613, 0.15968, 0.20419,
    0.16532, 0.16417, 0.20881,
    0.16935, 0.16866, 0.21344,
    0.17338, 0.17315, 0.21806,
    0.17805, 0.17796, 0.22291,
    0.18273, 0.18277, 0.22775,
    0.1874, 0.18758, 0.2326,
    0.19208, 0.19238, 0.23744,
    0.19675, 0.19719, 0.24229,
    0.20156, 0.20224, 0.24758,
    0.20636, 0.20729, 0.25286,
    0.21117, 0.21234, 0.25815,
    0.21597, 0.21739, 0.26344,
    0.22078, 0.22244, 0.26872,
    0.2261, 0.22806, 0.27423,
    0.23143, 0.23367, 0.27974,
    0.23675, 0.23928, 0.28524,
    0.24208, 0.24489, 0.29075,
    0.2474, 0.2505, 0.29626,
    0.25299, 0.25651, 0.3022,
    0.25857, 0.26253, 0.30815,
    0.26416, 0.26854, 0.3141,
    0.26974, 0.27455, 0.32004,
    0.27532, 0.28056, 0.32599,
    0.28156, 0.28697, 0.33238,
    0.28779, 0.29339, 0.33877,
    0.29403, 0.2998, 0.34515,
    0.30026, 0.30621, 0.35154,
    0.30649, 0.31263, 0.35793,
    0.3126, 0.31904, 0.36388,
    0.3187, 0.32545, 0.36982,
    0.32481, 0.33186, 0.37577,
    0.33091, 0.33828, 0.38172,
    0.33701, 0.34469, 0.38767,
    0.34325, 0.3511, 0.39361,
    0.34948, 0.35752, 0.39956,
    0.35571, 0.36393, 0.40551,
    0.36195, 0.37034, 0.41145,
    0.36818, 0.37675, 0.4174,
    0.37429, 0.38317, 0.42313,
    0.38039, 0.38958, 0.42885,
    0.38649, 0.39599, 0.43458,
    0.3926, 0.4024, 0.44031,
    0.3987, 0.40882, 0.44604,
    0.40494, 0.41523, 0.45198,
    0.41117, 0.42164, 0.45793,
    0.4174, 0.42806, 0.46388,
    0.42364, 0.43447, 0.46982,
    0.42987, 0.44088, 0.47577,
    0.43623, 0.44689, 0.48128,
    0.4426, 0.45291, 0.48678,
    0.44896, 0.45892, 0.49229,
    0.45532, 0.46493, 0.4978,
    0.46169, 0.47094, 0.5033,
    0.46792, 0.47695, 0.50837,
    0.47416, 0.48297, 0.51344,
    0.48039, 0.48898, 0.5185,
    0.48662, 0.49499, 0.52357,
    0.49286, 0.501, 0.52863,
    0.49805, 0.50701, 0.53392,
    0.50325, 0.51303, 0.53921,
    0.50844, 0.51904, 0.54449,
    0.51364, 0.52505, 0.54978,
    0.51883, 0.53106, 0.55507,
    0.52442, 0.53667, 0.55969,
    0.53, 0.54228, 0.56432,
    0.53558, 0.5479, 0.56894,
    0.54117, 0.55351, 0.57357,
    0.54675, 0.55912, 0.57819,
    0.55182, 0.56433, 0.58304,
    0.55688, 0.56954, 0.58789,
    0.56195, 0.57475, 0.59273,
    0.56701, 0.57996, 0.59758,
    0.57208, 0.58517, 0.60242,
    0.57675, 0.58998, 0.60639,
    0.58143, 0.59479, 0.61035,
    0.5861, 0.5996, 0.61432,
    0.59078, 0.60441, 0.61828,
    0.59545, 0.60922, 0.62225,
    0.60065, 0.61403, 0.62709,
    0.60584, 0.61884, 0.63194,
    0.61104, 0.62365, 0.63678,
    0.61623, 0.62846, 0.64163,
    0.62143, 0.63327, 0.64648,
    0.62584, 0.63808, 0.65088,
    0.63026, 0.64289, 0.65529,
    0.63468, 0.6477, 0.65969,
    0.63909, 0.65251, 0.6641,
    0.64351, 0.65731, 0.6685,
    0.64857, 0.66132, 0.67269,
    0.65364, 0.66533, 0.67687,
    0.6587, 0.66934, 0.68106,
    0.66377, 0.67335, 0.68524,
    0.66883, 0.67735, 0.68943,
    0.67273, 0.68136, 0.69361,
    0.67662, 0.68537, 0.6978,
    0.68052, 0.68938, 0.70198,
    0.68442, 0.69339, 0.70617,
    0.68831, 0.69739, 0.71035,
    0.69221, 0.7022, 0.7141,
    0.6961, 0.70701, 0.71784,
    0.7, 0.71182, 0.72159,
    0.7039, 0.71663, 0.72533,
    0.70779, 0.72144, 0.72907,
    0.71169, 0.72505, 0.73348,
    0.71558, 0.72866, 0.73789,
    0.71948, 0.73226, 0.74229,
    0.72338, 0.73587, 0.7467,
    0.72727, 0.73948, 0.7511,
    0.73247, 0.74349, 0.75507,
    0.73766, 0.74749, 0.75903,
    0.74286, 0.7515, 0.763,
    0.74805, 0.75551, 0.76696,
    0.75325, 0.75952, 0.77093,
    0.75714, 0.76393, 0.77599,
    0.76104, 0.76834, 0.78106,
    0.76494, 0.77275, 0.78612,
    0.76883, 0.77715, 0.79119,
    0.77273, 0.78156, 0.79626,
    0.77792, 0.78677, 0.80132,
    0.78312, 0.79198, 0.80639,
    0.78831, 0.79719, 0.81145,
    0.79351, 0.8024, 0.81652,
    0.7987, 0.80762, 0.82159,
    0.80519, 0.81283, 0.82687,
    0.81169, 0.81804, 0.83216,
    0.81818, 0.82325, 0.83744,
    0.82468, 0.82846, 0.84273,
    0.83117, 0.83367, 0.84802,
    0.83636, 0.83888, 0.85286,
    0.84156, 0.84409, 0.85771,
    0.84675, 0.8493, 0.86256,
    0.85195, 0.85451, 0.8674,
    0.85714, 0.85972, 0.87225,
    0.86364, 0.86613, 0.87819,
    0.87013, 0.87255, 0.88414,
    0.87662, 0.87896, 0.89009,
    0.88312, 0.88537, 0.89604,
    0.88961, 0.89178, 0.90198,
    0.8961, 0.8986, 0.90947,
    0.9026, 0.90541, 0.91696,
    0.90909, 0.91222, 0.92445,
    0.91558, 0.91904, 0.93194,
    0.92208, 0.92585, 0.93943,
    0.92857, 0.93307, 0.94493,
    0.93506, 0.94028, 0.95044,
    0.94156, 0.94749, 0.95595,
    0.94805, 0.95471, 0.96145,
    0.95455, 0.96192, 0.96696,
    0.96364, 0.96954, 0.97357,
    0.97273, 0.97715, 0.98018,
    0.98182, 0.98477, 0.98678,
    0.99091, 0.99238, 0.99339,
    1, 1, 1 ), ncol=3,byrow=TRUE)

lms0 <- c( 34.918538957799996, 19.314796676499999, 0.585610818500000 )

m <- c( 18.32535,  44.60077,   7.46216, 4.09544,  28.20135,   6.66066, 0.02114,   0.10325,   1.05258 )

fac <- 1.0 / (lms0[1] + lms0[2])
mm <- c( sqrt(3.0)*fac, sqrt(3.0)*fac, 0.0, sqrt(lms0[1]*lms0[1]+lms0[2]*lms0[2])/lms0[1]*fac, 
-sqrt(lms0[1]*lms0[1]+lms0[2]*lms0[2])/lms0[2]*fac, 0.0, -fac, -fac, (lms0[1] + lms0[2]) / lms0[3] * fac )

# do a lookup RGB -> rgb:
if ( max(RGB)<=1 ){
  RGB <- ceiling( RGB*255 )
  RGB[which(RGB<1)] <- 1
  RGB[which(RGB>256)] <- 256
} else   RGB <- RGB + 1

aa1 <- lutRGB[ RGB[,1] , 1 ]
aa2 <- lutRGB[ RGB[,2] , 2 ]
aa3 <- lutRGB[ RGB[,3] , 3 ]

# now convert to LMS:
lms1 <- m[1] * aa1 + m[2] * aa2 + m[3] * aa3 - lms0[1]
lms2 <- m[4] * aa1 + m[5] * aa2 + m[6] * aa3 - lms0[2]
lms3 <- m[7] * aa1 + m[8] * aa2 + m[9] * aa3 - lms0[3]

# finally to DKL:
dkl1 <- mm[1] * lms1 + mm[2] * lms2 + mm[3] * lms3
dkl2 <- mm[4] * lms1 + mm[5] * lms2 + mm[6] * lms3
dkl3 <- mm[7] * lms1 + mm[8] * lms2 + mm[9] * lms3

# finally to DKLn:
dkl <- c( dkl1 * 0.5774, dkl2 * 2.7525, dkl3 * 0.4526 )
dkl
}

dkl2dklCart<- function(dklMatrix)
{ # Convert from spherical coords to cartesian coords
if (!is.matrix(dklMatrix)) dklMatrix <- matrix(dklMatrix, ncol=3,byrow=TRUE)
elev = dklMatrix[,1]
azim = dklMatrix[,2]
radius = dklMatrix[,3]
cbind(radius * cos(pi/180*(elev))*sin(pi/180*(azim)), radius * cos(pi/180*(elev))*cos(pi/180*(azim)), radius * sin(pi/180*(elev)) )
}

dklCart2dkl <- function(dklMatrix)
{ # Convert from cartesian coords to spherical coords
if (!is.matrix(dklMatrix)) dklMatrix <- matrix(dklMatrix, ncol=3,byrow=TRUE)
n = dklMatrix[,1]
m = dklMatrix[,2]
p = dklMatrix[,3]
radius <- sqrt(n^2+m^2+p^2)
elevation <- asin(p/radius)
azimuth <- atan2(n,m)
elevation <- elevation *180/pi
azimuth <- azimuth *180/pi
cbind(elevation,azimuth,radius)
}

#CIE 1976 Luv to u', v' CIE 1976
Luv2Yuv<-function(Luvmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment())) 
{
if (!is.matrix(Luvmatrix)) Luvmatrix <- matrix(Luvmatrix, ncol=3,byrow=TRUE)
XYZ2Yuv(Luv2XYZ(Luvmatrix,illuminant,observer,RefWhite))
}
#CIE u', v' CIE 1976 to CIE 1976 Luv
Yuv2Luv<-function(Yu.v.matrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment())) 
{
if (!is.matrix(Yu.v.matrix)) Yu.v.matrix <- matrix(Yu.v.matrix, ncol=3,byrow=TRUE)
XYZ2Luv(Yuv2XYZ(Yu.v.matrix),illuminant,observer,RefWhite)
}
# chromaticity coordinates u', v' CIE 1976
# source: CIE 1976 color space, From Wikipedia, the free encyclopedia http://en.wikipedia.org/wiki/CIE_1976_color_space
XYZ2Yuv<-function(XYZmatrix) {
if (!is.matrix(XYZmatrix)) XYZmatrix <- matrix(XYZmatrix, ncol=3,byrow=TRUE)
cbind(Y=XYZmatrix[,2],u76=4*XYZmatrix[,1]/(XYZmatrix[,1]+15*XYZmatrix[,2]+3*XYZmatrix[,3]), v76=9*XYZmatrix[,2]/(XYZmatrix[,1]+15*XYZmatrix[,2]+3*XYZmatrix[,3]))
}
Yxy2Yuv<-function(Yxymatrix){
if (!is.matrix(Yxymatrix)) Yxymatrix <- matrix(Yxymatrix, ncol=3,byrow=TRUE)
cbind(Y=Yxymatrix[,1], u76=4*Yxymatrix[,2]/(-2*Yxymatrix[,2]+12*Yxymatrix[,3]+3), v76=9*Yxymatrix[,3]/(-2*Yxymatrix[,2]+12*Yxymatrix[,3]+3)) 
}
Yuv2XYZ<-function(Yu.v.matrix) {
if (!is.matrix(Yu.v.matrix)) Yu.v.matrix <- matrix(Yu.v.matrix, ncol=3,byrow=TRUE)
cbind(X=9*Yu.v.matrix[,1]*Yu.v.matrix[,2]/(4*Yu.v.matrix[,3]), Y=Yu.v.matrix[,1], Z=3*Yu.v.matrix[,1]*(4-Yu.v.matrix[,2])/(4*Yu.v.matrix[,3])-5*Yu.v.matrix[,1])
}
Yuv2xy<-function(Yu.v.matrix){
if (!is.matrix(Yu.v.matrix)) Yu.v.matrix <- matrix(Yu.v.matrix, ncol=3,byrow=TRUE)
cbind(Y=Yu.v.matrix[,1], x=9*Yu.v.matrix[,2]/(6*Yu.v.matrix[,2]-16*Yu.v.matrix[,3]+12), y=4*Yu.v.matrix[,3]/(6*Yu.v.matrix[,2]-16*Yu.v.matrix[,3]+12))
}
# CIE-1976  u', v' to CIE-1931 x, y
# http://www.color-theory-phenomena.nl/10.03.htm
CIE1976uv2CIE1931xy<-function(uvmatrix) {
if (!is.matrix(uvmatrix)) uvmatrix <- matrix(uvmatrix, ncol=2,byrow=TRUE)
cbind( x = ( 27 * uvmatrix[,1] / 4 ) / ( ( 9 * uvmatrix[,1] / 2) - 12 * uvmatrix[,2] + 9 ), y = ( 3 * uvmatrix[,2]) / ( ( 9 * uvmatrix[,1] / 2) - 12 * uvmatrix[,2] + 9 )) 
}
#CIE-1976  u', v' to  CIE-1960 u, v
CIE1976uv2CIE1960uv<-function(uvmatrix){
if (!is.matrix(uvmatrix)) uvmatrix <- matrix(uvmatrix, ncol=2,byrow=TRUE)
cbind(u = uvmatrix[,1], v = ( 2 * uvmatrix[,2] / 3 ))
}
#CIE-1931 X, Y, Z to CIE-1931 x, y, z
CIE1931XYZ2CIE1931xyz<-function(XYZmatrix) {
if (!is.matrix(XYZmatrix)) XYZmatrix <- matrix(XYZmatrix, ncol=3,byrow=TRUE)
cbind(x = XYZmatrix[,1] / ( XYZmatrix[,1] + XYZmatrix[,2] + XYZmatrix[,3] ), 
y = XYZmatrix[,2] / ( XYZmatrix[,1] + XYZmatrix[,2] + XYZmatrix[,3] ), z = 1 - ( XYZmatrix[,1] / ( XYZmatrix[,1] + XYZmatrix[,2] + XYZmatrix[,3] ) + XYZmatrix[,2] / ( XYZmatrix[,1] + XYZmatrix[,2] + XYZmatrix[,3] ) ) )
}

#CIE-1931 X, Y, Z to CIE-1960 u, v
CIE1931XYZ2CIE1960uv<-function(XYZmatrix) {
if (!is.matrix(XYZmatrix)) XYZmatrix <- matrix(XYZmatrix, ncol=3,byrow=TRUE)
cbind(u = 4 * XYZmatrix[,1] / ( XYZmatrix[,1] + 15 * XYZmatrix[,2] + 3 * XYZmatrix[,3] ), 
v = 6 * XYZmatrix[,2] / ( XYZmatrix[,1] + 15 * XYZmatrix[,2] + 3 * XYZmatrix[,3] ))
}
#CIE-1931 x, y to CIE-1960 u, v
CIE1931xy2CIE1960uv<-function(xymatrix){
if (!is.matrix(xymatrix)) xymatrix <- matrix(xymatrix, ncol=2,byrow=TRUE)
cbind(u = 4 * xymatrix[,1] / ( -2 * xymatrix[,1] + 12 * xymatrix[,2] + 3 ), v = 6 * xymatrix[,2] / ( -2 * xymatrix[,1] + 12 * xymatrix[,2] + 3 ))
}
#CIE-1931 X, Y, Z to CIE-1976 u', v'
CIE1931XYZ2CIE1976uv<-function(XYZmatrix){
if (!is.matrix(XYZmatrix)) XYZmatrix <- matrix(XYZmatrix, ncol=3,byrow=TRUE)
cbind(u = 4 * XYZmatrix[,1] / ( XYZmatrix[,1] + 15 * XYZmatrix[,2] + 3 * XYZmatrix[,3] ), 
v = 9 * XYZmatrix[,2] / ( XYZmatrix[,1] + 15 * XYZmatrix[,2] + 3 * XYZmatrix[,3] ))
}
#CIE-1931 x, y to CIE-1976 u', v'
CIE1931xy2CIE1976uv<-function(xymatrix){
if (!is.matrix(xymatrix)) xymatrix <- matrix(xymatrix, ncol=2,byrow=TRUE)
cbind(u = 4 * xymatrix[,1] / ( -2 * xymatrix[,1] + 12 * xymatrix[,2] + 3 ), 
v = 9 * xymatrix[,2] / ( -2 * xymatrix[,1] + 12 * xymatrix[,2] + 3 ))
}
#CIE 1960 color space (CIE 1960 UCS)
# source: CIE 1964 color space, From Wikipedia, the free encyclopedia http://en.wikipedia.org/wiki/CIE_1960_color_space
Yxy2CIE1960UCS<-function(Yxymatrix){
if (!is.matrix(Yxymatrix)) Yxymatrix <- matrix(Yxymatrix, ncol=3,byrow=TRUE)
cbind(u=(0.4661*Yxymatrix[,1]+0.1593*Yxymatrix[,2])/(Yxymatrix[,2]-0.15735*Yxymatrix[,1]+0.2424), v=(0.6581**Yxymatrix[,2])/(Yxymatrix[,2]-0.15735*Yxymatrix[,1]+0.2424))
}
CIE1960UCS2xy<-function(uvMatrix){
if (!is.matrix(uvMatrix)) uvMatrix <- matrix(uvMatrix, ncol=2,byrow=TRUE)
cbind(x=3*uvMatrix[,1]/(2*uvMatrix[,1]-8*uvMatrix[,2]+4), y=2*uvMatrix[,2]/(2*uvMatrix[,1]-8*uvMatrix[,2]+4))
}
# CIE 1964 color space
# source: CIE 1964 color space, From Wikipedia, the free encyclopedia http://en.wikipedia.org/wiki/CIE_1964_color_space
CIE1960UCS2CIE1964<-function(uvYmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment())) 
{
if (is.null(dim(uvYmatrix))) if (length(uvYmatrix)>2) uvYmatrix<-matrix(uvYmatrix, ncol=3,byrow=TRUE)
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
xymatrix <- CIE1960UCS2xy(uvYmatrix[,1:2])
XYZmatrix <- xyY2XYZ(cbind(xymatrix,uvYmatrix[,3]))
xr <- XYZmatrix[,1] / Rx * 100
yr <- XYZmatrix[,2] / Ry * 100
zr <- XYZmatrix[,3] / Rz * 100
W <- 25*yr^(1/3)-17
U <- 13*W*(uvYmatrix[,1]-Rx)
V <- 13*uvYmatrix[,3]*(uvYmatrix[,2]-Ry)
cbind(U=U,V=V ,W=W )
}

XYZ2BVR <-function(XYZmatrix){# matrix for XYZ to Landolt BVR transformation 
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
XYZmatrix %*% matrix(c(0.200308,-0.017427,0.012508,0.013093,0.789019,-0.328110,-0.224010,-0.465060,1.381969),3,3,byrow=TRUE)
}

BVR2XYZ <-function(BVRmatrix){# matrix for Landolt BVR to XYZ transformation 
if (is.null(dim(BVRmatrix))) if (length(BVRmatrix)>2) BVRmatrix<-matrix(BVRmatrix, ncol=3,byrow=TRUE)
BVRmatrix %*% matrix(c(4.961444,0.293121,0.902866, 0.096640,1.479324,0.513487,-0.021961,0.348571,0.837347),3,3,byrow=TRUE)
}

RGB2LEF<-function(RGBmatrix)
{# based on: Kang, Henry R, 2006, Computational color technology, Spie Press Bellingham
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-t(matrix(RGBmatrix, ncol=3,byrow=TRUE))
matrix(c(2/3,2/3,2/3,1,-1/2,-1/2,0,sqrt(3)/2,-sqrt(3)/2),3,3,byrow=TRUE) %*% RGBmatrix
}

LEF2RGB<-function(LEFmatrix)
{# based on: Kang, Henry R, 2006, Computational color technology, Spie Press Bellingham
if (is.null(dim(LEFmatrix))) if (length(LEFmatrix)>2) LEFmatrix<-t(matrix(LEFmatrix, ncol=3,byrow=TRUE))
matrix(c(1/2,2/3,0,1/2,-1/3,1/sqrt(3),1/2,-1/3,-1/sqrt(3)),3,3,byrow=TRUE) %*% LEFmatrix
}

RGB2YPbPr<-function(RGBmatrix)
{# based on: ColorObject.pm Graphics/ColorObject version 0.5.0, Copyright 2003-2005 by Alex Izvorski, (Portions Copyright 2001-2003 by Alfred Reibenschuh)
# reference: http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.txt
# Y is [0..1], Pb and Pr are [-0.5..0.5]
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
matrix(c(0.299,0.587,0.114,-0.168736,-0.331264,0.5,0.5,-0.418688,-0.081312),3,3,byrow=TRUE) %*% t(RGBmatrix)
}

YPbPr2RGB<-function(YPbPrmatrix)
{# based on: ColorObject.pm Graphics/ColorObject version 0.5.0, Copyright 2003-2005 by Alex Izvorski, (Portions Copyright 2001-2003 by Alfred Reibenschuh)
# reference: http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.txt
# result is NTSC non-linear rgb
# Pb and Pr are [-0.5..0.5], Y is [0..1]
if (is.null(dim(YPbPrmatrix))) if (length(YPbPrmatrix)>2) YPbPrmatrix<-matrix(YPbPrmatrix, ncol=3,byrow=TRUE)
matrix(c(1.0,0.0,1.402,1.0,-0.344136,-0.714136,1.0,1.772,0.0),3,3,byrow=TRUE) %*% t(YPbPrmatrix)
}

RGB2YCbCr<-function(RGBmatrix)
{# based on: ColorObject.pm Graphics/ColorObject version 0.5.0, Copyright 2003-2005 by Alex Izvorski, (Portions Copyright 2001-2003 by Alfred Reibenschuh)
# reference: http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.txt
# input should be NTSC non-linear rgb
# Y is [16..235], Cb and Cr are [16..239]
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
matrix(c(65.481,128.553,24.966,-37.797,-74.203,112.0 ,112.0,-93.786,-18.214),3,3,byrow=TRUE) %*% t(RGBmatrix)
}

YCbCr2RGB<-function(YPbPrmatrix)
{# based on: ColorObject.pm Graphics/ColorObject version 0.5.0, Copyright 2003-2005 by Alex Izvorski, (Portions Copyright 2001-2003 by Alfred Reibenschuh)
# reference: http://www.poynton.com/notes/colour_and_gamma/ColorFAQ.txt
# Y is [16..235], Cb and Cr are [16..239]
if (is.null(dim(YPbPrmatrix))) if (length(YPbPrmatrix)>2) YPbPrmatrix<-matrix(YPbPrmatrix, ncol=3,byrow=TRUE)
matrix(c(0.00456621,0.0,0.00625893, 0.00456621,-0.00153632,-0.00318811,0.00456621, 0.00791071, 0.0),3,3,byrow=TRUE) %*% t(YPbPrmatrix)
}

RGB2YIQ<-function(RGBmatrix)
{
# assumes RGBmatrix in the range [0,1]
# returns YIQmatrix in the range [-1,1]
# based on Color space conversion by Sophie Kirschner, http://www.blitzbasic.com/codearcs/codearcs.php?code=2953
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
YIQmatrix<-matrix(0,dim(RGBmatrix)[1],3)
YIQmatrix[,1]<- 0.299000*RGBmatrix[,1]+0.587000*RGBmatrix[,2]+0.114000*RGBmatrix[,3]
YIQmatrix[,2]<- 0.595716*RGBmatrix[,1]-0.274453*RGBmatrix[,2]-0.321264*RGBmatrix[,3]
YIQmatrix[,3]<- 0.211456*RGBmatrix[,1]-0.522591*RGBmatrix[,2]+0.311350*RGBmatrix[,3]
YIQmatrix
}

YIQ2RGB<-function(YIQmatrix)
{
# assumes YIQmatrix in the range [-1,1]
# returns RGBmatrix in the range [0,1]
# based on Color space conversion by Sophie Kirschner, http://www.blitzbasic.com/codearcs/codearcs.php?code=2953
if (is.null(dim(YIQmatrix))) if (length(YIQmatrix)>2) YIQmatrix<-matrix(YIQmatrix, ncol=3,byrow=TRUE)
RGBmatrix<-matrix(0,dim(YIQmatrix)[1],3)
RGBmatrix[,1]<- YIQmatrix[,1]+0.9563*YIQmatrix[,2]+0.6210*YIQmatrix[,3]
RGBmatrix[,2]<- YIQmatrix[,1]-0.2721*YIQmatrix[,2]-0.6474*YIQmatrix[,3]
RGBmatrix[,3]<- YIQmatrix[,1]-1.1070*YIQmatrix[,2]+1.7046*YIQmatrix[,3]
RGBmatrix
}

RGB2YUV<-function(RGBmatrix)
{
# assumes RGBmatrix in the range [0,1]
# returns YUVmatrix in the range [-1,1]
# based on Color space conversion by Madk, Sophie Kirschner, http://www.blitzbasic.com/codearcs/codearcs.php?code=2953
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
yuv_uconst <- 1.0/0.436
yuv_vconst <- 1.0/0.615
YUVmatrix<-matrix(0,dim(RGBmatrix)[1],3)
YUVmatrix[,1]<- 0.299*RGBmatrix[,1]+0.587*RGBmatrix[,2]+0.114*RGBmatrix[,3]
YUVmatrix[,2]<- (0.492*(RGBmatrix[,3]-YUVmatrix[,1]))*yuv_uconst
YUVmatrix[,3]<- (0.877*(RGBmatrix[,1]-YUVmatrix[,1]))*yuv_vconst
YUVmatrix[,1]<- YUVmatrix[,1]*2-1
YUVmatrix
}

YUV2RGB<-function(YUVmatrix)
{
# assumes YUVmatrix in the range [-1,1]
# returns RGBmatrix in the range [0,1]
# based on Color space conversion by Madk, Sophie Kirschner, http://www.blitzbasic.com/codearcs/codearcs.php?code=2953
if (is.null(dim(YUVmatrix))) if (length(YUVmatrix)>2) YUVmatrix<-matrix(YUVmatrix, ncol=3,byrow=TRUE)
yuv_uconst <- 1.0/0.436
yuv_vconst <- 1.0/0.615
ly <- (YUVmatrix[,1]+1)/2.0 #YUVmatrix[,2]/yuv_uconst,YUVmatrix[,3]/yuv_vconst)
RGBmatrix<-matrix(0,dim(YUVmatrix)[1],3)
RGBmatrix[,1]<- ly +1.140*YUVmatrix[,3]/yuv_vconst
RGBmatrix[,2]<- ly -0.395*YUVmatrix[,2]/yuv_uconst - 0.581*YUVmatrix[,3]/yuv_vconst
RGBmatrix[,3]<- ly +2.032*YUVmatrix[,2]/yuv_uconst
RGBmatrix
}

RGB2HSL<-function(RGBmatrix)
{ # code adapted from: easyrgb Color conversion math and formulas http://www.easyrgb.com/
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
var.R <- ( RGBmatrix[,1] / 255 )                     #RGB values <- From 0 to 255
var.G <- ( RGBmatrix[,2] / 255 )
var.B <- ( RGBmatrix[,3] / 255 )
var.Min <- apply(RGBmatrix / 255,1,min)    #Min. value of RGB
var.Max <- apply(RGBmatrix / 255,1,max)    #Max. value of RGB
del.Max <- var.Max - var.Min             #Delta RGB value
L <- ( var.Max + var.Min ) / 2
H<-vector('numeric',length(L))
S<-H
LS1.max<-which((del.Max != 0) & (L < 0.5))
LS2.max<-which((del.Max != 0) & (L >= 0.5))
if (length(LS1.max)>0) S[LS1.max] <- del.Max[LS1.max] / ( var.Max[LS1.max] + var.Min[LS1.max] )
if (length(LS2.max)>0) S[LS2.max] <- del.Max[LS2.max] / ( 2 - var.Max[LS2.max] - var.Min[LS2.max] )
del.R <- ( ( ( var.Max - var.R ) / 6 ) + ( del.Max / 2 ) ) / del.Max
del.G <- ( ( ( var.Max - var.G ) / 6 ) + ( del.Max / 2 ) ) / del.Max
del.B <- ( ( ( var.Max - var.B ) / 6 ) + ( del.Max / 2 ) ) / del.Max
R.max<-which((del.Max != 0) & (var.R == var.Max))
G.max<-which((del.Max != 0) & (var.G == var.Max))
B.max<-which((del.Max != 0) & (var.B == var.Max))
if (length(R.max)>0) H[R.max] <- del.B[R.max] - del.G[R.max]
if (length(G.max)>0) H[G.max] <- ( 1 / 3 ) + del.R[G.max] - del.B[G.max]
if (length(B.max)>0) H[B.max] <- ( 2 / 3 ) + del.G[B.max] - del.R[B.max]
Hl0<-which((del.Max != 0) & (H < 0))
if (length(Hl0)>0) H[Hl0] <- H[Hl0] + 1
Hg1<-which((del.Max != 0) & (H > 1))
if (length(Hg1)>0) H[Hg1] <- H[Hg1] - 1
cbind(H=H,S=S,L=L)
}

HSL2RGB<-function(HSLmatrix)
{
if (is.null(dim(HSLmatrix))) if (length(HSLmatrix)>2) HSLmatrix<-matrix(HSLmatrix, ncol=3,byrow=TRUE)
R<-vector('numeric',dim(HSLmatrix)[2])
G<-R
B<-R
S0<-which(HSLmatrix[,2] == 0)
if (length(S0)>0) {
R[S0] <- HSLmatrix[S0,1] * 255
G[S0] <- HSLmatrix[S0,2] * 255
B[S0] <- HSLmatrix[S0,3] * 255
}
S1<-which(HSLmatrix[,2] != 0)
if (length(S1)>0) {
var.2<-vector('numeric',length(S1))
Ll05<-which((HSLmatrix[S1,3] < 0.5))
if (length(Ll05)>0) var.2 <- HSLmatrix[Ll05,3] * ( 1 + HSLmatrix[Ll05,2] )
Lg05<-which((HSLmatrix[S1,3] >= 0.5))
if (length(Lg05)>0) var.2 <- ( HSLmatrix[Lg05,2] + HSLmatrix[Lg05,3] ) - ( HSLmatrix[Lg05,2] * HSLmatrix[Lg05,3] )
var.1 <- 2 * HSLmatrix[S1,3] - var.2
R[S1] <- 255 * Hue.2.RGB( var.1, var.2, HSLmatrix[S1,1] + ( 1 / 3 ) ) 
G[S1] <- 255 * Hue.2.RGB( var.1, var.2, HSLmatrix[S1,1] )
B[S1] <- 255 * Hue.2.RGB( var.1, var.2, HSLmatrix[S1,1] - ( 1 / 3 ) )
}
HSLret<-cbind(R=R,G=G,B=B)
HSLret[1:dim(HSLmatrix)[1],]
}

Hue.2.RGB<-function(v1, v2, vH)
{
S0<-which(vH < 0)
vH[S0] <- vH[S0] + 1
S1<-which(vH > 0)
vH[S1] <- vH[S1] - 1
h6<-which(( 6 * vH ) < 1)
v1[h6] <- ( v1[h6] + ( v2[h6] - v1[h6] ) * 6 * vH[h6] )
h2<-which(( 2 * vH ) < 1)
v1[h2] <- ( v1[h2] + ( v2[h2] - v1[h2] ) * 6 * vH[h2] )
h3<-which(( 3 * vH ) < 2)
v1[h3] <- ( v1[h3] + ( v2[h3] - v1[h3] ) * 6 * vH[h3] )
v1
}

RGB2HSV<-function(RGBmatrix)
{ # code adapted from: easyrgb Color conversion math and formulas http://www.easyrgb.com/
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
var.R <- ( RGBmatrix[,1] / 255 )                     #RGB values <- From 0 to 255
var.G <- ( RGBmatrix[,2] / 255 )
var.B <- ( RGBmatrix[,3] / 255 )
var.Min <- apply(RGBmatrix / 255,1,min)    #Min. value of RGB
var.Max <- apply(RGBmatrix / 255,1,max)    #Max. value of RGB
del.Max <- var.Max - var.Min             #Delta RGB value
V <- var.Max
H<-vector('numeric',length(V))
S<- H
Sn0<-which(del.Max != 0)
if (length(Sn0)>0) S[Sn0]<- del.Max[Sn0] / var.Max[Sn0]
del.R <- ( ( ( var.Max - var.R ) / 6 ) + ( del.Max / 2 ) ) / del.Max
del.G <- ( ( ( var.Max - var.G ) / 6 ) + ( del.Max / 2 ) ) / del.Max
del.B <- ( ( ( var.Max - var.B ) / 6 ) + ( del.Max / 2 ) ) / del.Max
R.max<-which((del.Max != 0) & (var.R == var.Max))
G.max<-which((del.Max != 0) & (var.G == var.Max))
B.max<-which((del.Max != 0) & (var.B == var.Max))
if (length(R.max)>0) H[R.max] <- del.B[R.max] - del.G[R.max]
if (length(G.max)>0) H[G.max] <- ( 1 / 3 ) + del.R[G.max] - del.B[G.max]
if (length(B.max)>0) H[B.max] <- ( 2 / 3 ) + del.G[B.max] - del.R[B.max]
Hl0<-which((del.Max != 0) & (H < 0))
if (length(Hl0)>0) H[Hl0] <- H[Hl0] + 1
Hg1<-which((del.Max != 0) & (H > 1))
if (length(Hg1)>0) H[Hg1] <- H[Hg1] - 1
cbind(H=H,S=S,V=V)
}

HSV2RGB<-function(HSVmatrix)
{
if (is.null(dim(HSVmatrix))) if (length(HSVmatrix)>2) HSVmatrix<-matrix(HSVmatrix, ncol=3,byrow=TRUE)
R<-vector('numeric',dim(HSVmatrix)[2])
G<-R
B<-R
S0<-which(HSVmatrix[,2] == 0)
if (length(S0)>0) {
R[S0] <- HSVmatrix[S0,1] * 255
G[S0] <- HSVmatrix[S0,2] * 255
B[S0] <- HSVmatrix[S0,3] * 255
}
Sn0<-which(HSVmatrix[,2] != 0)
if (length(Sn0)>0) {
var.h<-vector('numeric',dim(HSVmatrix)[2])
var.i<-vector('numeric',dim(HSVmatrix)[2])
var.1<-vector('numeric',dim(HSVmatrix)[2])
var.2<-vector('numeric',dim(HSVmatrix)[2])
var.3<-vector('numeric',dim(HSVmatrix)[2])
var.r<-vector('numeric',dim(HSVmatrix)[2])
var.g<-vector('numeric',dim(HSVmatrix)[2])
var.b<-vector('numeric',dim(HSVmatrix)[2])
var.h[Sn0] <- HSVmatrix[Sn0,1] * 6
var.i[Sn0] <- floor( var.h[Sn0] )
var.1[Sn0] <- HSVmatrix[Sn0,3] * ( 1 - HSVmatrix[Sn0,2] )
var.2[Sn0] <- HSVmatrix[Sn0,3] * ( 1 - HSVmatrix[Sn0,2] * ( var.h[Sn0] - var.i[Sn0] ) )
var.3[Sn0] <- HSVmatrix[Sn0,3] * ( 1 - HSVmatrix[Sn0,2] * ( 1 - ( var.h[Sn0] - var.i[Sn0] ) ) )
ieq0<-which((HSVmatrix[,2] != 0) & (var.i == 0))
if (length(ieq0)>0) {
var.r[ieq0] <- HSVmatrix[ieq0]
var.g[ieq0] <- var.3[ieq0]
var.b[ieq0] <- var.1[ieq0]
}
ieq1<-which((HSVmatrix[,2] != 0) & (var.i == 1))
if (length(ieq1)>0) {
var.r[ieq1] <- var.2[ieq1]
var.g[ieq1] <- HSVmatrix[ieq1]
var.b[ieq1] <- var.1[ieq1]
}
ieq2<-which((HSVmatrix[,2] != 0) & (var.i == 2))
if (length(ieq2)>0) {
var.r[ieq2] <- var.1[ieq2]
var.g[ieq2] <- HSVmatrix[ieq2]
var.b[ieq2] <- var.3[ieq2]
}
ieq3<-which((HSVmatrix[,2] != 0) & (var.i == 3))
if (length(ieq3)>0) {
var.r[ieq3] <- var.1[ieq3]
var.g[ieq3] <- var.2[ieq3]
var.b[ieq3] <- HSVmatrix[ieq3]
}
ieq4<-which((HSVmatrix[,2] != 0) & (var.i == 4))
if (length(ieq4)>0) {
var.r[ieq4] <- var.3[ieq4]
var.g[ieq4] <- var.1[ieq4]
var.b[ieq4] <- HSVmatrix[ieq4]
}
ieq5<-which((HSVmatrix[,2] != 0) & !any(var.i %in% 1:4))
if (length(ieq5)>5) {
var.r[ieq5] <- HSVmatrix[ieq5]
var.g[ieq5] <- var.1[ieq5]
var.b[ieq5] <- var.2[ieq5]
}
R[Sn0] <- var.r[Sn0] * 255                  #RGB results <- From 0 to 255
G[Sn0] <- var.g[Sn0] * 255
B[Sn0] <- var.b[Sn0] * 255
}
HSVret<- cbind(R=R,G=G,B=B)
HSVret[1:dim(HSVmatrix)[1],]
}

RGB2CMY<-function(RGBmatrix)
{#RGB values <- From 0 to 255
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
C <- 1 - ( RGBmatrix[,1] / 255 )
M <- 1 - ( RGBmatrix[,2] / 255 )
Y <- 1 - ( RGBmatrix[,3] / 255 )
cbind(C=C,M=M,Y=Y)
}

CMY2RGB<-function(CMYmatrix)
{#CMY values <- From 0 to 1
if (is.null(dim(CMYmatrix))) if (length(CMYmatrix)>2) CMYmatrix<-matrix(CMYmatrix, ncol=3,byrow=TRUE)
R <- ( 1 - CMYmatrix[,1] ) * 255
G <- ( 1 - CMYmatrix[,2] ) * 255
B <- ( 1 - CMYmatrix[,3] ) * 255
cbind(R=R,G=G,B=B)
}

CMY2CMYK<-function(CMYmatrix)
{#CMY values <- From 0 to 1
if (is.null(dim(CMYmatrix))) if (length(CMYmatrix)>2) CMYmatrix<-matrix(CMYmatrix, ncol=3,byrow=TRUE)
var.K <- rep(1,dim(CMYmatrix)[1])
CK<-which(CMYmatrix[,1] < var.K)
if (length(CK)>0) var.K[CK] <- CMYmatrix[CK,1]
CM<-which(CMYmatrix[,2] < var.K)
if (length(CM)>0) var.K[CM] <- CMYmatrix[CM,2]
CY<-which(CMYmatrix[,3] < var.K)
if (length(CY)>0) var.K[CY] <- CMYmatrix[CY,3]
C <- ( CMYmatrix[,1] - var.K ) / ( 1 - var.K )
M <- ( CMYmatrix[,2] - var.K ) / ( 1 - var.K )
Y <- ( CMYmatrix[,3] - var.K ) / ( 1 - var.K )
cbind(C=C,M=M,Y=Y,K=var.K)
}

CMYK2CMY<-function(CMYKmatrix)
{#CMYK values <- From 0 to 1
if (is.null(dim(CMYKmatrix))) if (length(CMYKmatrix)>2) CMYKmatrix<-matrix(CMYKmatrix, ncol=4,byrow=TRUE)
C <- ( CMYKmatrix[,1] * ( 1 - CMYKmatrix[,4] ) + CMYKmatrix[,4] )
M <- ( CMYKmatrix[,2] * ( 1 - CMYKmatrix[,4] ) + CMYKmatrix[,4] )
Y <- ( CMYKmatrix[,3] * ( 1 - CMYKmatrix[,4] ) + CMYKmatrix[,4] )
cbind(C=C,M=M,Y=Y)
}

XYZ2HunterLab<-function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{ # adapted from: easyrgb Color conversion math and formulas http://www.easyrgb.com/
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
xr <- XYZmatrix[,1] / Rx * 100
yr <- XYZmatrix[,2] / Ry * 100
zr <- XYZmatrix[,3] / Rz * 100
L<-10 * sqrt( yr )
a<-17.5 * ( ( ( 1.02 * xr ) - yr ) / L )
b<-7 * ( ( yr - ( 0.847 * zr ) ) / L )
cbind(L=L,a=a,b=b)
}

HunterLab2XYZ<-function(HunterLabmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{ # adapted from: easyrgb Color conversion math and formulas http://www.easyrgb.com/
if (is.null(dim(HunterLabmatrix))) if (length(HunterLabmatrix)>2) HunterLabmatrix<-matrix(HunterLabmatrix, ncol=3,byrow=TRUE)
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
Y<-HunterLabmatrix[,1] / 10
X<-HunterLabmatrix[,2] / 17.5 * Y
Z<-HunterLabmatrix[,3] / 7 * Y
Y<-Y^2
X<- ( X + Y ) / 1.02
Z<- -( Z - Y ) / 0.847
X <- X * Rx
Y <- Y * Ry
Z <- Z * Rz
cbind(X=X,Y=Y,Z=Z)
}

NickersonColorDifference<-function(MunsellHVC1, MunsellHVC2) {
# deltaE = 2/5CdH + 6dV + 3dC
# COLOR TECHNOLOGY in the textile industry Second Edition
m1<-MunsellHVC(MunsellHVC1)
m2<-MunsellHVC(MunsellHVC2)
C<- as.numeric(m1[,'C'])
dC<- C - as.numeric(m2[,'C'])
dV<- as.numeric(m1[,'V']) - as.numeric(m2[,'V'])
dH<-huedegree(m1[,'H']) - huedegree(m2[,'H'])
as.numeric(2*C*dH/5 + 6*dV + 3*dC)
}

CCT2XYZ<-function (CCTmatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
C1 <- 2.0 * pi * 6.626176 * 2.99792458 * 2.99792458
C2 <- (6.626176 * 2.99792458) / 1.380662
nmSeq<-seq(360,830,by=5)
dWavelengthM <- nmSeq * 1.0e-3
dWavelengthM5 <- dWavelengthM * dWavelengthM * dWavelengthM * dWavelengthM * dWavelengthM
blackbody <- C1 / (dWavelengthM5 * 1.0e-12 * (exp(C2 / (CCTmatrix * dWavelengthM * 1.0e-3)) - 1.0))
XYZmatrix<-blackbody * get("ciexyz31", envir = environment())[nmSeq-359,c('xbar','ybar','zbar')]
XYZmatrix[,1]<-XYZmatrix[,1] / XYZmatrix[,2]
XYZmatrix[,3]<-XYZmatrix[,3] / XYZmatrix[,2]
XYZmatrix[,2]<-1.0
}

Lab2LCHab<-function(LabMatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(LabMatrix))) if (length(LabMatrix)>2) LabMatrix<-matrix(LabMatrix, ncol=3,byrow=TRUE)
LCHabmatrix<-LabMatrix
LCHabmatrix[,2]<-sqrt(LabMatrix[,2]^2 + LabMatrix[,3]^2)
LCHabmatrix[,3]<-180.0 * atan2(LabMatrix[,3], LabMatrix[,2]) / pi
if (LCHabmatrix[,3] < 0.0) LCHabmatrix[,3]<-LCHabmatrix[,3]+360
LCHabmatrix
}

LCHab2Lab<-function(LCHabmatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(LCHabmatrix))) if (length(LCHabmatrix)>2) LCHabmatrix<-matrix(LCHabmatrix, ncol=3,byrow=TRUE)
LabMatrix<-LCHabmatrix
LabMatrix[,2]<-LCHabmatrix[,2] * cos(LCHabmatrix[,3] * pi / 180.0)
LabMatrix[,3]<-LCHabmatrix[,2] * sin(LCHabmatrix[,3] * pi / 180.0)
LabMatrix
}

Luv2LCHuv<-function(LuvMatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(LuvMatrix))) if (length(LuvMatrix)>2) LuvMatrix<-matrix(LuvMatrix, ncol=3,byrow=TRUE)
LCHuvmatrix<-LuvMatrix
LCHuvmatrix[,2]<-sqrt(LuvMatrix[,2]^2 + LuvMatrix[,3]^2)
LCHuvmatrix[,3]<-180.0 * atan2(LuvMatrix[,3], LuvMatrix[,2]) / pi
if (LCHuvmatrix[,3] < 0.0) LCHuvmatrix[,3]<-LCHuvmatrix[,3]+360
LCHuvmatrix
}

LCHuv2Luv<-function(LCHuvmatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(LCHuvmatrix))) if (length(LCHuvmatrix)>2) LCHuvmatrix<-matrix(LCHuvmatrix, ncol=3,byrow=TRUE)
LuvMatrix<-LCHuvmatrix
LuvMatrix[,2]<-LCHuvmatrix[,2] * cos(LCHuvmatrix[,3] * pi / 180.0)
LuvMatrix[,3]<-LCHuvmatrix[,2] * sin(LCHuvmatrix[,3] * pi / 180.0)
LuvMatrix
}

DIN6167.YellownessIndex<-function(XYZmatrix,illuminant='C',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment())) 
{ # source: Basic equations for optical properties, 
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
RxRyRz<-XYZ2RxRyRz(XYZmatrix,illuminant,observer,RefWhite)
100*( RxRyRz[,1] - RxRyRz[,3] ) / RxRyRz[,2]
}

ASTM.D1925.YellownessIndex<-function(XYZmatrix) 
{
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
(100*(1.28*XYZmatrix[,1]-1.06*XYZmatrix[,3]))/XYZmatrix[,2]
}
#ASTM D 1925 Yellowness Index for Plastics
#X, Y and Z are the tri-stimulus values for the calculated for illuminant C
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012, pp. 6

ASTM.E313.YellownessIndex<-function(XYZmatrix)
{
#ASTM D 1925 Yellowness Index
#X, Y and Z are the tri-stimulus values for the calculated for illuminant C
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012, pp. 6
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
(100*(1-0.847*XYZmatrix[,3]))/XYZmatrix[,2]
}

ASTM.E313.Whiteness<-function(XYZmatrix)
{
#ASTM D 1925 Whiteness Index
#X, Y and Z are the tri-stimulus values for the calculated for illuminant C
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012, pp. 6
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
3.388*XYZmatrix[,3]-3*XYZmatrix[,2]
}

CIE.Whiteness<-function(xyYmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# strictly for D65 and 2 or 10 deg observer
#Values bigger than 100 indicate a bluish white
#Values smaller than 100 indicate a yellowish white
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(xyYmatrix))) if (length(xyYmatrix)>2) xyYmatrix<-matrix(xyYmatrix, ncol=3,byrow=TRUE)
Rrgbwhitergb<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
xr<-Rrgbwhitex / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
yr<-Rrgbwhitey / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
xyYmatrix[,3]+800*(xr-xyYmatrix[,1])+1700*(yr-xyYmatrix[,2])
}

Berger59.Whiteness<-function(xyYmatrix,illuminant='C',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(xyYmatrix))) if (length(xyYmatrix)>2) xyYmatrix<-matrix(xyYmatrix, ncol=3,byrow=TRUE)
Rrgbwhitergb<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
xr<-Rrgbwhitex / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
zr<-Rrgbwhitez / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
0.333*xyYmatrix[,2]+125*(xyYmatrix[,3]/zr)-125*(xyYmatrix[,1]/xr)
}

Stensby68.Whiteness<-function(LabHunterMatrix) {
if (is.null(dim(LabHunterMatrix))) if (length(LabHunterMatrix)>2) LabHunterMatrix<-matrix(LabHunterMatrix, ncol=3,byrow=TRUE)
LabHunterMatrix[,1]-3*LabHunterMatrix[,3]+3*LabHunterMatrix[,2]
}
#L, a and b are Hunter Color Coordinates
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012

Taube60.Whiteness<-function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
Rrgbwhitergb<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
400*XYZmatrix[,3]/Rrgbwhitez-3*XYZmatrix[,2]
}

Hunter60.WhitenessIndex<-function(LabHunterMatrix) {
#L, a and b are Hunter Color Coordinates
# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(LabHunterMatrix))) if (length(LabHunterMatrix)>2) LabHunterMatrix<-matrix(LabHunterMatrix, ncol=3,byrow=TRUE)
LabHunterMatrix[,1]-3*LabHunterMatrix[,3]
}

GanzGrieser.Whiteness<-function(xyYmatrix)
{# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(xyYmatrix))) if (length(xyYmatrix)>2) xyYmatrix<-matrix(xyYmatrix, ncol=3,byrow=TRUE)
P <- -1868.322
Q <- -3695.690
C <- 1809.441
xyYmatrix[,3]+P*xyYmatrix[,1]+Q*xyYmatrix[,2]+C
}

GanzGrieser.Tint<-function(xyYmatrix)
{# Color iQC and Color iMatch Color Calculations Guide Version 8.0 July 2012
if (is.null(dim(xyYmatrix))) if (length(xyYmatrix)>2) xyYmatrix<-matrix(xyYmatrix, ncol=3,byrow=TRUE)
m <- -1001.223
n <- 748.366
k <- 68.261
m*xyYmatrix[,1]+n*xyYmatrix[,2]+k
}

CIETint<-function(xymatrix,illuminant='D65',observer=2)
{#Tint indices: CIE Tint and ASTM E313 Tint
#Tint E313 = Tint CIE = Tx (x_n - x) - 650 (y_n - y)
#CIE Publication 15:2004, "Colorimetry"
#ASTM E313, "Standard Practice for Calculating Yellowness and Whiteness Indices from Instrumentally Measured Color Coordinates"
if (is.null(dim(xymatrix))) if (length(xymatrix)>2) xymatrix<-matrix(xymatrix, ncol=3,byrow=TRUE)
x<-array(matrix(c(1000,1000,1000,900,900,900,0.3101,0.3457,0.3127,0.3104,0.3477,0.3138,0.3161,0.3585,0.3290,0.3191,0.3595,0.3310),3,6),c(3,2,3))
dimnames(x)<-list(c('C','D50','D65'), c('2','10'), c('T','x','y'))
Txy<-x[illuminant,observer,]
Txy[1]*(Txy[2]-xymatrix[,1])-650*(Txy[3]-xymatrix[,2])
}

DominantWavelength<-function(xyYmatrix, illuminant='D65',observer=2,RefWhiteIllum=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(xyYmatrix))) if (length(xyYmatrix)>2) xyYmatrix<-matrix(xyYmatrix, ncol=3,byrow=TRUE)
Rrgbwhitergb<-RefWhiteIllum[which(RefWhiteIllum[["Illuminant"]]==illuminant ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
xr<-Rrgbwhitex / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
yr<-Rrgbwhitey / (Rrgbwhitex + Rrgbwhitey + Rrgbwhitez)
count <- 1
tArray <- vector('numeric',2)# t
wArray <- vector('numeric',2)# wavelength
cArray <- vector('numeric',2)# cycle
a <- xyYmatrix[,1] - xr
b <- xyYmatrix[,2] - yr
if ((a >= -0.000001) & (a <= 0.000001) & (b >= -0.000001) & (b <= 0.000001)) stop('Error, (x, y) must be different from (xr, yr) at least by 1e-6')
for (nm in seq(360, 830, by=5))
{
i1 <- (nm - 360) / 5
i2 <- ifelse(nm == 830, 0, i1 + 1)
nm2 <- 5 * i2 + 360
x1 <- get("ciexyz31", envir = environment())[i1*5 + 1,'xbar'] / (get("ciexyz31", envir = environment())[i1*5 + 1,'xbar']+get("ciexyz31", envir = environment())[i1*5 + 1,'ybar']+get("ciexyz31", envir = environment())[i1*5 + 1,'zbar'])
y1 <- get("ciexyz31", envir = environment())[i1*5 + 1,'ybar'] / (get("ciexyz31", envir = environment())[i1*5 + 1,'xbar']+get("ciexyz31", envir = environment())[i1*5 + 1,'ybar']+get("ciexyz31", envir = environment())[i1*5 + 1,'zbar'])
x2 <- get("ciexyz31", envir = environment())[i2*5 + 1,'xbar'] / (get("ciexyz31", envir = environment())[i2*5 + 1,'xbar']+get("ciexyz31", envir = environment())[i2*5 + 1,'ybar']+get("ciexyz31", envir = environment())[i2*5 + 1,'zbar'])
y2 <- get("ciexyz31", envir = environment())[i2*5 + 1,'ybar'] / (get("ciexyz31", envir = environment())[i2*5 + 1,'xbar']+get("ciexyz31", envir = environment())[i2*5 + 1,'ybar']+get("ciexyz31", envir = environment())[i2*5 + 1,'zbar'])
C <- x1 - xr
d <- y1 - yr
E <- x2 - x1
f <- y2 - y1
s <- (a * d - b * C) / (b * E - a * f)
#cat(nm,i1,i2,nm2,x1,y1,x2,y2,C,d,E,f,s,'\n')
if ((s < 0.0) | (s >= 1.0)) next
Tval <- ifelse(abs(a) >= abs(b),((E * s + C) / a), ((f * s + d) / b))
tArray[count] <- Tval
cArray[count] <- nm
wArray[count] <- (nm2 - nm) * s + nm
count <- count + 1
}
if ((cArray[2] == 830) & (tArray[2] > 0.0)) dominantWavelength <- -wArray[1] else dominantWavelength <- 
ifelse(tArray[1] >= 0.0, wArray[1], wArray[2])
dominantWavelength
}

XYZ2RGB<-function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()),RGBModel='sRGB',RefWhiteRGB=get("whitepointsRGB", envir = environment()),gamma=NA,RefWhiteIllum=get("XYZperfectreflectingdiffuser", envir = environment()),CAT='Bradford',CATarray=get("ChromaticAdaptation", envir = environment())) 
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(XYZmatrix))) if (length(XYZmatrix)>2) XYZmatrix<-matrix(XYZmatrix, ncol=3,byrow=TRUE)
CATmatrix<-CATarray[, , CAT, 'direct']
Rrgb<-RefWhiteRGB[which(RefWhiteRGB[["description"]]==RGBModel ),]
Rillum<-unlist(Rrgb[ 'whitepointilluminant' ])
Rgamma<-unlist(Rrgb['gamma'])
Rrgbwhitergb<-RefWhiteIllum[which(RefWhiteIllum[["Illuminant"]]==Rillum ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
if (is.na(gamma)) gamma<-Rgamma
if (RGBModel=='sRGB') gamma<- -gamma
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
xr<-unlist(Rrgb['xRed'])
yr<-unlist(Rrgb['yRed'])
xg<-unlist(Rrgb['xGreen'])
yg<-unlist(Rrgb['yGreen'])
xb<-unlist(Rrgb['xBlue'])
yb<-unlist(Rrgb['yBlue'])
m<-matrix(c(xr/yr,xg/yg,xb/yb,1.0,1.0,1.0,(1.0-xr-yr)/yr,(1.0-xg-yg)/yg,(1.0-xb-yb)/yb),3,3,byrow=TRUE)
mi<-solve(m)
sr<- mi %*% matrix(c(Rrgbwhitex,Rrgbwhitey,Rrgbwhitez),3,1)
MtxXYZ2RGB<- t(matrix(sr,3,3,byrow=TRUE) * m)
MtxXYZ2RGB<- solve(MtxXYZ2RGB)
Asz<-CATmatrix %*% matrix(c(Rx,Ry,Rz),3,1,byrow=T)
Adz<-CATmatrix %*% matrix(c(Rrgbwhitex,Rrgbwhitey,Rrgbwhitez),3,1,byrow=T)
XYZ<-XYZmatrix %*% t(CATmatrix)
xyz1<- XYZ * c(Adz/Asz)
xyz2<- xyz1 %*% t(solve(CATmatrix))
xyz3<- xyz2 %*% MtxXYZ2RGB
XYZ<-apply(xyz3,1:2,function(x) Compand(x, gamma))
XYZ
}

RGB2XYZ<-function(RGBmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()),RGBModel='sRGB',RefWhiteRGB=get("whitepointsRGB", envir = environment()),gamma=NA,
RefWhiteIllum=get("XYZperfectreflectingdiffuser", envir = environment()),CAT='Bradford',CATarray=get("ChromaticAdaptation", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
CATmatrix<-CATarray[, , CAT, 'direct']
Rrgb<-RefWhiteRGB[which(RefWhiteRGB[["description"]]==RGBModel ),]
Rillum<-unlist(Rrgb[ 'whitepointilluminant' ])
Rgamma<-unlist(Rrgb['gamma'])
Rrgbwhitergb<-RefWhiteIllum[which(RefWhiteIllum[["Illuminant"]]==Rillum ),]
Rrgbwhitex<-unlist(Rrgbwhitergb[paste('X',observer,sep='')])
Rrgbwhitey<-unlist(Rrgbwhitergb[paste('Y',observer,sep='')])
Rrgbwhitez<-unlist(Rrgbwhitergb[paste('Z',observer,sep='')])
if (is.na(gamma)) gamma<-Rgamma
if (RGBModel=='sRGB') gamma<- -gamma
RGBmatrix<-apply(RGBmatrix,1:2,function(x) InvCompand(x, gamma))
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
xr<-unlist(Rrgb['xRed'])
yr<-unlist(Rrgb['yRed'])
xg<-unlist(Rrgb['xGreen'])
yg<-unlist(Rrgb['yGreen'])
xb<-unlist(Rrgb['xBlue'])
yb<-unlist(Rrgb['yBlue'])
m<-matrix(c(xr/yr,xg/yg,xb/yb,1.0,1.0,1.0,(1.0-xr-yr)/yr,(1.0-xg-yg)/yg,(1.0-xb-yb)/yb),3,3,byrow=TRUE)
mi<-solve(m)
sr<- mi %*% matrix(c(Rrgbwhitex,Rrgbwhitey,Rrgbwhitez),3,1)
MtxRGB2XYZ<- t(matrix(sr,3,3,byrow=TRUE) * m)
XYZ<-RGBmatrix %*% MtxRGB2XYZ
Adz<-CATmatrix %*% matrix(c(Rx,Ry,Rz),3,1,byrow=T)
Asz<-CATmatrix %*% matrix(c(Rrgbwhitex,Rrgbwhitey,Rrgbwhitez),3,1,byrow=T)
xyz2<- t(apply(XYZ,1, function(x) rowSums( matrix(x,3,3,byrow=T) * CATmatrix)))
xyz2<- xyz2 *  matrix(c(Adz/Asz),dim(XYZ)[1],3,byrow=T)
XYZ<-t(apply(xyz2,1, function(x) c(x) %*% t(solve(CATmatrix))))
XYZ
}

PreucilAngle<-function(RGBmatrix) {
# The Reproduction of Colour  By R. W. G. Hunt
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
apply(RGBmatrix,1,function(x) {
R<-x[1]
G<-x[2]
B<-x[3]
if ((R >= G) & (G >= B)) return(60*(G-B)/(R-B))
if ((G > R) & (R >= B)) return(60*(2-(R-B)/(G-B)))
if ((G >= B) & (B > R)) return(60*(2+(B-R)/(G-R)))
if ((B > G) & (G > R)) return(60*(4-(G-R)/(B-R)))
if ((B > R) & (R >= G)) return(60*(4+(R-G)/(B-G)))
if ((R >= B) & (B > G)) return(60*(6-(B-G)/(R-G)))
})
}

PreucilPercentGreyness<-function(RGBmatrix) { # The Reproduction of Colour  By R. W. G. Hunt pp. 513
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
apply(RGBmatrix,1,function(x) {
M<-max(x)
L<-min(x)
H<-x[which(x != c(M,L))]
100*L/H
})
}

PreucilPercentHueError<-function(RGBmatrix) { # The Reproduction of Colour  By R. W. G. Hunt
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
apply(RGBmatrix,1,function(x) {
M<-max(x)
L<-min(x)
H<-x[which(x != c(M,L))]
100*(M-L)/(H-L)
})
}

RGB2hue<-function(RGBmatrix) {
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
atan2(sqrt(3)*(RGBmatrix[,2]-RGBmatrix[,3]), 2*RGBmatrix[,1]-RGBmatrix[,2]-RGBmatrix[,3])
}

CIE1976chroma<-function(CIELMatrix){
# CIELab and CIELuv input
if (is.null(dim(CIELMatrix))) if (length(CIELMatrix)>2) CIELMatrix<-matrix(CIELMatrix, ncol=3,byrow=TRUE)
sqrt(CIELMatrix[,2]^2+CIELMatrix[,3]^2)
}
CIE1976hueangle<-function(CIELMatrix){
#  CIELab and CIELuv input
if (is.null(dim(CIELMatrix))) if (length(CIELMatrix)>2) CIELMatrix<-matrix(CIELMatrix, ncol=3,byrow=TRUE)
#atan(CIELMatrix[,3]/CIELMatrix[,2])*180/pi
atan2(CIELMatrix[,3],CIELMatrix[,2])*180/pi
}

CIE1976uvSaturation<-function(uvMatrix, whitepoint){
# CIELuv input
if (is.null(dim(uvMatrix))) if (length(uvMatrix)>2) uvMatrix<-matrix(uvMatrix, ncol=3,byrow=TRUE)
13*sqrt((uvMatrix[,2]-whitepoint[1])^2 + (uvMatrix[,3]-whitepoint[2])^2)
}

WestlandBlacknessIndex<-function(CIELabMatrix){
# (Westland, et al., 2006) blackness index
if (is.null(dim(CIELabMatrix))) if (length(CIELabMatrix)>2) CIELabMatrix<- matrix(CIELabMatrix, ncol=3,byrow=TRUE)
9.6523-0.3118*CIELabMatrix[,1]-0.0054*CIELabMatrix[,2]^2-0.0003*CIELabMatrix[,3]^2
}

RGB2LSLM<-function(RGBmatrix) {#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon
if (is.null(dim(RGBmatrix))) if (length(RGBmatrix)>2) RGBmatrix<-matrix(RGBmatrix, ncol=3,byrow=TRUE)
RGBmatrix<-RGBmatrix^2.2
cbind(0.209*(RGBmatrix[,1] - 0.5) + 0.715*(RGBmatrix[,2] - 0.5) + 0.076*(RGBmatrix[,3] - 0.5), 0.209*(RGBmatrix[,1] - 0.5) + 0.715*(RGBmatrix[,2] - 0.5) - 0.924*(RGBmatrix[,3] - 0.5),
3.14*(RGBmatrix[,1] - 0.5) - 2.799*(RGBmatrix[,2] - 0.5) - 0.349*(RGBmatrix[,3] - 0.5))
}

LSLM2RGB<-function(LSLMmatrix) {#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon
if (is.null(dim(LSLMmatrix))) if (length(LSLMmatrix)>2) LSLMmatrix<-matrix(LSLMmatrix, ncol=3,byrow=TRUE)
cbind(LSLMmatrix[,1] - 0.013*LSLMmatrix[,2] + 0.252*LSLMmatrix[,3] + 0.5, LSLMmatrix[,1] + 0.110*LSLMmatrix[,2] - 0.074*LSLMmatrix[,3] + 0.5, LSLMmatrix[,1] - LSLMmatrix[,2] + 0.5)
}

XYZ2LMS<-function(XYZmatrix)   XYZmatrix %*% (matrix(c(0.15514,0.54312,-0.03286,-0.15514,0.45684,0.03286,0,0,0.01608),3,3,TRUE))
#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon

LMS2XYZ<-function(LMSmatrix)   LMSmatrix %*% matrix(c(2.944813,-3.500978,13.17218,1.000040,1.000040,0.00000,0.000000,0.000000,62.18905),3,3,TRUE)
#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon

RGB2LMS<-function(RGBmatrix)
{#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon
if (any(RGBmatrix>1)) RGBmatrix<-RGBmatrix/255
RGBmatrix<-RGBmatrix^2.2
RGBmatrix %*% matrix(c(17.8824, 43.5161, 4.11935, 3.45565, 27.1554, 3.86714, 0.0299566, 0.184309, 1.46709),3,3,TRUE)
}

LMS2RGB<-function(LMSmatrix)
{#Digital Video Colourmaps for Checking the Legibility of Displays by Dichromats Francoise Vienot, Hans Brettel,John D. Mollon
l<-LMSmatrix %*% matrix(c(0.0809444479,-0.130504409,0.116721066,-0.0102485335,0.0540193266,-0.113614708,-0.000365296938,-0.00412161469,0.693511405),3,3,TRUE)
matrix(as.integer(l^(1/2.2)*255),ncol=3,byrow=F)
}

spectra2lux<-function(spectraIn=NA, ciexyzIn=NA,wlIn=NA, wlInterval=NA)
{# convert spectral data to illuminance
#spectraIn = spectra from spectrophotometer
#ciexyzIn = photopic response
#wlIn = range of output wavelengths
#wlInterval = arbitrary wavelength interval to be applied to all series through interpolation
#require(Hmisc)
if (any(is.na(spectraIn))) stop('<<spectraIn>> must be a numeric array nx2')
if (is.null(dim(spectraIn))) stop('<<spectraIn>> must be a numeric array nx2')
if (dim(spectraIn)[2] != 2) stop('<<spectraIn>> must be a numeric array nx2')
if (!is.numeric(spectraIn)) stop('<<spectraIn>> must be a numeric array nx2')
if (any(is.na(wlIn))) wlIn<-c(min(spectraIn[,1]), max(spectraIn[,1]))
if (any(is.na(wlInterval))) wlInterval<-spectraIn[2,1]-spectraIn[1,1]
if (any(is.na(ciexyzIn))) ciexyzIn<-get("ciexyz31", envir = environment())
wlMin<- min(wlIn)
wlMax<- max(wlIn)
wlSeq<-seq(wlMin,wlMax,wlInterval)
if ((min(spectraIn[,1])<=wlMin) & (max(spectraIn[,1])>=wlMax)) spectraIn <- spectraIn[which(spectraIn[,1] >=wlMin & spectraIn[,1] <=wlMax ),] else {
# extrapolate 
f <- approxfun(spectraIn[,1], spectraIn[,2])
xtrapolated<-approxExtrap(spectraIn[,1], spectraIn[,2], wlSeq)$y
xtrapolated[which(xtrapolated<0)]<-0
#curve(f(x), wlMin,wlMax, col = "green2")
#points(spectraIn)
spectraIn<-cbind(wlSeq,xtrapolated)
#curve(f(x), wlMin,wlMax, col = "green2",ylim=c(0,40))
#points(spectraIn)
}
if ((min(ciexyzIn[,1])<=wlMin) & (max(ciexyzIn[,1])>=wlMax)) ciexyzIn <- ciexyzIn[which(ciexyzIn[,1] %in% wlSeq),] else {
# extrapolate 
f <- approxfun(ciexyzIn[,1], ciexyzIn[,2])
xtrapolated1<-approxExtrap(ciexyzIn[,1], ciexyzIn[,2], wlSeq)$y
xtrapolated1[which(xtrapolated<0)]<-0
f <- approxfun(ciexyzIn[,1], ciexyzIn[,3])
xtrapolated2<-approxExtrap(ciexyzIn[,1], ciexyzIn[,3], wlSeq)$y
xtrapolated2[which(xtrapolated<0)]<-0
f <- approxfun(ciexyzIn[,1], ciexyzIn[,4])
xtrapolated3<-approxExtrap(ciexyzIn[,1], ciexyzIn[,4], wlSeq)$y
xtrapolated3[which(xtrapolated<0)]<-0
ciexyzIn<-cbind(wlSeq,xtrapolated1,xtrapolated2,xtrapolated3)
}
683 * sum( spectraIn[,2] * ciexyzIn[,3])  # multiply spectral distribution curve * photopic response curve
}

spectra2XYZ<-function(spectraIn=NA, illuminantIn=NA, ciexyzIn=NA,wlIn=NA, wlInterval=NA)
{# convert spectral data to tristimulus values
#spectraIn = spectra from spectrophotometer
#illuminantIn = illuminant
#ciexyzIn = photopic response
#wlIn = range of output wavelengths
#wlInterval = arbitrary wavelength interval to be applied to all series through interpolation
#require(Hmisc)
if (any(is.na(spectraIn))) stop('<<spectraIn>> must be a numeric array nx2')
if (is.null(dim(spectraIn))) stop('<<spectraIn>> must be a numeric array nx2')
if (dim(spectraIn)[2] != 2) stop('<<spectraIn>> must be a numeric array nx2')
if (!is.numeric(spectraIn)) stop('<<spectraIn>> must be a numeric array nx2')
if (any(is.na(wlIn))) wlIn<-c(min(spectraIn[,1]), max(spectraIn[,1]))
if (any(is.na(wlInterval))) wlInterval<-spectraIn[2,1]-spectraIn[1,1]
if (any(is.na(illuminantIn))) illuminantIn<-get("illuminantD65", envir = environment())[which(get("illuminantD65", envir = environment())[,1] %in% seq(min(get("illuminantD65", envir = environment())[,1]), max(get("illuminantD65", envir = environment())[,1]), wlInterval)),]
if (any(is.na(ciexyzIn))) ciexyzIn<-get("ciexyz31", envir = environment())
wlMin<- min(wlIn)
wlMax<- max(wlIn)
wlSeq<-seq(wlMin,wlMax,wlInterval)
if ((min(illuminantIn[,1])<=wlMin) & (max(illuminantIn[,1])>=wlMax)) illuminantIn <- illuminantIn[which(illuminantIn[,1] %in% wlSeq),] else {
# extrapolate 
f <- approxfun(illuminantIn[,1], illuminantIn[,2])
xtrapolated<-approxExtrap(illuminantIn[,1], illuminantIn[,2], wlSeq)$y
xtrapolated[which(xtrapolated<0)]<-0
illuminantIn<-cbind(wlSeq,xtrapolated)
}
if ((min(spectraIn[,1])<=wlMin) & (max(spectraIn[,1])>=wlMax)) spectraIn <- spectraIn[which(spectraIn[,1] %in% wlSeq),] else {
# extrapolate 
f <- approxfun(spectraIn[,1], spectraIn[,2])
xtrapolated<-approxExtrap(spectraIn[,1], spectraIn[,2], wlSeq)$y
xtrapolated[which(xtrapolated<0)]<-0
#curve(f(x), wlMin,wlMax, col = "green2")
#points(spectraIn)
spectraIn<-cbind(wlSeq,xtrapolated)
#curve(f(x), wlMin,wlMax, col = "green2",ylim=c(0,40))
#points(spectraIn)
}
if ((min(ciexyzIn[,1])<=wlMin) & (max(ciexyzIn[,1])>=wlMax)) ciexyzIn <- ciexyzIn[which(ciexyzIn[,1] %in% wlSeq),] else {
# extrapolate 
f <- approxfun(ciexyzIn[,1], ciexyzIn[,2])
xtrapolated1<-approxExtrap(ciexyzIn[,1], ciexyzIn[,2], wlSeq)$y
xtrapolated1[which(xtrapolated<0)]<-0
f <- approxfun(ciexyzIn[,1], ciexyzIn[,3])
xtrapolated2<-approxExtrap(ciexyzIn[,1], ciexyzIn[,3], wlSeq)$y
xtrapolated2[which(xtrapolated<0)]<-0
f <- approxfun(ciexyzIn[,1], ciexyzIn[,4])
xtrapolated3<-approxExtrap(ciexyzIn[,1], ciexyzIn[,4], wlSeq)$y
xtrapolated3[which(xtrapolated<0)]<-0
ciexyzIn<-cbind(wlSeq,xtrapolated1,xtrapolated2,xtrapolated3)
}
r<-illuminantIn[,2] * spectraIn[,2] # multiply illuminant by spectra
xv<-sum(r * ciexyzIn[,2])/100.0 # sum the product of illuminant, spectra and CIE XYZ Color Matching Functions
yv<-sum(r * ciexyzIn[,3])/100.0
zv<-sum(r * ciexyzIn[,4])/100.0
k<-sum(illuminantIn[,2] * ciexyzIn[,3]) # divide by the photopic response
xv<-xv/k
yv<-yv/k
zv<-zv/k
c(xv,yv,zv)
}

createIsoTempLinesTable <- function(SPD=NA,CIETable = get("ciexyz31", envir = environment()), TCS = get("TCSdata", envir = environment())){
# generate data for isotemperature lines needed for calculating correlated color temperature
# Light source SPD
# reference data values CIETable
# The CIE 1931 2 degree Standard Colorimetric Observer
# wavelength(nm)   xbar    ybar zbar
# The spectral reflectance data of 14 color test samples for CRI TCS
# wavelength (nm) TCS1 TCS2 TCS3 ... TCS14

  wavelength_spd <- SPD[,1]
  spd <- SPD[,2]

        wavelength <- CIETable[,1]
        xbar <- CIETable[,2]
        ybar <- CIETable[,3]
        zbar <- CIETable[,4]

# Data for isotemperature lines needed for calculating correlated color temperature

# The following provides a table of isotemperature lines for use with the Robertson Method
# (Robertson, 1968) to interpolate isotemperature lines from the CIE 1960 UCS.
# The spacing of the isotemp lines is very small (1 1/MK) so very little
# interpolation is actually needed for determining CCT. The latest (2002)
# recommended values for the physical constants determining blackbody
# radiation spectra are used

dwave = wavelength[2]-wavelength[1] # wavelength increment = 1 nm

ubar = (2/3)*xbar
vbar = ybar
wbar = -0.5*xbar + (3/2)*ybar + 0.5*zbar

# 2002 CODATA recommended values
h = 6.6260693e-34
c2002 = 299792458
k = 1.3806505e-23
c1 = 2*pi*h*c2002^2
c2 = h*c2002/k
MrecpK <- c(0.01, 1:600) # mega reciprocal Kelvin values of isotemperature lines
TisotempLines <- 1/(MrecpK*1e-6)


u=v=sl=m=vector(mode='integer',0)
for (i in 1:length(TisotempLines)){
   spdref = c1 * (1e-9*wavelength)^-5 / (exp(c2/(TisotempLines[i]* 1e-9*wavelength)) - 1)
   spdref = spdref/max(spdref)
   wave = wavelength*1e-9
   
   # Equations from Wyszecki and Sitles, Color Science, 2nd ed. 1982, page
   # 226 and 227
   U = sum(spdref*ubar)
   V = sum(spdref*vbar)
   W = sum(spdref*wbar)
   R = U+V+W
   u[i] = U/R
   v[i] = V/R
   
   Uprime = c1*c2*(TisotempLines[i])^-2*sum(wave^-6*ubar*exp(c2/(wave*TisotempLines[i]))*(exp(c2/(wave*(TisotempLines[i])))-1)^-2)*dwave
   Vprime = sum(c1*c2*TisotempLines[i]^-2*wave^-6*vbar*exp(c2/(wave*TisotempLines[i]))*(exp(c2/(wave*(TisotempLines[i])))-1)^-2)*dwave
   Wprime = sum(c1*c2*TisotempLines[i]^-2*wave^-6*wbar*exp(c2/(wave*TisotempLines[i]))*(exp(c2/(wave*(TisotempLines[i])))-1)^-2)*dwave
   Rprime = Uprime+Vprime+Wprime
   
   sl[i] = (Vprime*R-V*Rprime)/(Uprime*R-U*Rprime)
   m[i] = -1/sl[i]
}
cbind(T=TisotempLines, u=u, v=v, m=m) # isoTempLinesTable
}


spectra2CCT <- function(SPD=NA, isoTempLinesTable=NA,CIETable = get("ciexyz31", envir = environment()), TCS = get("TCSdata", envir = environment())){
# Correlated Color Temperature CCT
if(any(is.na(isoTempLinesTable))) isoTempLinesTable=createIsoTempLinesTable(SPD)
m <- isoTempLinesTable[,"m"]
TisotempLines <- isoTempLinesTable[,"T"]

wavelength_spd <- SPD[,1]
spd <- SPD[,2]
wavelength <- CIETable[,1]
xbar <- CIETable[,2]
ybar <- CIETable[,3]
zbar <- CIETable[,4]

# Correlated Color Temperature (CCT)
# Interpolate CIE functions to spd increments

xbar <- approx(wavelength,xbar,wavelength_spd)$y
xbar[which(is.na(xbar))] <- 0.0
ybar <- approx(wavelength,ybar,wavelength_spd)$y
ybar[which(is.na(ybar))] <- 0.0
zbar <- approx(wavelength,zbar,wavelength_spd)$y
zbar[which(is.na(zbar))] <- 0.0
# Calculate Chromaticity Coordinates
X <- trapz(wavelength_spd,spd*xbar)
Y <- trapz(wavelength_spd,spd*ybar)
Z <- trapz(wavelength_spd,spd*zbar)
x <- X/(X+Y+Z)
y <- Y/(X+Y+Z)
u <- 4*x/(-2*x+12*y+3)
v <- 6*y/(-2*x+12*y+3)
cat('x <- ',x,'\ty <- ',y,'\n')
ut <- isoTempLinesTable[,"u"]
vt <- isoTempLinesTable[,"v"]
tt <- m #isoTempLinesTable$m

# Find adjacent lines to (us, vs) 
n <- length (TisotempLines)
index <- 0
d1 <- ((v-vt[1]) - tt[1]*(u-ut[1]))/sqrt(1+tt[1]*tt[1])
for (i in 2:n){
    d2 <- ((v-vt[i]) - tt[i]*(u-ut[i]))/sqrt(1+tt[i]*tt[i])
    if (d1/d2 < 0){
        index <- i
        break
   } else d1 <- d2
}
if (index == 0) {
    Tc <- -1 # Not able to calculate CCT, u, v coordinates outside range.
    stop('Not able to calculate CCT, u, v coordinates outside range.')
} else {
    # Calculate CCT by interpolation between isotemperature lines
    Tc <- 1/(1/TisotempLines[index-1]+d1/(d1-d2)*(1/TisotempLines[index]-1/TisotempLines[index-1]))
    return(Tc)
}
}

spectra2CRIGAIFSCI <- function(SPD=NA, isoTempLinesTable=NA, CCT=NA, CIETable = get("ciexyz31", envir = environment()), TCS = get("TCSdata", envir = environment())){
# CRI, GAI and FSCI
# Color Rendering Index CRI
# Gamut Area Index GAI
# full spectrum index FSCI
if(any(is.na(isoTempLinesTable))) isoTempLinesTable=createIsoTempLinesTable(SPD)
if(any(is.na(CCT))) CCT=spectra2CCT(SPD, isoTempLinesTable)

wavelength_spd <- SPD[,1]
spd <- SPD[,2]
wavelength <- CIETable[,1]
xbar <- CIETable[,2]
ybar <- CIETable[,3]
zbar <- CIETable[,4]
# Interpolate CIE functions to spd increments

xbar <- approx(wavelength,xbar,wavelength_spd)$y
xbar[which(is.na(xbar))] <- 0.0
ybar <- approx(wavelength,ybar,wavelength_spd)$y
ybar[which(is.na(ybar))] <- 0.0
zbar <- approx(wavelength,zbar,wavelength_spd)$y
zbar[which(is.na(zbar))] <- 0.0

# calculate the Color Rendering Indices (CRI and its 14 indices)
# Calculate Reference Source Spectrum, spdref. 
if (CCT < 5000){
    c1 = 3.7418e-16;
    c2 = 1.4388e-2;
    spdref = c1 * (1e-9*wavelength_spd)^-5 / (exp(c2/(CCT* 1e-9*wavelength_spd)) - 1)
} else {
    if (CCT <= 25000){
        #load('CIEDaySn','wavelength','S0','S1','S2');
        S0 <- get("daylightcomponents", envir = environment())[["S0"]]
        S1 <- get("daylightcomponents", envir = environment())[["S1"]]
        S2 <- get("daylightcomponents", envir = environment())[["S2"]]
        if (CCT <= 7000){
            xd = -4.6070e9 / CCT^3 + 2.9678e6 / CCT^2 + 0.09911e3 / CCT + 0.244063
        } else {
            xd = -2.0064e9 / CCT^3 + 1.9018e6 / CCT^2 + 0.24748e3 / CCT + 0.237040
        }
        yd = -3.000*xd*xd + 2.870*xd - 0.275
        M1 = (-1.3515 - 1.7703*xd + 5.9114*yd) / (0.0241 + 0.2562*xd - 0.7341*yd)
        M2 = (0.0300 - 31.4424*xd + 30.0717*yd) / (0.0241 + 0.2562*xd - 0.7341*yd)
        spdref = S0 + M1*S1 + M2*S2
        spdref = approx(wavelength,spdref,wavelength_spd)
        spdref[which(is.na(spdref))] = 0.0
    } else {
        return(NA)
        }
    }

# Interpolate TCS values from 5 nm to spd nm increments
TCS_1 = matrix(0, length(wavelength_spd),14)
for (i in 1:14) TCS_1[,i] = approx(TCS[,1],TCS[,i+1],wavelength_spd,'linear',0)$y
TCS_1[which(is.na(TCS_1))] = 0

# Calculate u, v chromaticity coordinates of samples under test illuminant, uk, vk and
# reference illuminant, ur, vr.
Yki = vector(mode='integer',length=14)
Yri = vector(mode='integer',length=14)
uki = vector(mode='integer',length=14)
vki = vector(mode='integer',length=14)
uri = vector(mode='integer',length=14)
vri = vector(mode='integer',length=14)
X = trapz(wavelength_spd,spd * xbar)
Y = trapz(wavelength_spd,spd * ybar)
Z = trapz(wavelength_spd,spd * zbar)
Yknormal = 100 / Y
Yk = Y*Yknormal
uk = 4*X/(X+15*Y+3*Z)
vk = 6*Y/(X+15*Y+3*Z)
X = trapz(wavelength_spd,spdref * xbar)
Y = trapz(wavelength_spd,spdref * ybar)
Z = trapz(wavelength_spd,spdref * zbar)
Yrnormal = 100 / Y
Yr = Y*Yrnormal
ur = 4*X/(X+15*Y+3*Z)
vr = 6*Y/(X+15*Y+3*Z)
for (i in 1:14){
	X = trapz(wavelength_spd,spd * TCS_1[,i] * xbar)
	Y = trapz(wavelength_spd,spd * TCS_1[,i] * ybar)
	Z = trapz(wavelength_spd,spd * TCS_1[,i] * zbar)
	Yki[i] = Y*Yknormal
	uki[i] = 4*X/(X+15*Y+3*Z)
	vki[i] = 6*Y/(X+15*Y+3*Z)
	X = trapz(wavelength_spd,spdref * TCS_1[,i] * xbar)
	Y = trapz(wavelength_spd,spdref * TCS_1[,i] * ybar)
	Z = trapz(wavelength_spd,spdref * TCS_1[,i] * zbar)
	Yri[i] = Y*Yrnormal
	uri[i] = 4*X/(X+15*Y+3*Z)
	vri[i] = 6*Y/(X+15*Y+3*Z)
}

# Check tolerance for reference illuminant
DC = sqrt((uk-ur)^2 + (vk-vr)^2)

# Apply adaptive (perceived) color shift.
ck = (4 - uk - 10*vk) / vk
dk = (1.708*vk + 0.404 - 1.481*uk) / vk
cr = (4 - ur - 10*vr) / vr
dr = (1.708*vr + 0.404 - 1.481*ur) / vr

ukip = vector(mode='integer',length=14)
vkip  = vector(mode='integer',length=14)
for (i in 1:14){
	cki = (4 - uki[i] - 10*vki[i]) / vki[i]
	dki = (1.708*vki[i] + 0.404 - 1.481*uki[i]) / vki[i]
	ukip[i] = (10.872 + 0.404*cr/ck*cki - 4*dr/dk*dki) / (16.518 + 1.481*cr/ck*cki - dr/dk*dki);
	vkip[i] = 5.520 / (16.518 + 1.481*cr/ck*cki - dr/dk*dki);
}

#  Transformation into 1964 Uniform space coordinates.
Wstarr = vector(mode='integer',length=14)
Ustarr = vector(mode='integer',length=14)
Vstarr = vector(mode='integer',length=14)
Wstark = vector(mode='integer',length=14)
Ustark = vector(mode='integer',length=14)
Vstark = vector(mode='integer',length=14)
for (i in 1:14){
	Wstarr[i] = 25*Yri[i]^.333333 - 17
	Ustarr[i] = 13*Wstarr[i]*(uri[i] - ur)
	Vstarr[i] = 13*Wstarr[i]*(vri[i] - vr)
	
	Wstark[i] = 25*Yki[i]^.333333 - 17
	Ustark[i] = 13*Wstark[i]*(ukip[i] - ur) # after applying the adaptive color shift, u'k = ur
	Vstark[i] = 13*Wstark[i]*(vkip[i] - vr) # after applying the adaptive color shift, v'k = vr
}

# Determination of resultant color shift, delta E.
deltaE = vector(mode='integer',length=14)
R = vector(mode='integer',length=14)
for (i in 1:14){
	deltaE[i] = sqrt((Ustarr[i] - Ustark[i])^2 + (Vstarr[i] - Vstark[i])^2 + (Wstarr[i] - Wstark[i])^2)
	R[i] = 100 - 4.6*deltaE[i]
}
Ra = sum(R[1:8])/8

# fourth, calculate the gamut area formed by the 8 CIE standard color samples
ukii=c(uki[1:8],uki[1])
vkii=1.5*c(vki[1:8],vki[1])
Ga=polyarea(ukii,vkii)
# Normalize gamut area to equal energy source 
Ga=Ga/0.00728468*100
#fprintf(1,'Gamut Area Index = %.1f\n',Ga)

# Fifth, calculate the FSI (full spectrum index)
# Calculates the Full-spectrum Index

# Interpolate to wavelength interval of 1nm from 380nm to 730nm
numWave = 351
t9=380:730
wavelength_spd <- SPD[,1]
spd <- SPD[,2]
spd=interp1(wavelength_spd,spd,t9,method='spline')
spd[which(is.na(spd))] = 0.0
spd = spd/sum(spd) # Normalize the relative spd so that the total power equals 1
#Equal energy cumulative spd
EEcum=seq(1/numWave,1,by=1/numWave)
#Calculate FSI

sumSqrDiff <- vector(mode='integer',length=numWave)
for (j in 1:numWave ){
cum = cumsum(spd) # A MatLab function for cumulative sums 
sqrDiff = (cum-EEcum)^2
sumSqrDiff[j]=sum(sqrDiff)
spd=circshift(spd,1)
}
FSI=mean(sumSqrDiff)
FSCI=100-5.1*FSI

return(cbind(CRIra = Ra, GamutAreaIndex = Ga, FSCI = FSCI))
}

CIELabtoDIN99<-function(Lab){
# Conversion from CIELAB color space to DIN99 coordinates
# input: c(L,a,b)
# http://de.wikipedia.org/w/index.php?title=Diskussion:DIN99-Farbraum
L <- Lab[1]
a <- Lab[2]
b <- Lab[3]
kE <- 1.0 #brightness factor, 1.0 for CIE reference conditions
kCH <- 1.0 #chroma and hue factor, 1.0 for CIE reference conditions
ang <- 2*pi/360*26 #26 degree angle in radians
L99f <- 100/log(139/100.) #L99 scaling factor - method 1
#L99f <- 303.67 #original L99 factor
L99o <- L99f/kE*log(1+.0039*L) #Lightness correction kE
if ((a==0.0) && (b==0.0)) {
a99o <- 0.0
b99o <- 0.0
} else {
eo <- a*cos(ang) + b*sin(ang)       #a rotation
     fo <- 0.83*(b*cos(ang) - a*sin(ang)) #b rotation/stretching
     Go <- sqrt(eo^2 + fo^2)         #chroma
     C99o <- log(1.0 + 0.075 * Go) / (0.0435 * kCH*kE)  #factor for chroma compression and viewing conditions
    heofo <- atan2(fo,eo) #arctan in four quadrants
     h99o <- heofo+ang  #hue rotation
    a99o <- C99o * cos (h99o)
    b99o <- C99o * sin (h99o)
}
cbind(L99o,a99o,b99o)
}

DIN99toCIELab<-function(Lab99o){
# Conversion from DIN99 coordinates to CIELAB color space
# input: c(L99o,a99o,b99o)
# http://de.wikipedia.org/w/index.php?title=Diskussion:DIN99-Farbraum
L99o <- Lab99o[1]
a99o <- Lab99o[2]
b99o <- Lab99o[3]
kE <- 1.0 #brightness factor, 1.0 for CIE reference conditions
kCH <- 1.0 #chroma and hue factor, 1.0 for CIE reference conditions
ang <- 2*pi/360*26 #26 degrees in radians
L99f <- 100/log(139/100.0) #corrected L99 factor
#L99f <- 303.67 #original L99 factor
L <- (exp(L99o*kE/L99f)-1.0)/0.0039
h99ef <- atan2(b99o,a99o) #arctan in four quadrants
heofo <- h99ef-ang #backwards hue rotation
C99 <- sqrt(a99o^2+b99o^2)     #DIN99 chroma
G <- (exp(0.0435 * kE * kCH * C99)-1.0)/0.075 #factor for chroma decompression and viewing conditions
e <- G*cos(heofo)
f <- G*sin(heofo)
a <- (e * cos(ang) - (f/0.83) * sin(ang)) #rotation by 26 degrees
b <- (e * sin(ang) + (f/0.83) * cos(ang)) #rotation by 26 degrees
cbind(L=L,a=a,b=b)
}

Y2MunsellVtable1D1535<-function(Y){
# CIE XYZ "Y" to Munsell "V"
# NLSQ regression for obtaining similar results to table 1 from ASTM Standard D1535-08
Y*100-9.876522e-06 *(Y-1.749992e+02)^2 + 6.139636e-03*(Y-1.749992e+02) + 1.370407e+00
}

MunsellV2Y<-function(V) 1.1914*V-0.22533*V ^ 2 +0.23352*V ^ 3-0.020484*V ^ 4 +0.00081939*V ^ 5
# Munsell "V" to CIE XYZ "Y" 
# ASTM Standard D1535-08

Y2MunsellV<-function(Y){
# CIE XYZ "Y" to Munsell "V" 
# ASTM Standard D1535-08
a = 2.49268
b = 1.5614
C = 0.985
d = 0.1073
e = 3.084
f = 7.54
g = 0.0133
H = 2.3
j = 0.0084
k = 4.1
m = 0.0221
n = 0.39
P = 0.037
Q = 0.44
S = 1.28
T = 0.53
U = 0.87455
W = 0.9967
if (Y <= 0.9) V <- U * Y ^ W else V <- a * (Y ^ (1 / 3)) - b - C / (((d * Y - e) ^ 2) + f) + g / (Y ^ H)  + j *
sin((k * (Y ^ (1 / 3)) + 1)/ 180 * pi) + m / Y * sin((n * (Y - 2))/ 180 * pi) + P / Q / Y * sin((S * (Y - T))/ 180 * pi)
V
}

XYZ2Lab <- function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(XYZmatrix)) XYZmatrix<-matrix(XYZmatrix,ncol=3,byrow=TRUE)
kE <- 216.0 / 24389.0
kK <- 24389.0 / 27.0
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
xr <- XYZmatrix[,1] / Rx #* 100
yr <- XYZmatrix[,2] / Ry #* 100
zr <- XYZmatrix[,3] / Rz #* 100
fx <- ifelse((xr > kE) , (xr^(1.0 / 3.0)) , ((kK * xr + 16.0) / 116.0))
fy <- ifelse((yr > kE) , (yr^(1.0 / 3.0)) , ((kK * yr + 16.0) / 116.0))
fz <- ifelse((zr > kE) , (zr^(1.0 / 3.0)) , ((kK * zr + 16.0) / 116.0))
L <- 116.0 * fy - 16.0
a <- 500.0 * (fx - fy)
b <- 200.0 * (fy - fz)
cbind(L=L,a=a,b=b)
}

XYZ2Luv <- function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(XYZmatrix)) XYZmatrix<-matrix(XYZmatrix,ncol=3,byrow=TRUE)
kE <- 216.0 / 24389.0
kK <- 24389.0 / 27.0
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])

U = ( 4 * XYZmatrix[,1] ) / ( XYZmatrix[,1] + ( 15 * XYZmatrix[,2] ) + ( 3 * XYZmatrix[,3] ) )
V = ( 9 * XYZmatrix[,2] ) / ( XYZmatrix[,1] + ( 15 * XYZmatrix[,2] ) + ( 3 * XYZmatrix[,3] ) )
vY = XYZmatrix[,2] #/ 100
vY = ifelse (( vY > kE ), vY ^ ( 1/3 ), ( kK * vY  + 16) / 116 )
Ru = ( 4 * Rx ) / ( Rx + ( 15 * Ry ) + ( 3 * Rz ) )
Rv = ( 9 * Ry ) / ( Rx + ( 15 * Ry ) + ( 3 * Rz ) )
L = ( 116 * vY ) - 16
u = 13 * L* ( U - Ru )
v = 13 * L * ( V - Rv )

cbind(L=L,u=u,v=v)
}

xyY2XYZ<-function(xyYmatrix)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(xyYmatrix)) xyYmatrix<-matrix(xyYmatrix,ncol=3,byrow=TRUE)
if (dim(xyYmatrix)[1]>1) {
z<-apply(xyYmatrix,1,function(r) if (r[1] < 0.000001) return(c(X=0,Y=0,Z=0)) else return(c(X=(r[1] * r[3]) / r[2],Y=r[3],Z= ((1.0 - r[1] - r[2]) * r[3]) / r[2])) )
XYZ<-t(z)
} else if (xyYmatrix[3] < 0.000001) XYZ<-c(X=0,Y=0,Z=0) else XYZ<-cbind(X=(xyYmatrix[1] * xyYmatrix[3]) / xyYmatrix[2],Y=xyYmatrix[3],Z= ((1.0 - xyYmatrix[1] - xyYmatrix[2]) * xyYmatrix[3]) / xyYmatrix[2])
colnames(XYZ)<-c('X','Y','Z')
XYZ
}

XYZ2xyY<-function(XYZmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(XYZmatrix)) XYZmatrix<-matrix(XYZmatrix,ncol=3,byrow=TRUE)
Den <- rowSums(XYZmatrix)
DenG0<-which(Den > 0.0)
xyYmatrix<-XYZmatrix
xyYmatrix[DenG0,1:2]<- XYZmatrix[DenG0,1:2]/Den
xyYmatrix[DenG0,3]<- XYZmatrix[DenG0,2]
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
x <- Rx / (Rx + Ry + Rz)
y <- Ry / (Rx + Ry + Rz)
xyYmatrix[-DenG0,1]<-x
xyYmatrix[-DenG0,2]<-y
xyYmatrix
}

Lab2XYZ<-function(Labmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(Labmatrix)) Labmatrix<-matrix(Labmatrix,ncol=3,byrow=TRUE)
L<-Labmatrix[,1]
a<-Labmatrix[,2]
b<-Labmatrix[,3]
kE <- 216.0 / 24389.0
kK <- 24389.0 / 27.0
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
fy <- (L + 16.0) / 116.0
fx <- 0.002 * a + fy
fz <- fy - 0.005 * b
fx3 <- fx * fx * fx
fz3 <- fz * fz * fz
xr <- ifelse((fx3 > kE) , fx3 , ((116.0 * fx - 16.0) / kK))
yr <- ifelse((L > 8.0) , (((L + 16.0) / 116.0)^3.0) , (L / kK))
zr <- ifelse((fz3 > kE) , fz3 , ((116.0 * fz - 16.0) / kK))
X <- xr * Rx
Y <- yr * Ry
Z <- zr * Rz
cbind(X=X,Y=Y,Z=Z)
}

Luv2XYZ<-function(Luvmatrix,illuminant='D65',observer=2,RefWhite=get("XYZperfectreflectingdiffuser", envir = environment()))
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
if (!is.matrix(Luvmatrix)) Luvmatrix<-matrix(Luvmatrix,ncol=3,byrow=TRUE)
L<-Luvmatrix[,1]
u<-Luvmatrix[,2]
v<-Luvmatrix[,3]
kE <- 216.0 / 24389.0
kK <- 24389.0 / 27.0
R<-RefWhite[which(RefWhite[["Illuminant"]]==illuminant ),]
Rx<-unlist(R[paste('X',observer,sep='')])
Ry<-unlist(R[paste('Y',observer,sep='')])
Rz<-unlist(R[paste('Z',observer,sep='')])
Y <- ifelse((L > 8.0) , (((L + 16.0) / 116.0)^3.0) , (L / kK))
u0 <- (4.0 * Rx) / (Rx + 15.0 * Ry + 3.0 * Rz)
v0 <- (9.0 * Ry) / (Rx + 15.0 * Ry + 3.0 * Rz)
a <- (((52.0 * L) / (u + 13.0 * L * u0)) - 1.0) / 3.0
b <- -5.0 * Y
c <- -1.0 / 3.0
d <- Y * (((39.0 * L) / (v + 13.0 * L * v0)) - 5.0)
X <- (d - b) / (a - c)
Z <- X * a + b
cbind(X=X,Y=Y,Z=Z)
}

Compand<-function(linearV, gammaV)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html

if (gammaV > 0.0) companded <- ifelse((linearV >= 0.0) , (linearV^(1.0 / gammaV)) , -(-linearV^(1.0 / gammaV))) else if (gammaV < 0.0)
{
# sRGB
signV <- 1.0
if (linearV < 0.0)
{
signV <- -1.0
linearV <- -linearV
}
companded <- ifelse((linearV <= 0.0031308), (linearV * 12.92) , (1.055 * (linearV^(1.0 / 2.4)) - 0.055)) 
companded <- companded *signV
}
else
{
# L*
signV <- 1.0
if (linearV < 0.0)
{
signV <- -1.0
linearV <- -linearV
}
companded <- ifelse((linearV <= (216.0 / 24389.0)), (linearV * 24389.0 / 2700.0) , (1.16 * (linearV^(1.0 / 3.0)) - 0.16)) 
companded <- companded *signV
}
return(companded)
}

InvCompand<-function(companded, gammaV)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html

if (gammaV > 0.0) linearV <- ifelse((companded >= 0.0), (companded^gammaV) , -(-companded^gammaV)) else if (gammaV < 0.0) 
{
# sRGB
signV <- 1.0
if (companded < 0.0)
{
signV <- -1.0
companded <- -companded
}
linearV <- ifelse((companded <= 0.04045), (companded / 12.92) , (((companded + 0.055) / 1.055)^2.4)) 
linearV <- linearV *signV
}
else
{
# L*
signV <- 1.0
if (companded < 0.0)
{
signV <- -1.0
companded <- -companded
}
linearV <- ifelse((companded <= 0.08), (2700.0 * companded / 24389.0) , ((((1000000.0 * companded + 480000.0) * companded + 76800.0) * companded + 4096.0) / 1560896.0)) 
linearV <- linearV *signV
}
return(linearV)
}

compuphaseDifferenceRGB<-function(RGB1, RGB2)
{ # difference metric for RGB
# http://www.compuphase.com/cmetric.htm
  rmean <- mean( RGB1[1] + RGB2[1] )
  r <- RGB1[1] - RGB2[1]
  g <- RGB1[2] - RGB2[2]
  b <- RGB1[3] - RGB2[3]
  sqrt((((512+rmean)*r^2)/256) + 4*g^2 + (((767-rmean)*b^2)/256))
}

deltaE1976 <- function(Lab1, Lab2)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html

delL <- Lab1[1] - Lab2[1]
dela <- Lab1[2] - Lab2[2]
delb <- Lab1[3] - Lab2[3]
sqrt(delL * delL + dela * dela + delb * delb)
}

deltaE1994 <- function (Lab1, Lab2, textiles=FALSE)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
k1 <- ifelse((textiles == TRUE),0.048, 0.045)
k2 <- ifelse((textiles == TRUE),0.014, 0.015)
kL <- ifelse((textiles == TRUE),2.0, 1.0)
kC <- 1.0
kH <- 1.0
C1 <- sqrt(Lab1[2] * Lab1[2] + Lab1[3] * Lab1[3])
C2 <- sqrt(Lab2[2] * Lab2[2] + Lab2[3] * Lab2[3])
delA <- Lab1[2] - Lab2[2]
delB <- Lab1[3] - Lab2[3]
delC <- C1 - C2
delH2 <- delA * delA + delB * delB - delC * delC
delH <- ifelse((delH2 > 0.0),sqrt(delH2), 0.0) 
delL <- Lab1[1] - Lab2[1]
sL <- 1.0
sC <- 1.0 + k1 * C1
sH <- 1.0 + k2 * C1
vL <- delL / (kL * sL)
vC <- delC / (kC * sC)
vH <- delH / (kH * sH)
if (textiles == TRUE) return(sqrt(vL * vL + vC * vC + vH * vH)) else return(sqrt(vL * vL + vC * vC + vH * vH))
}

deltaE2000 <- function(Lab1, Lab2)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
kL <- 1.0
kC <- 1.0
kH <- 1.0
lBarPrime <- 0.5 * (Lab1[1] + Lab2[1])
c1 <- sqrt(Lab1[2] * Lab1[2] + Lab1[3] * Lab1[3])
c2 <- sqrt(Lab2[2] * Lab2[2] + Lab2[3] * Lab2[3])
cBar <- 0.5 * (c1 + c2)
cBar7 <- cBar^7
g <- 0.5 * (1.0 - sqrt(cBar7 / (cBar7 + 6103515625.0)))# 6103515625 = 25^7
a1Prime <- Lab1[2] * (1.0 + g)
a2Prime <- Lab2[2] * (1.0 + g)
c1Prime <- sqrt(a1Prime * a1Prime + Lab1[3] * Lab1[3])
c2Prime <- sqrt(a2Prime * a2Prime + Lab2[3] * Lab2[3])
cBarPrime <- 0.5 * (c1Prime + c2Prime)
h1Prime <- (atan2(Lab1[3], a1Prime) * 180.0) / pi
if (h1Prime < 0.0)
h1Prime <- h1Prime +360.0
h2Prime <- (atan2(Lab2[3], a2Prime) * 180.0) / pi
if (h2Prime < 0.0)
h2Prime <- h2Prime +360.0
hBarPrime <- ifelse((abs(h1Prime - h2Prime) > 180.0),(0.5 * (h1Prime + h2Prime + 360.0)), (0.5 * (h1Prime + h2Prime))) 
t <- 1.0 -
0.17 * cos(pi * (      hBarPrime - 30.0) / 180.0) +
0.24 * cos(pi * (2.0 * hBarPrime       ) / 180.0) +
0.32 * cos(pi * (3.0 * hBarPrime +  6.0) / 180.0) -
0.20 * cos(pi * (4.0 * hBarPrime - 63.0) / 180.0)
if (abs(h2Prime - h1Prime) <= 180.0) dhPrime <- h2Prime - h1Prime else dhPrime <- ifelse((h2Prime <= h1Prime),(h2Prime - h1Prime + 360.0), (h2Prime - h1Prime - 360.0)) 
dLPrime <- Lab2[1] - Lab1[1]
dCPrime <- c2Prime - c1Prime
dHPrime <- 2.0 * sqrt(c1Prime * c2Prime) * sin(pi * (0.5 * dhPrime) / 180.0)
sL <- 1.0 + ((0.015 * (lBarPrime - 50.0) * (lBarPrime - 50.0)) / sqrt(20.0 + (lBarPrime - 50.0) * (lBarPrime - 50.0)))
sC <- 1.0 + 0.045 * cBarPrime
sH <- 1.0 + 0.015 * cBarPrime * t
dTheta <- 30.0 * exp(-((hBarPrime - 275.0) / 25.0) * ((hBarPrime - 275.0) / 25.0))
cBarPrime7 <- cBarPrime * cBarPrime * cBarPrime * cBarPrime * cBarPrime * cBarPrime * cBarPrime
rC <- sqrt(cBarPrime7 / (cBarPrime7 + 6103515625.0))
rT <- -2.0 * rC * sin(pi * (2.0 * dTheta) / 180.0)
sqrt(
(dLPrime / (kL * sL)) * (dLPrime / (kL * sL)) +
(dCPrime / (kC * sC)) * (dCPrime / (kC * sC)) +
(dHPrime / (kH * sH)) * (dHPrime / (kH * sH)) +
(dCPrime / (kC * sC)) * (dHPrime / (kH * sH)) * rT)
}

deltaECMC <- function(Lab1, Lab2, L=1, C=1)
{# Based on: Bruce Justin Lindbloom 2013 http:#www.brucelindbloom.com/index.html?ColorCalculator.html
c1 <- sqrt(Lab1[2] * Lab1[2] + Lab1[3] * Lab1[3])
c2 <- sqrt(Lab2[2] * Lab2[2] + Lab2[3] * Lab2[3])
sl <- ifelse((Lab1[1] < 16.0),(0.511) ,((0.040975 * Lab1[1]) / (1.0 + 0.01765 * Lab1[1]))) 
sc <- (0.0638 * c1) / (1.0 + 0.0131 * c1) + 0.638
h1 <- ifelse((c1 < 0.000001),0.0 ,((atan2(Lab1[3], Lab1[2]) * 180.0) / pi)) 
while (h1 < 0.0)
h1 <- h1 +360.0
while (h1 >= 360.0)
h1 <- h1 -360.0
t <- ifelse((h1 >= 164.0) && (h1 <= 345.0),(0.56 + abs(0.2 * cos((pi * (h1 + 168.0)) / 180.0))), (0.36 + abs(0.4 * cos((pi * (h1 + 35.0)) / 180.0)))) 
c4 <- c1^4
f <- sqrt(c4 / (c4 + 1900.0))
sh <- sc * (f * t + 1.0 - f)
delL <- Lab1[1] - Lab2[1]
delC <- c1 - c2
delA <- Lab1[2] - Lab2[2]
delB <- Lab1[3] - Lab2[3]
dH2 <- delA * delA + delB * delB - delC * delC
v1 <- delL / (L * sl)
v2 <- delC / (C * sc)
v3 <- sh
sqrt(v1 * v1 + v2 * v2 + (dH2 / (v3 * v3)))
}

MunsellV2relativeLuminanceY<-function(V) 1.2219*V - 0.23111*V^2 + 0.23951*V^3 - 0.021009*V^4 + 0.0008404*V^5
# convert Munsell value V to relative luminance Y
# Color Appearance Models, Mark D. Fairchild pp. 96

CIEluminanceY2NCSblackness<-function(Y) 100 - 156*Y / (56 + Y )
# approximated NCS blackness s by the CIE luminance factor Y
# Introduction to Color Imaging Science, Lee pp. 366

huedegree<-function(MunIn){
# convert Munsell hue to degree
# based on Takahiro Onodera 2010 Color-Model-Munsell-Util http://annocpan.org/dist/Color-Model-Munsell-Util
# 10.0RP will return 0
# sapply(MunIn[,'H']
m<-sapply(MunIn , function(x) {if (any(is.na(x))) return(NA) else {
sidN<-gregexpr('^(\\d{1,2}\\.\\d{1,2}|\\d{1,2})(RP|YR|Y|GY|G|BG|B|PB|P|R)', x, perl=TRUE)
h1<-substr(x,attr(sidN[[1]], "capture.start")[1,][1],attr(sidN[[1]], "capture.start")[1,][1]+attr(sidN[[1]], "capture.length")[1,][1]-1)
hM<-substr(x,attr(sidN[[1]], "capture.start")[1,][2],attr(sidN[[1]], "capture.start")[1,][2]+attr(sidN[[1]], "capture.length")[1,][2]-1)
h2<-which(hM == get("MunsellHues", envir = environment()))
tmp <- (h2-1)*10+as.numeric(h1)
if (tmp>=100) tmp<-tmp %% 100
return(as.numeric(tmp))
} })
m
}

Adjust<-function(myColor, Factor, Gamma,IntensityMax=1.0)# HIDDEN
{# part of heuristic.wlnm2RGB
if (length(myColor)==1) {
if (myColor == 0.0) return(