このブログポストはレイトレ合宿9アドベントカレンダーの1日目の投稿です。

はじめに

このブログポストでは、MCレイトレーシングの"Spectral Rendering"がなぜ重要であるのか、また、Spectralなデータが手元にない場合でも"Spectral Rendering"を行う方法について説明します。

下準備

今後、色に関する処理を多く行うため、Pythonの便利な色に関するライブラリであるcolour-scienceライブラリをインストールします。

In [1]:
!pip install colour-science -q

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import *
import colour

with colour.utilities.suppress_warnings(
    colour_runtime_warnings = True,
    colour_usage_warnings = True,
    colour_warnings = True,
    python_warnings = True ):
    from colour import *
    from colour.plotting import *
    from colour.models import *
    from colour.recovery import *
    from colour.colorimetry import *
    from colour.utilities import *

     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 2.3/2.3 MB 61.3 MB/s eta 0:00:00

RGBレンダリング

まずはRGBレンダリングについて説明します。特にSpectralを考慮しないレンダラーでは、次のようにRGBで計算が進みます。

vec3 lightColor = calcLightColor(hitPoint,ns)
vec3 refColor   = calcSurfaceRef(wo,wi)
pixelColor = lightColor * baseColor

このRGBの世界では、「白い光源から出た光がシアン色（=青緑）の壁に当たり、その反射光で紫色(赤青)の床が青色になる」ことになります。この表現は非常に直観的で、一見間違いがないように思えますが、実はこれは完全には正確ではありません。この理由を説明するために、まずスペクトル分布について説明します。

スペクトル分布(Spectrum Distribution, SD)

人間は特定の波長の電磁波を光として感じ、電磁波の波長ごとの強度の違いを色として認識します。しかしながら、人間が直接に電磁波の強度の違いを感じ取っているわけではなく、それを三つの感覚的な強度まで情報を減らしています。これを三刺激値と呼びます。例えば、カラーチェッカーのパッチの反射スペクトルは次のような分布となりますが、これを人間が直接そのまま認識するわけではないのです。


Image source: ColorChecker (Wikipedia)

In [2]:
with colour.utilities.suppress_warnings(colour_runtime_warnings = True, colour_usage_warnings = True, colour_warnings = True, python_warnings = True ):
    plot_multi_sds(
        [SDS_COLOURCHECKERS['ColorChecker N Ohta'][color].copy()
        for color in ['red', 'orange', 'purple', 'green', 'cyan', 'magenta']])
    plt.show()

上述のプロットは反射の強度の分布で、これを反射スペクトルと呼びます。横軸が波長で、縦軸がその波長の強度を示します。このような波長毎の強度の分布を一般にスペクトル分布(Spectrum Distribution, SD)と呼びます。反射スペクトルの場合、この強度は0から1の範囲に収まります。反射されるエネルギーが入射エネルギーを超えることはないという直感から、この制約は納得しやすいでしょう。

この反射スペクトルを人間が感じる色に変換するためには、使用している光源の情報と、どの波長がどの程度人間の感覚に感じやすいかを示す関数が必要となります。後者の関数を等色関数(Color Matching Function, CMF)と呼びます。等色関数はSDを直接RGBに変換するものもありますが、一度XYZ色空間に変換したのち、色変換行列を使用してsRGBやRec. 2020など、出力するための色空間に変換することが一般的です。

In [3]:
plot_single_cmfs("CIE 1931 2 Degree Standard Observer")
plt.show()

Spectral Renderingの必要性

現実の世界における光は電磁波であり、人間がそれを直接知覚することはできませんが、それが一体CGのレンダリングとどのように関連するのでしょうか？ CGのレンダリングには次の二つの観点から関連します。

  1. 波長に依存する現象をシミュレーションできない
  2. 大域照明や多重散乱の計算結果が変化する

これらについて一つずつ説明します。

1. 波長に依存する現象を直接取り扱えない

現実世界では、波長に依存する現象が多く存在します。例えば、ガラスのプリズムによる分光現象は、ガラスとの界面での屈折が波長に依存するために起こるのです。


Image source: Dispersive prism (Wikipedia)

これをシミュレーションするためには、レンダリングのコード内に波長が記述されていなければなりません。

# 波長をサンプル
lambda = sampleLambda()
# BRDFが波長依存
r = brdf(wo,wi,lambda)

他にも構造色、薄膜干渉、遊色、回折など現実には波長に依存する現象が数多くあります。これらを再現するためには、Spectral Renderingが必要となるでしょう。

