スペクトルレンダリング

![joined.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/69648c1a-75a8-d3e4-3a7d-65fbdcc4e99f.png)

# はじめまして
この記事は[レイトレアドベントカレンダー2021](https://qiita.com/advent-calendar/2021/raytracing)の記事として作成されました。
初めてレイトレアドベントカレンダーを書かせていただきます。質問や間違いがありましたらコメントにてお願いします。

## 動機
![image](https://upload.wikimedia.org/wikipedia/commons/6/66/Chromatic_aberration_%28comparison%29.jpg)
> Stan Zurek, CC BY-SA 3.0 <http://creativecommons.org/licenses/by-sa/3.0/>, via Wikimedia Commons

上図の様な色収差をレイトレーシングで計算したいなと思い、本記事でリアルなカメラレンズをシミュレーションしました。

#概要
光の周波数を考慮できるようにレンダリング方程式を拡張し、ガラスでの屈折率を周波数に依存させたりと他にも周波数に依存する物理現象を計算可能にします。おまけとして、計算の効率化についても解説します。
 
#レンダリング方程式
$$L_0(x, \vec{\omega_0}) = L_e(x, \vec{\omega_0}) + \int_{S^2} f_s(x, \vec{\omega_i}, \vec{\omega_0} ) L_i(x, \vec{\omega_i}) |\vec{\omega_i} \cdot \vec{n}| d\vec{\omega_i})$$
この有名なレンダリング方程式に、光の周波数を考慮させるようにします。
$$L_0(x, \vec{\omega_0}, \lambda) = L_e(x, \vec{\omega_0}, \lambda) + \int_{S^2} f_s(x, \vec{\omega_i}, \vec{\omega_0}, \lambda) L_i(x, \vec{\omega_i}, \lambda) |\vec{\omega_i} \cdot \vec{n}| d\vec{\omega_i})$$
ただ、周波数$\lambda$引数を追加しただけです。
これで特定の周波数の輝度が分かるようになります。

| |レンダリング方程式|光の周波数を考慮されたレンダリング方程式|
|-|-|-|
|入力|$x, \omega_o$|$x, \omega_o, \lambda$|
|出力|輝度(RGB)|特定の周波数の輝度|


## スペクトルから色情報に変換
周波数を考慮しないレンダリング方程式は計算結果をそのままRGBとして扱えますが、周波数を考慮した方程式は任意の周波数輝度スペクトラムをRGBに変換する必要があります。ここで`スペクトラム`から`RGB`に変換について説明します。
しかし、スペクトラムから直接RGBに変換すると被積分式が負の値を持つことがあるので少し複雑になってしまいます。なので、`スペクトラム`から`XYZ`に一旦変換した後に`RGB`に変換します。これによって非負な被積分式を扱うことができます。

