分光反射率から色を計算してみよう [Python]

分光反射率計はあるけれど、測色計がない！という場面がありました。
まあ原理は同じだしな…ということで計算で何とかすることにしました。（他部署まで都度借りにいくのは面倒ですし）

いろいろなウェブサイトにエッセンスは載っていますが、このような目的のためにまとめられたものは見当たらず、結構苦戦したので備忘録もかねて記事にしておきます。

具体的には分光反射率のスペクトルデータから XYZ、L*a*b*、RGB を計算します。前置きはこれくらいにして早速中身に入ります。

必要なデータ

1. サンプルの分光反射率スペクトル：
波長 λ (nm) とそれに対応する分光反射率 R (%)
※必ず「全光線反射測定」を行いましょう。鏡面反射測定や拡散反射測定ではダメです。

2. CIE等色関数：
波長 λ (nm) とそれに対応する X, Y, Z (三刺激値)

3. D65光源の相対強度スペクトル：
波長 λ (nm) とそれに対応する光源の相対強度 S [a.u.]

それぞれ同じ波長の範囲・間隔のデータを準備してください。
波長 380 - 780 nm ( 5 or 10 nm 刻み) であれば十分でしょう。
2 と 3 は本記事の最後にも載せていますし、探せばいくらでも出てきます。 

分光反射率は光源の種類に依存しないので、ハロゲンランプなど、D65 光源以外で取得したデータでも全くかまいません。

---

計算の方法

細かい理屈は省きますが、簡単に計算方法を示します。
python のコードを載せているので、同じようにコーディングしてもらえれば数式については理解せずともなんとかなります。

分光スペクトルデータから XYZ への変換までがこの記事の本質ですが、どうせなら L*a*b*、RGB まで変換できたほうが使い勝手がよいので、そこまでを ゴールとします。

【Step 1】XYZ の算出
分光反射率 R(λ)、光源の強度分布 S(λ)、CIE 等色関数 $${\={x}(\lambda)}$$, $${\={y}(\lambda)}$$, $${\={z}(\lambda)}$$から下記のように算出できます。

$$
X = \frac{\Sigma_\lambda R(\lambda)S(\lambda)\={x}(\lambda)}{\Sigma_\lambda S(\lambda)\={y}(\lambda)} 
$$

$$
Y = \frac{\Sigma_\lambda R(\lambda)S(\lambda)\={y}(\lambda)}{\Sigma_\lambda S(\lambda)\={y}(\lambda)}
$$

$$
Z = \frac{\Sigma_\lambda R(\lambda)S(\lambda)\={z}(\lambda)}{\Sigma_\lambda S(\lambda)\={y}(\lambda)}
$$

【Step 2】L*a*b* の算出
上記で算出した XYZ は D65 光源下での色を示しているため、D65 の白色点の三刺激値 $${X_n, Y_n, Z_n}$$ で割りつけて相対値に変換してください。

$$
X_n = 95.047,    Y_n = 100.000,    Z_n = 108.883
$$

L*a*b* は下式で計算することができます。

$$
f(t) = \begin{cases} t^{1/3} & \ t >(6/29)^3 = 0.008856…\text{のとき} \\
\frac{1}{3}(\frac{29}{6})^2t + \frac{4}{29} & \text{それ以外} \end{cases}
$$

$$
L^* = 116 f\left( \frac{Y}{Y_n} \right) -16
$$

$$
a^* = 500 f\left( \frac{X}{X_n} - \frac{Y}{Y_n} \right)
$$

$$
b^* = 200 f\left( \frac{Y}{Y_n} - \frac{Z}{Z_n} \right)
$$

【Step 3】sRGB の算出
XYZ から sRGB への変換を行います。これが非常に面倒くさい。
今回は D65 光源の白色点を使うので、下記の変換式を使ってまずは linear sRGB に変換します。