2. 大域照明や多重散乱の計算結果が変化する

先に挙げた現象は確かに現実に存在するものの、日常で頻繁に目にする現象というよりは、特殊な現象と感じられるかもしれません。レンダリングするシーン中にそのような現象が存在しない場合、Spectral Renderingを行う必要がないように思えます。

しかし、Spectral Renderingを行わないと、現実により近いモデルからの乖離が大きくなり、MCレイトレーシングにとって大きな問題が生じます。それが大域照明(Global Illumination, GI)や多重散乱の計算結果が変わる ということです。

例えば、「白い光源から出た光がシアン色（=青緑）の壁に当たり、その反射光で紫色(赤青)の床が青色になる」という解釈は大まかには合っていますが、厳密には一致しません。「シアン色の壁」「紫色の床」という色の情報だけでは、床がどのような色になるかは厳密には判断できません。

RGBレンダリングでは大域照明や多重散乱の計算において、三刺激値それぞれを反射を繰り返す度に乗算していきます。一方、Spectral Renderingでは反射を繰り返す度にスペクトル分布(Spectrum Distribution)の乗算を繰り返し、最終的に目に入る際に三刺激値に変換します。つまり、RGBレンダリングでは、最初に畳み込まれて三刺激値になったものを乗算するのに対し、Spectral Renderingではスペクトル分布を乗算した結果が最終的に畳み込まれます。一般的に「乗算の畳み込み」と「畳み込みの乗算」は異なる結果をもたらします。

sd   = SpectrumDistribution('cyan')
assert(sd_to_xyz(sd * sd) == sd_to_xyz(sd) * sd_to_xyz(sd)) # 失敗！

百聞は一見に如かず、具体的にこの計算がどれほど変わってしまうのか、計算してみましょう。

Spectrum Distribution -> sRGB

まずは下準備から始めます。Spectrum DistributionをsRGBに変換します。前に書いた通り、この変換は次の手順で行います。

  1. スペクトルに等色関数をかけて、XYZ色空間に変換
  2. XYZ色空間からsRGB色空間への色変換行列を掛ける

これを実装してみましょう。

In [4]:
def sd_to_XYZ_(sd):
    illuminant = colour.SDS_ILLUMINANTS["D65"].copy()
    return colour.colorimetry.sd_to_XYZ_integration(sd, illuminant=illuminant)/100.0

def sRGB_to_XYZ_(sRGB): return colour.sRGB_to_XYZ(sRGB)
def XYZ_to_sRGB_(XYZ): return colour.XYZ_to_sRGB(XYZ)
def sd_to_sRGB_(sd): return XYZ_to_sRGB_(sd_to_XYZ_(sd))

# カラーチェッカーの"red"
sd  = colour.SDS_COLOURCHECKERS['ColorChecker N Ohta']['red'].copy()
rgb = sd_to_sRGB_(sd)
print(rgb)

[ 0.69868675  0.18455924  0.22679089]

"red"のパッチが赤っぽく表示されているので、成功したようです。 次に、大域照明を実装してみましょう。

大域照明

XYZやsRGBへの変換が可能になったので、次に大域照明を試してみます。レンダラーに組み込むのではなく、白色光源下で、その色の反射を繰り返すケースを考えて、大域照明の結果がどうなるのかを見てみましょう。次に二つの方法が考えられます。

  1. SDのまま乗算を繰り返して、最後にsRGBに変換する。
  2. SDをsRGBへ変換した後に、何度もその色空間で乗算する

1の方がより現実世界に近いモデル化です。これをコードにしてみましょう。

In [5]:
#
def drawColors(colorsList, names, name):
    numColors = len(colorsList)
    numBounces = len(colorsList[0])
    h = 1.0/numColors
    w = 1.0/numBounces
    fig, ax = plt.subplots(figsize=(numBounces*1.0, numColors*1.0))
    for ci, colors in enumerate(colorsList):
        cii = ci #numColors - ci - 1
        #print(numColors,ci,cii)
        plt.text(0.005, h*3/4+h*cii, names[ci],c='white')
        for idx, color in enumerate(colors):
            color = np.clip(np.array(color),0.0, 1.0)
            rect = patches.Rectangle((w*idx, h*cii), w*(idx+1), h*(cii+1), facecolor=color)
            ax.add_patch(rect)
    ax.axis('off')
    plt.title(name)
    plt.show()