下の式でスペクトラムからXYZに変換できます。$P(\lambda)$はスペクトラム、$vis$は可視光域、$\bar{x}(\lambda),\bar{y}(\lambda),\bar{z}(\lambda)$はxyzの等色関数です。等色関数の近似式は[Wikiページ](https://en.wikipedia.org/wiki/CIE_1931_color_space#Analytical_approximation)を参照
$$y = \int_{vis} \bar{y}(\lambda) d\lambda$$
$$X = \frac{1}{y} \int_{vis} P(\lambda) \cdot \bar{x}(\lambda) d\lambda$$
$$Y = \frac{1}{y} \int_{vis} P(\lambda) \cdot \bar{y}(\lambda) d\lambda$$
$$Z = \frac{1}{y} \int_{vis} P(\lambda) \cdot \bar{z}(\lambda) d\lambda$$

次に積分結果のXYZをRGBに変換します。この変換式はディスプレイのRGB色空間によって変わるので気を付ける必要があります。
sRGB(D65)の場合の変換式は以下です。

```math
\begin{bmatrix}
R \\
G \\
B
\end{bmatrix}
=
\begin{bmatrix}
3.2410 & -1.5374 & -0.4986 \\
-0.9692 & 1.8760 & 0.0416 \\
0.0556 & -0.2040 & 1.0507
\end{bmatrix}
\begin{bmatrix}
X \\
Y \\
Z
\end{bmatrix}
```

他の色空間ヘの変換は[このサイト](http://k-ichikawa.blog.enjoy.jp/etc/HP/js/CIEXYZ2/cie3.html)が参考になります。
これで周波数を考慮したレンダリング方程式の計算結果をRGBを得る事ができるようになりました。

## 導出
前の章で計算の流れは分かったと思いますが、実際にプログラムに書くには$L_0(x, \vec{\omega_0}, \lambda)$の具体的な式が未定なので定積分で解くことはできません。よってモンテカルロ法を用いて数値計算で解きます。

$$X(x, \vec{\omega_0}) = \frac{1}{y} \int_{vis} L_0(x, \vec{\omega_0}, \lambda) \cdot \bar{x}(\lambda) d\lambda$$

一様に可視光の領域をサンプルするようなモンテカルロ法を実行します。
まず、一様なPDFを準備します。$C$は定数です。
$$p(\lambda) := C$$
定義域内で積分して1にならないといけないので、

$$\int_{vis} p(\lambda) d\lambda = 1$$

$$\int_{vis} C d\lambda = 1$$

$$C = \frac{1}{\int_{vis} d\lambda}$$

$vis \in [380, 780]$より

$$C = \frac{1}{400}$$

$$p(\lambda) = \frac{1}{400}$$

$p(\lambda)$のCDF$P(\lambda)$ は

$$P(\lambda) = \int_{380}^\lambda \frac{1}{400} d\lambda'$$
$$= \frac{\lambda - 380}{400}$$

逆関数法より
$\xi_\lambda$は[0, 1)範囲の一様乱数

$$P^{-1}(\xi_\lambda) = 400\xi_\lambda + 380$$

$$X(x, \vec{\omega_0}) = \frac{1}{y} \int_{vis} \frac{L_0(x, \vec{\omega_0}, \lambda) \cdot \bar{x}(\lambda)}{p(\lambda)} p(\lambda) d\lambda$$

モンテカルロ法を用いて
$$X(x, \vec{\omega_0}) \approx \frac{1}{y} \frac{1}{N} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{x}(\lambda_i)}{p(\lambda_i)} $$

$$X(x, \vec{\omega_0}) \approx \frac{1}{y} \frac{400}{N} \sum_{i=1}^N L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{x}(\lambda_i) $$

これで数値計算でXの期待値を得る事ができ、Y,Zも同様にします。

$$Y(x, \vec{\omega_0}) \approx \frac{1}{y} \frac{400}{N} \sum_{i=1}^N L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{y}(\lambda_i) $$
$$Z(x, \vec{\omega_0}) \approx \frac{1}{y} \frac{400}{N} \sum_{i=1}^N L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{z}(\lambda_i) $$

## プログラム
```cpp
struct Vec3{
    double x,y,z;
};
Vec3 xbybzbFromWavelength(double lambda) {
    Vec3 xyz;
    xyz.x = 1.056*gausse(lambda, 599.8, 37.9, 31.0) + 0.362*gausse(lambda, 442.0, 16.0, 26.7) - 0.065*gausse(lambda, 501.1, 20.4, 26.2);
    xyz.y = 0.821*gausse(lambda, 568.8, 46.9, 40.5) + 0.286*gausse(lambda, 530.9, 16.3, 31.1);
    xyz.z = 1.217*gausse(lambda, 437.0, 11.8, 36.0) + 0.681*gausse(lambda, 459.0, 26.0, 13.8);
    return xyz;
}

Vec3 SpectrumPathTracing(Vec3 point, Vec3 direction) {
    double integral_y = 106.91973463815504343620398716858; // =int_{vis} \bar{y}(\lambda) d\lambda
    
    Vec3 spectrumRadiance;
    for (int s=0;s<samples;++s) {
        double wavelength = sampleXYZ();
        Vec3 xyz = xbybzbFromWavelength(wavelength);
        double incomingRadiance = PathTracing(wavelength);
        spectrumRadiance.x += incomingRadiance * 400 * xyz.x / integral_y / samples ;
        spectrumRadiance.y += incomingRadiance * 400 * xyz.y / integral_y / samples;
        spectrumRadiance.x += incomingRadiance * 400 * xyz.z / integral_y / samples ;
    }
    return spectrumRadiance;
}
```

## 結果

|![result_RGB.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/3be56394-eb95-bf5a-9814-d47a447f2d73.jpeg)|
|:-:|
|通常レンダリング|


|![result_spectrum.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/74eca1c6-12ba-16f9-e504-781e62203988.jpeg)|
|:-:|
|スペクトルレンダリング|
上図は、周波数依存で屈折率が変わるガラスオブジェクト、カメラレンズでスペクトルレンダリングした結果です。

|スペクトルのサンプル数|パストレーシングのサンプル数|スーパーサンプリング|
|-|-|-|
|200|100|4|

|![result_spectrum_uponLight.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/5446f116-80dd-7215-a9f6-bcf336879c87.png)|![result_uponLight.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/d6e34c06-4af3-4184-0a39-4515e843212f.png)|
|:-:|:-:|
|スペクトルレンダリング|通常レンダリング|

上図では、白色のライトを拡大した画像です。`スペクトルレンダリング`では白ライトの青い周波数帯が屈折のずれで外側にずれてることが分かります。

|![result_spectrum.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/145a035e-cff5-ce07-54e2-f5e4e1a68490.jpeg)|![result_sphere.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/07b07c72-67eb-eff6-9825-ba42e9f1dd79.png)|
|:-:|:-:|
|スペクトルレンダリング|通常レンダリング|
上図では、ガラスオブジェクトを拡大した画像です。ここでも色収差が確認できます。

# インポータンスサンプリング
前の章でモンテカルロ法を用いてXYZの期待値を求めることができるようになりました。しかし、完全ランダムで任意の周波数をサンプルするよりも等色関数$\bar{x},\bar{y},\bar{z}$に従ってサンプルした方が収束が早くなります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/6c4b3a6b-4b69-9cbc-bde1-090319ef1e79.png)
上の図は等色関数$\bar{x},\bar{y},\bar{z}$です。
例えば、610nmと700nmでのサンプルはXYZの期待値の影響具合が違います。影響具合が大きい周波数を優先してサンプルすると収束が早くなります。この方法をインポータンスサンプリングと言います。

