ぎるばーとの日記

もっともっと遠くへ行きたい 空が広く見える場所まで

色温度からRGB(sRGB)に変換(完成版)

何をするの?

 こういうの見たことありませんか? 色温度に対応するRGB(sRGB)値を求めるコードを実装します。

f:id:zwxadz:20170502053830p:plain

計算した値の一部

 相関色温度4000K〜12000KのCIE昼光がどんな色をしているか計算しました。
 100K刻みで表にしましたが、実装コードではその間の値も計算できます。また、最高25000Kまで入力OKです。
 特定の色温度のカラーコード値を知るのに役立つかもしれません。

4000K: #FFD6A1 6000K: #FFF9F0 8000K: #E0ECFF 10000K: #CADDFF
4100K: #FFD9A6 6100K: #FFFBF3 8100K: #DFEBFF 10100K: #C9DCFF
4200K: #FFDBAB 6200K: #FFFCF6 8200K: #DDEAFF 10200K: #C8DCFF
4300K: #FFDDAF 6300K: #FFFDF9 8300K: #DCE9FF 10300K: #C7DBFF
4400K: #FFDFB4 6400K: #FFFEFC 8400K: #DBE8FF 10400K: #C6DBFF
4500K: #FFE2B9 6500K: #FFFFFF 8500K: #D9E7FF 10500K: #C6DAFF
4600K: #FFE4BD 6600K: #FCFEFF 8600K: #D8E6FF 10600K: #C5DAFF
4700K: #FFE6C1 6700K: #FAFCFF 8700K: #D7E5FF 10700K: #C4D9FF
4800K: #FFE7C5 6800K: #F7FAFF 8800K: #D5E5FF 10800K: #C3D9FF
4900K: #FFE9C9 6900K: #F5F9FF 8900K: #D4E4FF 10900K: #C3D8FF
5000K: #FFEBCD 7000K: #F3F7FF 9000K: #D3E3FF 11000K: #C2D8FF
5100K: #FFEDD1 7100K: #F1F6FF 9100K: #D2E2FF 11100K: #C1D7FF
5200K: #FFEED5 7200K: #EEF5FF 9200K: #D1E2FF 11200K: #C1D7FF
5300K: #FFF0D9 7300K: #ECF4FF 9300K: #D0E1FF 11300K: #C0D6FF
5400K: #FFF1DC 7400K: #EBF2FF 9400K: #CFE0FF 11400K: #C0D6FF
5500K: #FFF3E0 7500K: #E9F1FF 9500K: #CEE0FF 11500K: #BFD6FF
5600K: #FFF4E3 7600K: #E7F0FF 9600K: #CDDFFF 11600K: #BED5FF
5700K: #FFF6E7 7700K: #E5EFFF 9700K: #CCDEFF 11700K: #BED5FF
5800K: #FFF7EA 7800K: #E3EEFF 9800K: #CBDEFF 11800K: #BDD4FF
5900K: #FFF8ED 7900K: #E2EDFF 9900K: #CADDFF 11900K: #BDD4FF

参考にしたWeb上の情報

 各温度の黒体の色を求めています。等色関数を介してXYZ色空間にマッピングするので手間がかかりますが、広い範囲の温度について計算できます。

 色の用語、特にCIE関連について詳しい解説があります。

  • JIS Z 8720:2012

 表題は「測色用の標準イルミナント(標準の光)及び標準光源」。ググると法的に灰色(?)っぽいサイトがヒットして見られます。リンクは貼りません。

 sRGBへの変換はWikipediaを参考にしました。

前提条件

CIE昼光(daylight)

 黒体(black body)ではなく、CIE昼光です。昼光とほぼ等しい光と考えてください。
 JISでは、「主として、北空についての多くの分光測定値から統計的手法によって各々の相関色温度における昼光の代表としてCIEが定めた分光分布」です。昼光の測定データに基づいて定めたスペクトルの光ということですね。
 色温度というとき、普通は黒体、つまりプランクの法則に従う放射を考えます。しかし、黒体と異なるスペクトルをもつ光についても、知覚的に色の近い黒体の温度で表すことができ、これを相関色温度といいます。

白色点はD65

 約6500Kが白になります。より正確には約6504Kなのですが、最終的に8bit/channelで整数に丸めるので差が出ません。
 JIS Z 8720にあるxy色度座標への変換式を使用します。定義域の4000K〜25000Kから外れると計算できません。(というか、定義されない……?)

sRGB符号化

 だいぶ前からsRGBが標準になってると思うので。
 Webでよくある"#XXXXXX"表記はsRGB前提ですよね?(たぶん)

処理の大筋

  1. 色温度の範囲チェック
  2. 色温度→xy色度座標
  3. xyz→線形sRGB
  4. 線形sRGBを規格化
  5. 線形sRGB→非線形sRGB
  6. 非線形sRGBを量子化(0〜255の整数)