sd = SDS_COLOURCHECKERS['ColorChecker N Ohta']['blue flower'].copy()
colorsList = [[] for _ in range(2)]
for bounce in range(1,12):
    colorsList[1].append(sd_to_sRGB_(sd ** bounce))
    colorsList[0].append(eotf_sRGB(eotf_inverse_sRGB(sd_to_sRGB_(sd)) ** bounce))
drawColors(colorsList, ['sRGB','Spectrum'], sd.name)

それぞれの行が乗算を行う空間で、右に行くほど反射回数が増えていきます。 理想的にはSpectrumの行とsRGBの行が一致するはずですが、一番最初のパッチ（直接光に照らされている場合）を除き、すべて異なる色になっており、反射するたびに乖離が大きくなっています。

ここではsRGBをいわゆる「リニア」にしていますが、この問題が生じていることに注意してください。また、この問題はsRGB特有のものではなく、全ての色空間で起こります。

Spectral Upsampling

大域照明を計算する際はSDを用いた方が良いことは理解できました。しかし現実的には、手元にあるテクスチャはSDではなく、三刺激値RGBの値で保存されています。もともと分布という関数だったものが、三つの実数になってしまったので、情報が大幅に減ってしまっています 。

この大幅に減った情報を何とか補い、RGBからSDを作り出すのが、Spectral Upsamplingの目指すところです。ここではSpectral Upsamplingの一手法である"A Low-Dimensional Function Space for Efficient Spectral Upsampling,"について説明します。

A Low-Dimensional Function Space for Efficient Spectral Upsampling

まずSpectral Upsamplingの要件をリストアップします。

要件1. 元のRGBに戻ってほしい

RGBからアップサンプルしたSDを再度RGBにしたとき、元のRGBに戻ってほしい。完全には戻らなくても、その誤差はできるだけ小さくなってほしい。

要件2. 滑らかなスペクトルになってほしい

RGBからアップサンプルしたSDの形状は、滑らかな関数になっていることが望ましい。なぜなら、先に見たカラーチェッカーのように、現実世界の反射スペクトルは滑らかに変化する傾向があるからです。ただし、この論文では反射スペクトルのアップサンプルだけを対象にしていますが、発光のスペクトルである分光スペクトルについては、必ずしも滑らかであるとは限りません。例えばナトリウムランプは鋭いピークを持ちます。


Image source: Sodium-vapor lamp (Wikipedia)

では、これらの要件を満たす方法を見ていきましょう。

二次関数とシグモイド関数

次のような変換を行います。

$$ f(\lambda) = S(c_0 \lambda^2 + c_1 \lambda + c_2) \\ $$

ここで、$c0,c1,c2$は係数で、$\lambda$は波長を表します。$S$はシグモイド関数で、atan、双曲線関数、ロジスティック関数など何でもいいとされていますが、実行時の計算コストを考えて、以下の形のシグモイド関数が採用されています。

$$ S(x) = \frac{1}{2} + \frac{x}{2 \sqrt{1+x^2}} $$

これは単純な関数の組み合わせとなりますが、その意味を以下で説明します。

二次関数は山、谷、そして平坦な関数を表現でき、これにより大体の反射スペクトルは表現できるとされています。一方、シグモイド関数は定義域が$[-\infty,\infty]$、値域が$[0,1]$の関数です。反射スペクトルは$[0,1]$の範囲に収まってほしいですが、二次関数の値域は何にでもなり得ます。シグモイド関数を通すことで、強制的に$[0,1]$の範囲に制限することができます。二次関数を上下から押し潰すイメージです。二次関数もシグモイド関数も滑らかな関数なので、滑らかな関数になってほしいという要件を満たします。

係数$c_0,c_1,c_2$を求める

次に係数$c_0,c_1,c_2$の求め方ですが、これは以下のような非線形最適化問題となります。

$$ \hat{f}=\operatorname{argmin}\left\|b-T \int_{\Lambda} f(\lambda) W(\lambda) x y z(\lambda) d \lambda\right\| $$

ここで、$b$はRGBの値、$T$は色変換行列、$W$はD65のスペクトルエネルギー分布、$xyz$はCMFを表します。つまり、「D65光源で照らしたときの全波長のアップサンプルが、もっともRGBに戻したときの誤差の合計が小さい係数を選ぶ」ということになります。

この計算によって求めた係数と、いくつかの工夫により、64x64x64の係数のルックアップテーブルによって、sRGB色空間においては、RGB->SD->RGBとRGBに戻したときに、誤差が0になると論文では述べられています。これは誤差が少ない方がいいという要件1を満たします。

使ってみる