## 適切なPDF
適切なPDF用意して分散がより小さい式を採用します。
$X,Y,Z$に対して適切なPDFかを評価してみます。

$p_1(\lambda) \propto 1$
$p_2(\lambda) \propto \bar{x}(\lambda) \bar{y}(\lambda) \bar{z}(\lambda)$
$p_3(\lambda) \propto \bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda)$
上は比較するPDFになります。
(これよりもっと分散が小さくできるPDFがあるかもしれませんので探してみるのもいいかもしれません。)

$$X_{mean} = \int_{vis} \bar{x}(\lambda) d\lambda$$
$$Var(X)^2 = \int_{vis} (\lambda-X_{mean})^2 \bar{x}(\lambda) d\lambda$$
$$= \int_{vis} (\lambda-X_{mean})^2 \frac{\bar{x}(\lambda)}{p(\lambda)} p(\lambda) d\lambda$$
$P(\lambda) = \int_{380}^{\lambda} p(\lambda') d\lambda'$
$P$はCDFです。
$\lambda_i = P^{-1}(\xi) = 400 \xi + 380$
$\xi$は[0, 1)の一様乱数として、サンプル点$\lambda_i$を得ることができます。
$$Var(X)^2 = \frac{1}{N} \sum_i^N (\lambda_i-X_{mean})^2 \frac{\bar{x}(\lambda_i)}{p(\lambda_i)}$$

上の式より各$X$の分散を算出します。$Y,Z$も同様にします。

|PDF|Xの分散|Yの分散|Zの分散|
|-|-|-|-|
|$p_1(\lambda)$|132.264193|139.578475|211.571531|
|$p_2(\lambda)$|61817.248718|23832.438535|144.585084|
|$p_3(\lambda)$|73.824666 | 87.238794 | 124.222135|
※N = 1e7　の実行結果です。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/6143491d-8de4-eda6-8e8c-48901a7e7c6b.png)

数値、図から分かるように$p_3(\lambda)$が一番分散が小さいので、この式を使います。

## 乱数生成

今回は棄却法を使ってサンプル点を生成する関数を使います。
棄却法は複雑な式などでも扱えます。今回扱う$\bar{x}, \bar{y} \bar{z}$のような解析的に解くことが難しい式に対しても適応することができます。
デメリットとしては、次元が増えるとサンプル点１つを出すのに計算時間がかかる事があるので気を付けましょう。`次元の呪い`と言うようですね。今回は入力変数が1つで良かったのですが、2つ,3つ,...となると棄却法から別な方法に変える事を考える必要があるかもしれません。