$$
\begin{bmatrix} R_{sRGB}\\  G_{sRGB}\\ B_{sRGB} \end{bmatrix}
= \begin{bmatrix}  3.2406 & -1.5372 & -0.4986\\  -0.9689 & 1.8758 & 0.0415\\0.0557 & -0.2040 & 1.0570 \end{bmatrix}
\begin{bmatrix} X/100\\  Y/100\\ Z/100 \end{bmatrix}
$$

続いて、γ = 2.4 としてガンマ補正を行います。

$$
R'_{sRGB} = \begin{cases} 12.92 \cdot R_{sRGB} & R_{sRGB} \leq 0.0031308 \text{のとき }\\
1.055 \cdot (R_{sRGB}^{\frac{1}{2.4}} - 0.055 & \text{それ以外} \end{cases}
$$

$$
G'_{sRGB} = \begin{cases} 12.92 \cdot G_{sRGB} & G_{sRGB} \leq 0.0031308 \text{のとき }\\
1.055 \cdot (G_{sRGB}^{\frac{1}{2.4}} - 0.055 & \text{それ以外} \end{cases}
$$

$$
B'_{sRGB} = \begin{cases} 12.92 \cdot B_{sRGB} & B_{sRGB} \leq 0.0031308 \text{のとき }\\
1.055 \cdot (B_{sRGB}^{\frac{1}{2.4}} - 0.055 & \text{それ以外} \end{cases}
$$

この時点で 0 未満の値、1 を超える値は RGB で表すことができないため、それぞれ 0, 1 に置換 (クリップ) させてしまいます。

最後に、8ビットのデータ (0 - 255) となるよう変換します。

$$
R_{sRGB}(8) = 255\cdot R'_{sRGB}
$$

$$
G_{sRGB}(8) = 255\cdot G'_{sRGB}
$$

$$
B_{sRGB}(8) = 255\cdot B'_{sRGB}
$$

以上で目的の数値が計算できました。

---

Python での実行

頑張れば Excel でも計算できないことはないとは思いますが、非常に可読性が悪くなるので、Python で処理しましょう。以下がサンプルコードです。

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

#----------XYZを計算する----------

# CSVファイルの読み込みと列名の付与
sample_data = pd.read_csv('sample_data.csv', header=0, names=['wavelength', 'reflectance'])
color_matching_functions = pd.read_csv('color_matching_functions.csv', header=0, names=['wavelength', 'X', 'Y', 'Z'])
d65_data = pd.read_csv('D65_distribution.csv', header=0, names=['wavelength', 'intensity'])

# サンプルの分光スペクトルデータ（波長、分光反射率）
wavelengths = sample_data['wavelength'].values
reflectance = sample_data['reflectance'].values

# CIE等色関数（波長、X、Y、Z）
cie_x = color_matching_functions['X'].values
cie_y = color_matching_functions['Y'].values
cie_z = color_matching_functions['Z'].values

# D65光源の分光スペクトルデータ（波長、強度）
d65 = d65_data['intensity'].values

# XYZの計算
cie_y_sum = np.sum(d65 * cie_y)

X = np.sum(reflectance * d65 * cie_x) /cie_y_sum
Y = np.sum(reflectance * d65 * cie_y) /cie_y_sum
Z = np.sum(reflectance * d65 * cie_z) /cie_y_sum
print(f"X, Y, Z = {X, Y, Z}")

#----------L*a*b*を計算する----------

# D65の白色点
Xn, Yn, Zn = 95.047, 100.000, 108.883

# f(t)の定義
def f(t):
    delta = 6/29
    if t > delta**3:
        return t**(1/3)
    else:
        return t/(3*delta**2) + 4/29

# L*a*b*の計算
L = 116 * f(Y/Yn) - 16
a = 500 * (f(X/Xn) - f(Y/Yn))
b = 200 * (f(Y/Yn) - f(Z/Zn))

print(f"L*, a*, b* = {L, a, b}")

#----------RGBを計算する----------

def xyz_to_rgb(x, y, z):
    # XYZからRGBへの変換行列（D65基準）
    M = np.array([[ 3.2406, -1.5372, -0.4986],
                  [-0.9689,  1.8758,  0.0415],
                  [ 0.0557, -0.2040,  1.0570]])
    
    # XYZ値を行列と掛け算してRGB値を計算
    rgb = M @ np.array([x, y, z])
    
    # ガンマ補正
    rgb = np.where(rgb <= 0.0031308, rgb * 12.92, 1.055 * (rgb ** (1/2.4)) - 0.055)
    
    # RGBの範囲にクリッピング
    rgb = np.clip(rgb, 0, 1)
    
    return rgb

def display_color(rgb):
    # RGB値を0-255の範囲にスケーリング
    rgb = np.array(rgb) * 255
    
    # 色として表示
    plt.imshow([[rgb/255]])
    plt.axis('off')
    plt.show()

# XYZのスケーリング
X_norm = X / 100
Y_norm = Y / 100
Z_norm = Z / 100

# XYZからRGBに変換
rgb = xyz_to_rgb(X_norm, Y_norm, Z_norm)

# 色として表示
print(f"R, G, B = {rgb *255}")
display_color(rgb)

例として、りんごの分光反射率のデータから色を算出してみます。

実行結果

X, Y, Z = (36.740130682862116, 24.466222806971523, 10.524196512458953)
L*, a*, b* = (56.551595022603266, 51.50345755578145, 33.30369481376142)
R, G, B = [226.19815739  92.12650238  80.76438871]

無事においしそうなりんごの色が出力されました。
めでたし。

---

付録① D65光源の相対強度データ

$$
\begin{array}{c:c}
波長(nm) & 相対強度 \\
\hline\
380 & 50.614 \\ 385 & 52.608 \\ 390 & 54.602 \\ 395 & 68.649 \\ 400 & 82.696 \\ 405 & 87.060 \\ 410 & 91.424 \\ 415 & 92.398 \\ 420 & 93.373 \\ 425 & 90.002 \\ 430 & 86.632 \\ 435 & 95.723 \\ 440 & 104.814 \\ 445 & 110.886 \\ 450 & 116.958 \\ 455 & 117.362 \\ 460 & 117.767 \\ 465 & 116.294 \\ 470 & 114.822 \\ 475 & 115.355 \\ 480 & 115.889 \\ 485 & 112.336 \\ 490 & 108.783 \\ 495 & 109.057 \\ 500 & 109.332 \\ 505 & 108.557 \\ 510 & 107.783 \\ 515 & 106.280 \\ 520 & 104.777 \\ 525 & 106.229 \\ 530 & 107.681 \\ 535 & 106.040 \\ 540 & 104.400 \\ 545 & 104.221 \\ 550 & 104.043 \\ 555 & 102.022 \\ 560 & 100.000 \\ 565 & 98.168 \\ 570 & 96.336 \\ 575 & 96.065 \\ 580 & 95.793 \\ 585 & 92.242 \\ 590 & 88.692 \\ 595 & 89.354 \\ 600 & 90.016 \\ 605 & 89.814 \\ 610 & 89.612 \\ 615 & 88.662 \\ 620 & 87.713 \\ 625 & 85.509 \\ 630 & 83.305 \\ 635 & 83.512 \\ 640 & 83.718 \\ 645 & 81.882 \\ 650 & 80.046 \\ 655 & 80.141 \\ 660 & 80.237 \\ 665 & 81.270 \\ 670 & 82.303 \\ 675 & 80.306 \\ 680 & 78.309 \\ 685 & 74.026 \\ 690 & 69.743 \\ 695 & 70.688 \\ 700 & 71.633 \\ 705 & 73.003 \\ 710 & 74.372 \\ 715 & 67.998 \\ 720 & 61.623 \\ 725 & 65.765 \\ 730 & 69.906 \\ 735 & 72.508 \\ 740 & 75.109 \\ 745 & 69.360 \\ 750 & 63.611 \\ 755 & 55.022 \\ 760 & 46.432 \\ 765 & 56.629 \\ 770 & 66.825 \\ 775 & 65.113 \\ 780 & 63.401 \\
\end{array}
$$

---

付録② CIE 等色関数

$$
\begin{array}{c:cccc}
波長(nm) & X & Y & Z \\
\hline\
380 & 0.0014 & 0 & 0.0065 \\ 385 & 0.0022 & 0.0001 & 0.0105 \\ 390 & 0.0042 & 0.0001 & 0.0201 \\ 395 & 0.0077 & 0.0002 & 0.0362 \\ 400 & 0.0143 & 0.0004 & 0.0679 \\ 405 & 0.0232 & 0.0006 & 0.1102 \\ 410 & 0.0435 & 0.0012 & 0.2074 \\ 415 & 0.0776 & 0.0022 & 0.3713 \\ 420 & 0.1344 & 0.004 & 0.6456 \\ 425 & 0.2148 & 0.0073 & 1.0391 \\ 430 & 0.2839 & 0.0116 & 1.3856 \\ 435 & 0.3285 & 0.0168 & 1.623 \\ 440 & 0.3483 & 0.023 & 1.7471 \\ 445 & 0.3481 & 0.0298 & 1.7826 \\ 450 & 0.3362 & 0.038 & 1.7721 \\ 455 & 0.3187 & 0.048 & 1.7441 \\ 460 & 0.2908 & 0.06 & 1.6692 \\ 465 & 0.2511 & 0.0739 & 1.5281 \\ 470 & 0.1954 & 0.091 & 1.2876 \\ 475 & 0.1421 & 0.1126 & 1.0419 \\ 480 & 0.0956 & 0.139 & 0.813 \\ 485 & 0.058 & 0.1693 & 0.6162 \\ 490 & 0.032 & 0.208 & 0.4652 \\ 495 & 0.0147 & 0.2586 & 0.3533 \\ 500 & 0.0049 & 0.323 & 0.272 \\ 505 & 0.0024 & 0.4073 & 0.2123 \\ 510 & 0.0093 & 0.503 & 0.1582 \\ 515 & 0.0291 & 0.6082 & 0.1117 \\ 520 & 0.0633 & 0.71 & 0.0782 \\ 525 & 0.1096 & 0.7932 & 0.0573 \\ 530 & 0.1655 & 0.862 & 0.0422 \\ 535 & 0.2257 & 0.9149 & 0.0298 \\ 540 & 0.2904 & 0.954 & 0.0203 \\ 545 & 0.3597 & 0.9803 & 0.0134 \\ 550 & 0.4334 & 0.995 & 0.0087 \\ 555 & 0.5121 & 1 & 0.0057 \\ 560 & 0.5945 & 0.995 & 0.0039 \\ 565 & 0.6784 & 0.9786 & 0.0027 \\ 570 & 0.7621 & 0.952 & 0.0021 \\ 575 & 0.8425 & 0.9154 & 0.0018 \\ 580 & 0.9163 & 0.87 & 0.0017 \\ 585 & 0.9786 & 0.8163 & 0.0014 \\ 590 & 1.0263 & 0.757 & 0.0011 \\ 595 & 1.0567 & 0.6949 & 0.001 \\ 600 & 1.0622 & 0.631 & 0.0008 \\ 605 & 1.0456 & 0.5668 & 0.0006 \\ 610 & 1.0026 & 0.503 & 0.0003 \\ 615 & 0.9384 & 0.4412 & 0.0002 \\ 620 & 0.8544 & 0.381 & 0.0002 \\ 625 & 0.7514 & 0.321 & 0.0001 \\ 630 & 0.6424 & 0.265 & 0 \\ 635 & 0.5419 & 0.217 & 0 \\ 640 & 0.4479 & 0.175 & 0 \\ 645 & 0.3608 & 0.1382 & 0 \\ 650 & 0.2835 & 0.107 & 0 \\ 655 & 0.2187 & 0.0816 & 0 \\ 660 & 0.1649 & 0.061 & 0 \\ 665 & 0.1212 & 0.0446 & 0 \\ 670 & 0.0874 & 0.032 & 0 \\ 675 & 0.0636 & 0.0232 & 0 \\ 680 & 0.0468 & 0.017 & 0 \\ 685 & 0.0329 & 0.0119 & 0 \\ 690 & 0.0227 & 0.0082 & 0 \\ 695 & 0.0158 & 0.0057 & 0 \\ 700 & 0.0114 & 0.0041 & 0 \\ 705 & 0.0081 & 0.0029 & 0 \\ 710 & 0.0058 & 0.0021 & 0 \\ 715 & 0.0041 & 0.0015 & 0 \\ 720 & 0.0029 & 0.001 & 0 \\ 725 & 0.002 & 0.0007 & 0 \\ 730 & 0.0014 & 0.0005 & 0 \\ 735 & 0.001 & 0.0004 & 0 \\ 740 & 0.0007 & 0.0002 & 0 \\ 745 & 0.0005 & 0.0002 & 0 \\ 750 & 0.0003 & 0.0001 & 0 \\ 755 & 0.0002 & 0.0001 & 0 \\ 760 & 0.0002 & 0.0001 & 0 \\ 765 & 0.0001 & 0 & 0 \\ 770 & 0.0001 & 0 & 0 \\ 775 & 0.0001 & 0 & 0 \\ 780 & 0 & 0 & 0 \\
\end{array}
$$

---

付録③ りんごの分光反射率（要るか？？）

$$
\begin{array}{c:c}
波長(nm) & 反射率(％) \\
\hline\
380 & 4.832 \\ 385 & 5.025 \\ 390 & 5.219 \\ 395 & 5.413 \\ 400 & 5.607 \\ 405 & 5.801 \\ 410 & 6.1 \\ 415 & 6.779 \\ 420 & 7.463 \\ 425 & 8.178 \\ 430 & 8.914 \\ 435 & 9.649 \\ 440 & 10.335 \\ 445 & 10.909 \\ 450 & 11.616 \\ 455 & 11.63 \\ 460 & 10.91 \\ 465 & 10.047 \\ 470 & 9.254 \\ 475 & 8.526 \\ 480 & 8.047 \\ 485 & 7.726 \\ 490 & 7.512 \\ 495 & 7.565 \\ 500 & 7.619 \\ 505 & 7.678 \\ 510 & 7.566 \\ 515 & 7.566 \\ 520 & 7.566 \\ 525 & 7.405 \\ 530 & 7.192 \\ 535 & 7.085 \\ 540 & 7.085 \\ 545 & 7.113 \\ 550 & 7.907 \\ 555 & 10.528 \\ 560 & 13.601 \\ 565 & 17.404 \\ 570 & 21.127 \\ 575 & 26.643 \\ 580 & 33.514 \\ 585 & 44.287 \\ 590 & 50.604 \\ 595 & 54.815 \\ 600 & 58.138 \\ 605 & 61.006 \\ 610 & 63.282 \\ 615 & 65.077 \\ 620 & 66.476 \\ 625 & 67.753 \\ 630 & 68.805 \\ 635 & 69.793 \\ 640 & 70.575 \\ 645 & 71.202 \\ 650 & 71.835 \\ 655 & 72.254 \\ 660 & 72.532 \\ 665 & 72.851 \\ 670 & 73.064 \\ 675 & 73.277 \\ 680 & 73.492 \\ 685 & 73.546 \\ 690 & 73.599 \\ 695 & 73.867 \\ 700 & 74.19 \\ 705 & 74.188 \\ 710 & 74.23 \\ 715 & 74.203 \\ 720 & 74.142 \\ 725 & 74.234 \\ 730 & 74.142 \\ 735 & 74.112 \\ 740 & 74.121 \\ 745 & 74.286 \\ 750 & 74.108 \\ 755 & 74.206 \\ 760 & 74.149 \\ 765 & 74.198 \\ 770 & 74.183 \\ 775 & 74.171 \\ 780 & 74.143 \\
\end{array}
$$