--- 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) ```