```cpp
double GenerateSamplePoint() {
    // f関数に沿ったサンプル点を考える
    while(true) {
        xi_x = x_min + (x_max-x_min)*rand();
        xi_y = y_max*rand();
        if (xi_y < f(xi_x)) {
            return xi_x;
        }
    }
}
```

等色関数に沿ったサンプル点を求めることができたので、次はモンテカルロ法を用いたインポータンスサンプリングを行います。

## 導出
レイトレーシングした計算結果を$X, Y, Z$で使いまわすので、$X,Y,Z$それぞれに対して分散を小さくしてあげる必要があります。

等色関数のPDFを用意します。
PDFを$\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda)$に比例するように定義する。

$$p(\lambda) := C (\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda))$$
PDFより積分して1になるように$C$を求める

$$\int_{vis} p(\lambda) d\lambda = 1$$

$$\int_{vis} C (\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda)) d\lambda = 1$$

$$C = \frac{1}{ \int_{vis} \bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda) d\lambda }$$

定数$C$を数値計算で出しておきます。
$$\frac{1}{C} = \int_{vis} \bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda) d\lambda$$
$$p(\lambda) := \frac{1}{400}$$
$\xi_{\lambda}$は[0, 1)の一様な乱数とします。
$$\lambda_i = P^{-1}(\xi_{\lambda}) = 400 \xi_{\lambda} + 380$$
$$\frac{1}{C} = \int_{vis} \frac{\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda) }{p(\lambda)} p(\lambda) d\lambda$$
モンテカルロ法を用いて
$$\frac{1}{C} \approx \frac{1}{N} \sum_{i=1}^{N} \frac{\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda) }{p(\lambda)} $$
$$= \frac{400}{N} \sum_{i=1}^{N} ( \bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda) )$$

レンダリングする前に計算してプログラムに定数として書いておきます。($C$は定数なので、レンダリングする際に毎回計算する必要はありません。)

$\lambda_i := GenerateSamplePoint()$
$$X(x, \vec{\omega_0}) = \frac{1}{y} \int_{vis} \frac{L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{x}(\lambda)}{p(\lambda)} p(\lambda) d\lambda$$

モンテカルロ法を用いて
$$X(x, \vec{\omega_0}) \approx \frac{1}{y} \frac{1}{N} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{x}(\lambda_i)}{p(\lambda_i)} $$

$$ = \frac{1}{y} \frac{1}{N} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) \cdot \bar{x}(\lambda)}{ C(\bar{x}(\lambda) + \bar{y}(\lambda) + \bar{z}(\lambda)) } $$

$$ = \frac{1}{y N C} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) }{ 1 + \frac{\bar{y}(\lambda)}{\bar{x}(\lambda)} + \frac{\bar{z}(\lambda)}{\bar{x}(\lambda)} } $$

同様に他の$Y,Z$導出する。

$$Y(x, \vec{\omega_0}) \approx \frac{1}{y N C} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) }{ \frac{\bar{x}(\lambda)}{\bar{y}(\lambda)} + 1 + \frac{\bar{z}(\lambda)}{\bar{y}(\lambda)} } $$

$$Z(x, \vec{\omega_0}) \approx \frac{1}{y N C} \sum_{i=1}^N \frac{L_0(x, \vec{\omega_0}, \lambda_i) }{ \frac{\bar{x}(\lambda)}{\bar{z}(\lambda)} + \frac{\bar{y}(\lambda)}{\bar{z}(\lambda)} + 1 } $$

## プログラム