color scienceライブラリにはこのSpectral Upsampling手法が実装されているので、何も考えずにSpectral Upsamplingを実行することができます。実際に使ってみましょう。

In [6]:
# sRGBからSDへのアップサンプル
def upsample_sRGB_to_sd_(sRGB):
    return XYZ_to_sd(sRGB_to_XYZ_(sRGB), method="Jakob 2019")

# 試しにスペクトルをプロット
sd = colour.SDS_COLOURCHECKERS['ColorChecker N Ohta']['green'].copy()
sdUpsamp = upsample_sRGB_to_sd_(sd_to_sRGB_(sd))
sdUpsamp.display_name = "Upsample"
plot_multi_sds([sd,sdUpsamp])
plt.show()

アップサンプルしたSDが元のSDにそれなりにフィットしている様子がわかります。ただし注意すべき点としては、反射スペクトルに直接フィットしたわけではなく、これにD65光源とXYZ等色関数を掛けた後の世界でフィットしているので、このグラフはあまりフィットしているように見えないかもしれません。

最後に、アップサンプルの結果を使って大域照明の計算をしてみましょう。この計算ではアップサンプル後のSDを使用して乗算を行い、その結果をsRGBに変換します。

In [7]:
def plotComparison(sd):
    colorsList = [[] for _ in range(3)]
    for bounce in range(1,12):
        colorsList[2].append(sd_to_sRGB_(sd ** bounce))
        colorsList[1].append(sd_to_sRGB_(upsample_sRGB_to_sd_(sd_to_sRGB_(sd)) ** bounce))
        colorsList[0].append(eotf_sRGB(eotf_inverse_sRGB(sd_to_sRGB_(sd)) ** bounce))
    drawColors(colorsList, ['sRGB','Upsample','Spectrum'], sd.name)

plotComparison(colour.SDS_COLOURCHECKERS['ColorChecker N Ohta']['blue flower'].copy())

結果の乖離が明らかに減少したことが分かります。完全には一致していませんが、sRGBに比べて、値がSDにかなり近く保持されていることが確認できます。

改めて述べますが、これは本来知りうるSDを推測するというよりも、それらしい形状のSDをRGBから求める手段です。情報の乏しいRGBから出発している以上、完璧なSDを得ることは不可能です。しかしながら、それにもかかわらず、それらしいSDが得られているという事実は驚異的です。

まとめ

このブログポストでは、Spectral Upsamplingとその一手法について説明しました。特に、sRGBからSDへの変換とその後の大域照明計算の挙動について、みてきました。結果として、SDへのアップサンプリングを行うことで、反射を繰り返しても、sRGBに比べてSDにかなり近い値を維持できることが確認できました。完全に一致するまでは至っていませんが、明らかに結果の乖離が減少しました。

自分のレンダラーにもぜひSpectral Upsamplingを組み込んでみましょう!

参考文献

[1]W. Jakob and J. Hanika, “A Low-Dimensional Function Space for Efficient Spectral Upsampling,” Computer Graphics Forum (Proceedings of Eurographics), vol. 38, no. 2, Mar. 2019.

おまけ1. C++コードから使えるようにする

PythonからSpectrum Upsamplingを行う方法を説明してきましたが、Colour *Scienceライブラリには、内部の係数を出力する機能が付いています。これを使って係数のテーブルを出力して、使う方法があります。

In [8]:
LUT = LUT3D_Jakob2019()
LUT.generate(RGB_COLOURSPACE_sRGB, None, None, 3, lambda x: x) # 3を10にするとすごい時間がかかる...
LUT.write('sRGB.conf')

100%|██████████| 27/27 [00:23<00:00,  1.17it/s]
Out[8]:
True

またプロジェクトページのSupplemental materialに、係数テーブルとそれを使うC/C++のコードが含まれているのでこれを使うのがいいと思います。

おまけ2. カラーチェッカーの全ての色での比較

In [9]:
for k,sd in SDS_COLOURCHECKERS['ColorChecker N Ohta'].items():
    # sRGBのcyanが壊れるのでとりあえずスキップ
    if k == 'cyan': continue
    plotComparison(sd.copy())
In [ ]:
# ドライブをマウントし、Notebookがあるフォルダまで移動する
import google.colab as colab
colab.drive.mount('/content/drive/')
%cd "/content/drive/MyDrive/Colab Notebooks/"

# htmlに変換してダウンロード
!jupyter nbconvert --to html "spectral_upsampling.ipynb" --template=classic
from google.colab import files
files.download('spectral_upsampling.html')