sRGBに対応した色弱色差シミュレータを作る（車輪の再発明）

**偉大な先人達のおかげで、sRGBでの色弱シミュレータはブラウザのアドオンやスマホアプリなどで配布されています。**
**本記事で取り扱う変換方式や公式はそれらのごく一部であり、誤解を招く記載もありますが、「入門編」としてご理解ください。**

# [同じような](https://qiita.com/Kahiro-M/items/5d4cfeadc91073a8d278)の作ってなかった？
以前、「[色弱色差シミュレータ](https://qiita.com/Kahiro-M/items/5d4cfeadc91073a8d278)」を作ったが、対応している入力値がXYZ値（正確にはxy色度座標と輝度値）だった。

発光ダイオードの場合は製品の仕様値から、反射しているカラーチップなどは測色機による測定結果によって調べることができる。
しかし、ディスプレイは一般的に「現在表示している色がどんな光学的な値を出すか」は仕様などで公開されていない。
なので、sRGBに対応しているディスプレイに限定されるが、「RGB値を入力とした色弱色差シミュレータ」を作る。

# レシピ
### 材料
* sRGB->XYZの変換式
* XYZ->錐体反応値の変換式
* C型->P型,D型の変換式（いわゆる正常色覚から色弱への変換）
* Lab->色差の算出式

### 手順
1. RGB -> XYZ
2. XYZ -> LMS
3. LMS -> L'M'S'
4. L'M'S' -> X'Y'Z'
5. X'Y'Z' -> Lab
6. Lab -> 色差
7. 行列計算ライブラリ

# 作っていく
### まずは、何はなくともXYZにしよう
ディスプレイに表示するためのRGB値とは、ディスプレイが持っている三種類の光学素子を発光させる割合を示している。
Rは赤色発光させる素子の明るさを0〜100%の間で255段階に分けた値に該当する。同様にGは緑色、Bは青色の発光度合いを0〜255で表す。

各素子が発光し、複数の色を持った光線が同時に目にはいることで混色される。
この際、混ざる光の量が増えると「明るくなった」と感じる。「増える（足す）と明るくなる混色＝加法混色」という。

混色された発光色は「XYZ色空間」に変換することができる。ただし、RGB値をそのまま使うことはできず、割合で表現する必要がある。
なので、まずはRGB値を0~1に正規化する。

```math
R_{linear} = R / 255\\
G_{linear} = G / 255\\
B_{linear} = B / 255  
```

ディスプレイで表示されているRGBは必ずしも「入力した通りの出力（明るさ）になっている」とは限らない。
一般的には「低い入力値はより小さい出力値、大きい入力値はより大きな出力値」という物理的な特性があり、その補正[^1]をしなければならない。
[^1]: これをガンマ補正という

この補正には条件があり、極端に入力値が低い場合とそうではない場合で、場合分けの計算が必要となる。[sRGB - Wiki(英語)](https://en.wikipedia.org/wiki/SRGB#The_reverse_transformation)

```math
f(rgb) = \left\{
\begin{array}{ll}
\frac{rgb}{12.92} & (rgb \leq 0.04045) \\
\frac{rgb + 0.055}{1.055} & (otherwise)
\end{array}
\right.
```

これでようやくXYZと同じ単位系に揃うので、変換計算ができる。
変換式は連立方程式を解くのと同じ操作なので、行列で表現できる。（行列の方が逆の処理をする際に簡単に記述できるので、ここでは行列で進める。）

```math
\left[
 \begin{array}{r}
  X_{D65} \\
  Y_{D65} \\
  Z_{D65}
 \end{array}
\right]
=
\left[
 \begin{array}{rrr}
  0.4124 & 0.3576 & 0.1805\\
  0.2126 & 0.7152 & 0.0722\\
  0.0193 & 0.1192 & 0.9505\\
 \end{array}
\right]
\left[
 \begin{array}{r}
  R_{linear} \\
  G_{linear} \\
  B_{linear}
 \end{array}
\right]
```


JavaScript(math.js利用)で記載するなら、

``` .js
culcXYZfromRGBData(inRGBData){
	//https://en.wikipedia.org/wiki/SRGB#The_reverse_transformation
	var r = inRGBData.R / 255;
	var g = inRGBData.G / 255;
	var b = inRGBData.B / 255;

	r = r > 0.04045 ? Math.pow(((r + 0.055) / 1.055), 2.4) : (r / 12.92);
	g = g > 0.04045 ? Math.pow(((g + 0.055) / 1.055), 2.4) : (g / 12.92);
	b = b > 0.04045 ? Math.pow(((b + 0.055) / 1.055), 2.4) : (b / 12.92);

	var rgbMatrix = math.matrix([r, g, b]);
	var ret = math.multiply(rgb2xyz(),rgbMatrix);
	ret = math.multiply(ret,100);
	return ret;
}
rgb2xyz(){
	return math.matrix(
		[[0.4124, 0.3576, 0.1805],
		[0.2126, 0.7152, 0.0722],
		[0.0193, 0.1192, 0.9505]]);
}

```
これで0~1に正規化され、D65を基準としたXYZ値を得ることができた。

### あぁ、次はLMS値（錐体反応値）だ…

人間の眼は光を感じる為に網膜という器官がある。網膜には光を感じる細胞が大まかに二種類あり、桿体と錐体と呼ばれる。
桿体は光の量が少ない（周囲が暗い）時に反応する細胞で、一種類だけ存在している。一方、錐体は光の量が多い（周囲が明るい）時に反応する細胞で、三種類存在している。三種類の錐体それぞれが異なる波長に強く反応することで、人間は**三種類の原色**を感知することができる。
三種類の原色はそれぞれ**橙（赤）、緑、青**に対応している。（つまりほぼRGBと同じなのだが、ディスプレイはこれらの混色ですべての色を再現できるという理論に基づいて作られているので、当たり前といえば当たり前である。）
XYZ値と各錐体が反応する反応値を変換するには、各錐体がどのような波長に反応するかの分布から求めることができる。ここではVos(1978)の調査結果である錐体反応値を用いる。


```math
\left[
 \begin{array}{c}
  L\\
  M\\
  S
 \end{array}
\right]
=
\left[
 \begin{array}{rrr}
  0.15516 & 0.54307 & -0.03701\\
  -0.15516 & 0.45692 & 0.02969\\
  0 & 0 & 0.00732
 \end{array}
\right]
\left[
 \begin{array}{c}
  X_{D65}\\
  Y_{D65}\\
  Z_{D65}
 \end{array}
\right]
```

JSでは...

```.js
transXYZtoLMS(inXYZ){
	var Vos = VosMatrix();
	var LMS = math.multiply(Vos,inXYZ);
	return LMS;
}
VosMatrix(){
	//LMS錐体の反応値への変換行列
	//Vos(1978)の錐体分光感度を採用
	return math.matrix([
		[0.15516, 0.54307, -0.03701],
		[-0.15516, 0.45692, 0.02969],
		[0, 0, 0.00732]]);
};

```



### ようやく、色覚タイプ別に変換

色弱は錐体の反応が弱くなるor無くなることで起きる。「反応が弱くなる」にもさまざまな程度があるが、ここでは単純に「錐体が完全に反応しなくなった」場合を考える。
Hans Brettleら(1997)の報告[^2]でP型色弱者とD型色弱者は475nm(青)の光と575nm(黄色)の光の混色ですべての色を表現でき、T型色弱者は485nm(青緑)の光と660nm(赤)の光の混色ですべての色を表現できることがわかっている。また、黒、灰色、白といった無彩色はC型色覚(正常色覚、色弱ではない状態の人)と同じ見え方であることも分かっている。
なので、P型とD型の色の見え方を再現するには、C型のLMS値を475nmと白と575nmで構成された値に変換することで再現できる。各波長でのLMS値を算出するには、指定の波長でのXYZ値を先述の変換行列でLMS値にしてやればよい。ここでは、日本工業規格(JIS Z 8701 : 1995)で開示されているXYZ値を使用して算出した。


[^2]: Hans Brettle, Francoise Vienot, John Deutan. Mollon：Computerized simulation of color appearance for dichromats：J.Opt.Soc.Am.A, Vol.14, No.10, pp. 2647-2655 (1997)


```math
\left[
 \begin{array}{c}
  L_{Protan,Deutan}\\
  M_{Protan,Deutan}\\
  S_{Protan,Deutan}
 \end{array}
\right]
 = \left\{
\begin{array}{ll}
 \left[
 \begin{array}{c}
  0.044637199\\
  0.060334967\\
  0.007626708
 \end{array}
\right]
 & (波長475nm) \\
 \left[
 \begin{array}{c}
  0.627781960\\
  0.287595710\\
  0.000013176
 \end{array}
\right]
 & (波長575nm)
\end{array}
\right.
```
T型の色の見え方を再現するには、C型のLMS値を485nmと白と660nmで構成された値に変換することで再現できる。


```math

\left[
 \begin{array}{c}
  L_{Tritan}\\
  M_{Tritan}\\
  S_{Tritan}
 \end{array}
\right]
 = \left\{
\begin{array}{ll}
 \left[
 \begin{array}{c}
  0.078135469\\
  0.086652254\\
  0.004510584
 \end{array}
\right]
 & (波長485nm) \\
 \left[
 \begin{array}{c}
  0.058713154\\
  0.002286236\\
  0.000000000
 \end{array}
\right]
 & (波長660nm)
\end{array}
\right.
```

LMS値の各値を軸としたLMS空間で考えた場合、各波長と白は3本の線となり、それらを辺とした平面が作られる。
P型はL錐体が反応しない状態のため、空間上の任意の点はL軸に沿って平面に投影された点と同等の色となる。D型はM錐体が反応しないため、M軸に沿って投影された点と同等の色になる。

### 数学的な表現をするなら...
LMS色空間における変換前の値$Q(L_Q,M_Q,S_Q)$，変換後の値を$ Q'(L_Q′,M_Q′,S_Q′)$，白色を$ E(L_E,M_E,S_E)$，投影平面を構成する波長の成分を$A(L_A,M_A,S_A) $とした場合、

```math
(E \times A) \cdot Q' = 0\\
aL_Q′ + bM_Q′ + cS_Q′ = 0\\

```
ただし、

```math
a = M_ES_A − S_EM_A\\
b = S_EL_A − L_ES_A\\
c = L_EM_A − M_EL_A
```

となる。
P型の場合はL軸に沿って投影させるので

```math
L_Q′ = − \frac{bM_Q + cS_Q}{a}

```

D型の場合はM軸に沿って投影させるので

```math
M_Q′ = − \frac{aL_Q + cS_Q}{b}
```
となる。

T型の場合はS軸に沿って投影させるので

```math
SQ′ = − \frac{aL_Q + bM_Q}{c}
```
となる。

Protan変換のJSは

```.js
transLMStoProtan(inLMS){
	var eXYZ = eMatrix();
	var eLMS = transXYZtoLMS(eXYZ);
		
	//ProtanはLMS空間のL軸情報が喪失するのでMS平面に投影される。
	//入力値の傾きが等エネルギー白色より小さい場合は575nm、大きい場合は475nmに投影される。
	//Protanで投影される面(575nmとeLMSの平面 or 475nmとeLMSの平面)を調べる。
	var inTan = inLMS.get([2])/inLMS.get([1]);
	var eTan = eLMS.get([2])/eLMS.get([1]);
	var dichromatLMS;
	
	if(inTan<eTan){
		dichromatLMS = LMS575nm();
	}
	else{
		dichromatLMS = LMS475nm();
	}
	
	//二色型色覚の平面に投影するための係数を算出
	var abc = culcDichromatCoefficient(eLMS,dichromatLMS);
	
	//Protan変換の場合はL軸に沿って投影させる
	var protanL = (abc.get([1])*inLMS.get([1]) + abc.get([2])*inLMS.get([2]))/abc.get([0])*-1;		
	return math.matrix([ protanL, inLMS.get([1]), inLMS.get([2]) ]);
}

eMatrix(){
	//XYZ表色系における等エネルギー白色
	return math.matrix(	[0.333, 0.333, 0.333] );
}

LMS575nm(){
	return math.matrix(	[0.044637199, 0.060334967, 0.007626708] );
}
LMS475nm(){
	return math.matrix(	[0.627781960, 0.287595710, 0.000013176] );
}

culcDichromatCoefficient(inE,inDichro){
	var a = inE.get([1]) * inDichro.get([2]) - inE.get([2]) * inDichro.get([1]);
	var b = inE.get([2]) * inDichro.get([0]) - inE.get([0]) * inDichro.get([2]);
	var c = inE.get([0]) * inDichro.get([1]) - inE.get([1]) * inDichro.get([0]);
	return math.matrix([a,b,c]);
}

```


## XYZ、RGBに戻していく。
ここまでで色弱者が見ている色の数学的、光学的な値がわかったので、今度はRGBに戻していく。
この処理は今まで行ってきた処理の逆なので、行列では逆行列をかけることと同じである。
つまり、LMSをXYZに戻す計算は

```math
\left[
 \begin{array}{c}
  X_{D65}\\
  Y_{D65}\\
  Z_{D65}
 \end{array}
\right]
=
\left[
 \begin{array}{rrr}
  0.15516 & 0.54307 & -0.03701\\
  -0.15516 & 0.45692 & 0.02969\\
  0 & 0 & 0.00732
 \end{array}
\right]^{-1}
\left[
 \begin{array}{c}
  L_Q'\\
  M_Q'\\
  S_Q'
 \end{array}
\right]
```

となる。
つまり、

```.js
transLMStoXYZ(inLMS){
	var invVos = math.inv(VosMatrix());
	var mColor = math.multiply(invVos,inLMS);
	return mColor;
}
```
となる。行列で計算処理をするとコードの流用率が高まり、良い感じ。
次は、XYZをRGBに戻す。これも同様に逆行列をかけることになるので

```math
\left[
 \begin{array}{r}
  R_{linear} \\
  G_{linear} \\
  B_{linear}
 \end{array}
\right]
=
\left[
 \begin{array}{rrr}
  0.4124 & 0.3576 & 0.1805\\
  0.2126 & 0.7152 & 0.0722\\
  0.0193 & 0.1192 & 0.9505\\
 \end{array}
\right]^{-1}
\left[
 \begin{array}{r}
  X_{D65} \\
  Y_{D65} \\
  Z_{D65}
 \end{array}
\right]
```
となる。
このRGBは0~1に正規化され、かつ、ガンマ補正が掛けられた状態なので表示用のRBGにする必要がある。
まずは、ガンマ補正を戻すために

```math
f(rgb)^{-1} = \left\{
\begin{array}{ll}
 12.92 \times rgb & (rgb \leq 0.0032308) \\
 (1.055 \times rgb)^{ \frac{1}{2.4}} - 0.055 & (otherwise)
\end{array}
\right.
```

を行う。
最後に0~255にするために

```math
R = 255 \times R_{linear}\\
G = 255 \times G_{linear}\\
B = 255 \times B_{linear}
```
を処理することで各色覚タイプでのRGBとなる。


```.js
transXYZtoRGB(inXYZ){
	var normalizedXYZ = math.divide(inXYZ,100);
	var xyz2rgb = math.inv(rgb2xyz());
	var mColor = math.multiply(xyz2rgb,normalizedXYZ);

	var rli = mColor.get([0]);
	var gli = mColor.get([1]);
	var bli = mColor.get([2]);

	var r = rli>0.0032308 ? (1.055*math.pow(rli,1/2.4)-0.055) : rli*12.92;
	var g = gli>0.0032308 ? (1.055*math.pow(gli,1/2.4)-0.055) : gli*12.92;
	var b = bli>0.0032308 ? (1.055*math.pow(bli,1/2.4)-0.055) : bli*12.92;

	var ret = math.multiply(math.matrix([r,g,b]),255);
	var ret = math.round(ret);

	return ret;
}
```

これらを実装したのが[こちら](https://dichrosim.web.app/index_sRGB.html)(色差判定機能付き)。
色覚タイプなどの変換、Labへの変換に関してはそれなりの精度(丸め誤差程度)が出ているが、sRGB->XYZの誤差が大きいので随時改良したい。