入力チェックとxy色度座標への変換(1と2)

 JIS Z 8720にある、CIE昼光のxy色度座標を求める式は次の通り。

 相関色温度T_{\rm cp}が4000K以上7000K以下の場合の色度座標x_{\rm D}

x_{\rm D} = -4.6070\,\frac{10^9}{{T_{\rm cp}}^3} + 2.9678\,\frac{10^6}{{T_{\rm cp}}^2} + 0.09911\,\frac{10^3}{T_{\rm cp}} + 0.244063

 相関色温度T_{\rm cp}が7000Kを超え25000K以下の場合の色度座標x_{\rm D}

x_{\rm D} = -2.0064\,\frac{10^9}{{T_{\rm cp}}^3} + 1.9018\,\frac{10^6}{{T_{\rm cp}}^2} + 0.24748\,\frac{10^3}{T_{\rm cp}} + 0.237040

 色度座標y_{\rm D}

y_{\rm D} = -3.000{x_{\rm D}}^2 + 2.870x_{\rm D} - 0.275

 4000K未満、25000K超えは、とりあえず例外投げておけばいいかな……。

線形sRGBへの変換(3)

 まず、色度座標z_{\rm D}を求めます。

z_{\rm D} = 1 - x_{\rm D} - y_{\rm D}

 ここでは色度(色合い)だけに興味があって、輝度は気にしないでいいので、xyzをそのままXYZとして計算します。(xyzはx+y+z = 1となるよう、XYZをスケーリングしたもの。)

\begin{bmatrix}
R_\text{linear} \\ G_\text{linear} \\ B_\text{linear}
\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 \\ Y \\ Z
\end{bmatrix}

線形sRGBを規格化(4)

 赤緑青の最大値が1となるようにスケーリングします。

\begin{bmatrix}
R_\text{scaled} \\ G_\text{scaled} \\ B_\text{scaled}
\end{bmatrix} = \frac{1}{\max \{R_\text{linear}, G_\text{linear}, B_\text{linear}\}} \begin{bmatrix}
R_\text{linear} \\ G_\text{linear} \\ B_\text{linear}
\end{bmatrix}

 なお、xyzからY = 1と仮定してXYZに変換、さらに線形sRGBに変換して1を超えた値をクリッピングする方法も考えられます。クリッピングによって色度が変わってしまうのを避けるため、スケーリングを選択しました。

非線形sRGBへの変換(5)

 いわゆるガンマ補正です。次式でチャンネル毎に変換します。

C_\text{sRGB} = \begin{cases}
12.92 C_\text{linear} & \text{if $C_\text{linear} \leq 0.0031308$,} \\
1.055 {C_\text{linear}}^{1/2.4} - 0.055 & \text{otherwise.}
\end{cases}

8bit/channelで量子化(6)

 最後に、0〜1の浮動小数点数を0〜255の整数に丸めます。四捨五入なり偶数丸めなりお好みの関数で。(コードでは四捨五入にしました。)

C_\text{int-sRGB} = \text{round} (255 C_\text{sRGB})

 実装コードはQiitaの記事へどうぞ!

qiita.com

三角形内の一様乱数

 頂点ABCの三角形内で一様な乱数はどうすればいいでしょうか?

A new method to simulate the triangular distribution
 前回と同じ文献。「Random vectors in the plane」の部分で触れられてます。この論文より先に知られていた方法らしく。

u \sim \text{Uniform}(0,1),\ v \sim \text{Uniform}(0,1);
\text{MIN} = \text{MIN}(u,v),\ \text{MAX} = \text{MAX}(u,v)として、
\vec{a} \text{MIN}+ \vec{b} (1 - \text{MAX}) + \vec{c} (\text{MAX} - \text{MIN})

 なぜこれでOKなのかよく分かりませんが、2次元の場合でプロットしてみると確かに三角形になる……。ふむ。

Point2D randomPoint(double x1, double y1, double x2, double y2, double x3, double y3) {
    double u = Math.random(), v = Math.random();
    
    double m1 = Math.min(u, v);
    double m2 = 1.0 - Math.max(u, v);
    double m3 = Math.max(u, v) - Math.min(u, v);
    double x = m1 * x1 + m2 * x2 + m3 * x3;
    double y = m1 * y1 + m2 * y2 + m3 * y3;
    return new Point2D.Double(x, y);
}