```cpp
struct Vec3{
    double x,y,z;
};
Vec3 xbybzbFromWavelength(double lambda) {
    Vec3 xyz;
    xyz.x = 1.056*gausse(lambda, 599.8, 37.9, 31.0) + 0.362*gausse(lambda, 442.0, 16.0, 26.7) - 0.065*gausse(lambda, 501.1, 20.4, 26.2);
    xyz.y = 0.821*gausse(lambda, 568.8, 46.9, 40.5) + 0.286*gausse(lambda, 530.9, 16.3, 31.1);
    xyz.z = 1.217*gausse(lambda, 437.0, 11.8, 36.0) + 0.681*gausse(lambda, 459.0, 26.0, 13.8);
    return xyz;
}

double sampleXYZ() {
    while (true) {
        double xi_x = (color::MAX_WAVELENGTH - color::MIN_WAVELENGTH) * rand.RandomGenerate() + color::MIN_WAVELENGTH;
        auto max = 2.1655154441358; // max value of x+y+z
        double xi_y = rand() * max;
        auto point = color::xbybzbFromWavelength(xi_x);
        double y = point.x + point.y + point.z;
        if (xi_y < y) {
            return xi_x;
        }
    }
}


Vec3 SpectrumPathTracing(Vec3 point, Vec3 direction) {
    double integral_xyz = (106.76504616237932261548607350578+106.91973463815504343620398716858+106.82532490175775630470443615231); // =int_{vis} \bar{x}(\lambda) \bar{y}(\lambda) + \bar{z}(\lambda) d\lambda
    double integral_y = 106.91973463815504343620398716858; // =int_{vis} \bar{y}(\lambda) d\lambda
    
    Vec3 spectrumRadiance;
    for (int s=0;s<samples;++s) {
        double wavelength = sampleXYZ();
        Vec3 xyz = xbybzbFromWavelength(wavelength);
        double incomingRadiance = PathTracing(wavelength);
        spectrumRadiance.x += incomingRadiance / (1.0 + xyz.y/xyz.x + xyz.z/xyz.x) / integral_y * integral_xyz / samples ;
        spectrumRadiance.y += incomingRadiance / (xyz.x/xyz.y + 1.0 + xyz.z/xyz.x) / integral_y * integral_xyz / samples;
        spectrumRadiance.x += incomingRadiance / (xyz.x/xyz.z + xyz.y/xyz.z + 1.0) / integral_y * integral_xyz / samples ;
    }
    return spectrumRadiance;
}
```

## 評価
`完全なランダムサンプル`と`インポータンスサンプリング`の2つの手法でのレンダリングと十分にサンプルしたリファレンス画像との比較します。青線は完全なランダムサンプル、赤線はインポータンスサンプリングです。

![plot_diff_by_samples.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/38281b1d-7469-a4d1-c23e-7060b9bd8aeb.png)
上図は横軸はサンプル数、縦軸は差分です。
![plot_diff_by_elapsedTime.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/e2a7bce6-22e2-f51a-a44a-2bbfa1670523.png)

上図は横軸は計算時間、縦軸は差分です。


### レンダリング結果

|スペクトルサンプル数|パストレーシングのサンプル数|
|-|-|
|5|25|

![result5.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/303476e9-6f47-7707-a894-73af0b09a2e8.jpeg)
上図：完全ランダム

![result5.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/216139/378c5d3a-cc80-4b6a-9656-5ae390f59a98.jpeg)
上図：インポータンスサンプリング

図と画像から`インポータンスサンプリング`の方が収束が早いことが分かります。`インポータンスサンプリング`は棄却法の影響で1サンプルの計算時間が少しかかることが分かりますが、それ以上に収束に早めることができています。

# git
https://github.com/Penguin77jp/RaytracingMin/tree/feature/spectrum

# 参考文献
大田登『色彩工学』、東京電機大学出版局、2001年

https://en.wikipedia.org/wiki/CIE_1931_color_space

https://rayspace.xyz/CG/contents/Spectrum_RGB/

https://raytracing.github.io/books/RayTracingInOneWeekend.html#dielectrics/snell'slaw

# まとめ
一様乱数を使ったスペクトルレンダリングまではそんなに大変ではなかったのですが、インポータンスサンプリングが少し厄介でした。等色関数が複雑な式で構成されているので逆関数法ができず色々悩んでいました。そこで、初めて棄却法をレンダリングに取り入れたのですが、とても分かりやすくシンプルで強力な手法だと思いました。棄却法はちゃんと短所もあり、棄却法と逆関数法は両極端ですね。
## 今後
今回はガラスの屈折具合をあり得ないくらい分散させたので、次は現実的な屈折ができるように[Cauchy's equation](https://en.wikipedia.org/wiki/Cauchy%27s_equation)を使いたいと思います。
今回のカメラレンズは単純な凸レンズ(虫眼鏡)でレンダリングしたので、カメラレンズの構成についてもやってみたいなと思います。
色々やることが多そうですが、色収差はとても綺麗ですね！(*^▽^*)

<font color="LightGray">(最近、デノイザーとか面白そうだなーと浮気をしそうです。)</font>