---
title: 'RGB color correction for color blindess: protanopia, deuteranopia, tritanopia'
author: 'Martin Krzywinski'
web: http://mkweb.bcgsc.ca/colorblind
---

```{r}
gamma = 2.4

###############################################
# Linear RGB to XYZ
# https://en.wikipedia.org/wiki/SRGB
XYZ = matrix(c(0.4124564, 0.3575761, 0.1804375,
               0.2126729, 0.7151522, 0.0721750,
               0.0193339, 0.1191920, 0.9503041),
               byrow=TRUE,nrow=3)

SA = matrix(c(0.2126,0.7152,0.0722,
              0.2126,0.7152,0.0722,
              0.2126,0.7152,0.0722),byrow=TRUE,nrow=3)

###############################################
# XYZ to LMS, normalized to D65
# https://en.wikipedia.org/wiki/LMS_color_space
# Hunt, Normalized to D65
LMSD65  = matrix(c( 0.4002, 0.7076, -0.0808,
                   -0.2263, 1.1653,  0.0457,
                    0     , 0     ,  0.9182),
                   byrow=TRUE,nrow=3)
# Hunt, equal-energy illuminants
LMSEQ = matrix(c( 0.38971, 0.68898,-0.07868,
                 -0.22981, 1.18340, 0.04641,
                  0      , 0      , 1      ),
                byrow=TRUE,nrow=3)
# CIECAM97
SMSCAM97 = matrix(c(  0.8951,  0.2664, -0.1614,
                     -0.7502,  1.7135,  0.0367,
                      0.0389, -0.0685,  1.0296),
                  byrow=TRUE,nrow=3)
# CIECAM02
LMSCAM02 = matrix(c( 0.7328, 0.4296, -0.1624,
                    -0.7036, 1.6975,  0.0061,
                     0.0030, 0.0136,  0.9834),
                  byrow=TRUE,nrow=3)

###############################################
# Determine the color blindness correction in LMS space
# under the condition that the correction does not
# alter the appearance of white as well as 
# blue (for protanopia/deuteranopia) or red (for tritanopia).
# For achromatopsia, greyscale conversion is applied
# to the linear RGB values.
getcorrection = function(LMS,type="p",g=gamma) {
  red   = matrix(c(255,0,0),nrow=3)
  blue  = matrix(c(0,0,255),nrow=3)
  white = matrix(c(255,255,255),nrow=3)
  LMSr  = LMS %*% XYZ %*% apply(red,1:2,linearize,g)
  LMSb  = LMS %*% XYZ %*% apply(blue,1:2,linearize,g)
  LMSw  = LMS %*% XYZ %*% apply(white,1:2,linearize,g)
  if(type == "p") {
    x  = matrix(c(LMSb[2,1],LMSb[3,1],
                  LMSw[2,1],LMSw[3,1]),byrow=T,nrow=2)
    y = matrix(c(LMSb[1,1],LMSw[1,1]),nrow=2)
    ab = solve(x) %*% y
    C = matrix(c(0,ab[1,1],ab[2,1],0,1,0,0,0,1),byrow=T,nrow=3)
  } else if (type == "d") {
    x  = matrix(c(LMSb[1,1],LMSb[3,1],
                  LMSw[1,1],LMSw[3,1]),byrow=T,nrow=2)
    y = matrix(c(LMSb[2,1],LMSw[2,1]),nrow=2)
    ab = solve(x) %*% y
    C = matrix(c(1,0,0,ab[1,1],0,ab[2,1],0,0,1),byrow=T,nrow=3)
  } else if (type == "t") {
    x  = matrix(c(LMSr[1,1],LMSr[2,1],
                  LMSw[1,1],LMSw[2,1]),byrow=T,nrow=2)
    y = matrix(c(LMSr[3,1],LMSw[3,1]),nrow=2)
    ab = solve(x) %*% y
    C = matrix(c(1,0,0,0,1,0,ab[1,1],ab[2,1],0),byrow=T,nrow=3)
  } else if (type == "a" | type == "g") {
    C = matrix(c(0.2126,0.7152,0.0722,
                 0.2126,0.7152,0.0722,
                 0.2126,0.7152,0.0722),byrow=TRUE,nrow=3)
  }
  return(C)
}

# rgb is a column vector
convertcolor = function(rgb,LMS=LMSD65,type="d",g=gamma) {
  C = getcorrection(LMS,type)
  if(type == "a" | type == "g") {
    T = SA
  } else {
    M    = LMS %*% XYZ
    Minv = solve(M)
    T    = Minv %*% C %*% M
  }
  print(T)
  rgb_converted = T %*% apply(rgb,1:2,linearize,g)
  return(apply(rgb_converted,1:2,delinearize,g))
}

# This function implements the method by Vienot, Brettel, Mollon 1999.
# The approach is the same, just the values are different.
# http://vision.psychol.cam.ac.uk/jdmollon/papers/colourmaps.pdf
convertcolor2 = function(rgb,type="d",g=2.2) {
  xyz = matrix(c(40.9568, 35.5041, 17.9167,
                 21.3389, 70.6743, 7.98680,
                 1.86297, 11.4620, 91.2367),byrow=T,nrow=3)
  lms = matrix(c(0.15514, 0.54312, -0.03286,
                 -0.15514, 0.45684,0.03286,
                 0,0,0.01608),byrow=T,nrow=3)
  rgb = (rgb/255)**g
  if(type=="p") {
    S = matrix(c(0,2.02344,-2.52581,0,1,0,0,0,1),byrow=T,nrow=3)
    rgb = 0.992052*rgb+0.003974
  } else if(type=="d") {
    S = matrix(c(1,0,0,0.494207,0,1.24827,0,0,1),byrow=T,nrow=3)
    rgb = 0.957237*rgb+0.0213814
  } else {
    stop("Only type p,d defined for this function.")
  }
  M   = lms %*% xyz
  T   = solve(M) %*% S %*% M
  print(T)
  rgb = T %*% rgb
  rgb = 255*rgb**(1/g)
  return(rgb)
}

###############################################
# RGB to Lab
rgb2lab = function(rgb,g=gamma) {
  rgb = apply(rgb,1:2,linearize,g)
  xyz = XYZ %*% rgb
  delta = 6/29
  xyz = xyz / (c(95.0489,100,108.8840)/100)
  f = function(t) {
    if(t > delta**3) {
      return(t**(1/3))
    } else {
      return (t/(3*delta**2) + 4/29)
    }
  }
  L = 116*f(xyz[2]) - 16
  a = 500*(f(xyz[1]) - f(xyz[2]))
  b = 200*(f(xyz[2]) - f(xyz[3]))
  return(matrix(c(L,a,b),nrow=3))
}

# CIE76 (https://en.wikipedia.org/wiki/Color_difference)
deltaE = function(rgb1,rgb2) {
  lab1 = rgb2lab(rgb1)
  lab2 = rgb2lab(rgb2)
  return(sqrt(sum((lab1-lab2)**2)))
}

clip = function(v) {
  return(max(min(v,1),0))
}

###############################################
# RGB to/from linear RGB
#https://en.wikipedia.org/wiki/SRGB
linearize = function(v,g=gamma) {
  if(v <= 0.04045) {
    return(v/255/12.92)
  } else {
    return(((v/255 + 0.055)/1.055)**g)
  }
}

delinearize = function(v,g=gamma) {
  if(v <= 0.003130805) {
    return(255*12.92*clip(v))
  } else {
    return(255*clip(1.055*(clip(v)**(1/g))-0.055))
  }
}
pretty = function(x) {
  noquote(formatC(x,digits=10,format="f",width=9))
}

# a dark red
rgb1 = matrix(c(0,209,253),nrow=3)
# dark green
rgb2 = matrix(c(60,135,0),nrow=3)
# simulate deuteranopia
convertcolor(rgb1,type="d")
convertcolor(rgb2,type="d")
# get color distance before and after simulation
deltaE(rgb1,rgb2)
deltaE(convertcolor(rgb1,type="d"),convertcolor(rgb2,type="d"))
# transformation matrices for each color blindness type
M = LMSD65 %*% XYZ
pretty(solve(M) %*% getcorrection(LMSD65,"p") %*% M)
pretty(solve(M) %*% getcorrection(LMSD65,"d") %*% M)
pretty(solve(M) %*% getcorrection(LMSD65,"t") %*% M)
pretty(SA)
# method by Vienot, Brettel, Mollon, 1999
convertcolor2(rgb1,type="d",g=2.2)
convertcolor2(rgb2,type="d",g=2.2)
```