三角乱数

 マニアックな乱数記事を徐々に増やしていくキャンペーン。

 前回・前々回と、なぜ数学関数の使用を回避したいかを書いてなかったので、その補足から、、
 プログラムの数学関数、具体的にはサイン、コサイン、対数平方根などですが、処理が重いことが多いです。除法も加・減・乗法と比べれば若干重いとはいえ、これの回避は難しい……。
 数学関数を使わないかわりに、単位円内の一様乱数を作るための繰返しが入りました。その繰返し回数を見積もってみます。
 円と外接正方形の面積の比は\pi/4 \approx 0.785 \ldotsで、円内に点が入ったときに脱出。その逆数が繰返し回数の平均になります。幾何分布なので数回程度かかる場合もありますが、平均は1回ちょっとです。
 ちなみに、実際にサイン、コサインを使うより速いかは環境によりけりかも……?

 ここからが今日の本文。
 三角分布をPERTで使うという話は本当なのかどうか知りません。

Triangular distribution - Wikipedia, the free encyclopedia
 三角乱数、平方根を使う方法がWikipediaに書いてあります。が、もっといいやり方(感動的!)があるので紹介します。

A new method to simulate the triangular distribution
u \sim \text{Uniform}(0,1),\ v \sim \text{Uniform}(0,1)として、
(1-c)\text{MIN}(u,v) + c\text{MAX}(u,v)
 これで、三角乱数(a=0, b=1, cは最頻値で0以上1以下)になります。MINとMAXは比較だけなので、平方根の計算よりずっと軽い!

double randomTriangular(double c) {
    double u = Math.random(), v = Math.random();
    
    return (1.0 - c) * Math.min(u, v) + c * Math.max(u, v);
}

 この方法が紹介されているのをあまり見ませんが、美しいアルゴリズムの一つだと思います。

コーシー乱数

 連日の更新なんていつぶり……?

 コーシー分布について。ちなみにコーシーと孔子は別人です。(豆知識
 コーシー分布は自由度1のt分布、平均も分散も発散するたちが悪い分布として名を知られてますね。
 和に対して安定(和の分布もコーシー分布)なのは正規分布と同じです。が、正規分布はスケールの二乗(分散)が加法的で、コーシー分布はスケールが加法的になります。

 分散1の独立な正規乱数の和は分散2の正規乱数になって、スケール(標準偏差)で見ると√2倍。さらに標本平均の分散は、
\text{Var}(\bar X) = \text{Var}(\frac{\sum X_i}{n}) = \frac{1}{n^2}\text{Var}(\sum X_i) = \frac{n\sigma^2}{n^2} = \frac{\sigma^2}{n}
 標本サイズが大きくなることで標本平均の分散が小さくなり、母平均の推定の信頼性が高まるのでした。

 片や、スケール1の独立なコーシー乱数の和はスケール2のコーシー乱数になります。n個の和をとって、スケールnのコーシー乱数。nで割るとスケール1のコーシー乱数。と、標本平均で元の分布に戻ってしまいました。母確率分布がコーシー分布のとき、標本サイズを大きくしても母中心の推定の信頼性は高まりません。

Cauchy distribution - Wikipedia, the free encyclopedia
 コーシー乱数の生成には、逆関数法とか標準正規乱数の比をとるとかいろいろ見つかりますが、ここでは数学関数を使わない方法を……。

double randomCauchy() {
    double u, v;
    do {
        u = 2.0 * Math.random() - 1.0;
        v = Math.random();
    } while (u * u + v * v >= 1.0 || v == 0.0);
    
    return u / v;
}

 単位円内(上半分)からランダムに選んだ点で、X座標をY座標で割るとコーシー分布(0,1)の乱数が得られます。説明は割愛、、Box–Mullerとかで検索するといいんじゃないかな……。

逆正弦乱数

physnotes.jp

 逆正弦分布(アークサイン分布)という確率分布があるのを、高校物理の備忘録さんの記事で知りました。記事の後半の逆正弦法則の話もとても興味深いです。

 振り子の速度は最低点で最大、両端では0になります。なので、最低点付近にいる時間は短く、両端付近にいる時間が長くなり、ランダムなタイミングでチラッと見たときの位置は逆正弦分布に従います。

Arcsine distribution - Wikipedia, the free encyclopedia
 逆正弦分布(0,1)は、ベータ分布の特別な場合(α=β=0.5)の別名です。この分布の乱数を生成する方法(用途はまったく不明、、)を書き記しておきます。

 Wikipediaにある通り、一様乱数U \sim \text{Uniform}(-\pi,\pi)から、(\sin U)^2とすればいいわけですが、サイン関数を使わずに計算する方法もあります。正規乱数でおなじみのMarsaglia法を応用します。

double randomArcsine() {
    double u, v, s;
    do {
        u = Math.random();
        v = Math.random();
        s = u * u + v * v;
    } while (s >= 1.0 || s == 0.0);
    
    return (u * u) / s;
}

 第1象限で単位円内の点をランダムに選んで、原点からの距離で割ると円周上に分布します。座標を一つとって二乗すると目的の乱数が得られますね。正規乱数の場合のように二つ独立に作れるわけではないのに注意! (\sin x)^2 + (\cos x)^2 = 1なので……。