このブログポストはレイトレアドベントカレンダー2022の19日目の投稿です。
はじめに¶
テクスチャをGL_REPEATなどでタイリングすると、リピート感が出てしまい、とてもチープな見た目になってしまいます。
この問題を解決するために、入力画像と似たようなパターンを無限に作り出し、それをスムースに繋ぐという手法[1][2][3][4]があります。[1][2]の手法の核である「似たようなパターンを無限に作る」の部分は、手法の内容を知らないと魔法のように感じる部分ですが、根底にあるアイデアはとてもシンプルです。この投稿では、この核となるアイデアについて解説していきます。
それではまずは前準備からはじめてみましょう！
前準備¶
まずDr.Jitをインストールしします。 Dr.Jitとは、PythonコードをJITコンパイルしてSIMDやCudaで動くようにしたり、自動微分ができるライブラリです。微分可能レンダラーのMitsuba3のバックエンドとしても使われています。今回は自動微分機能は使わず、JITコンパイルをした上でCuda上で高速に動かすために使用します。Dr.Jitの詳細はgithubや、リファレンスを参照してください。このJupyterのNotebookはGoogle Colab上で作られているのですが、Colabの拡張機能としてターミナルのコマンドをNobtebook上で実行できるので、それを利用してpipを使ってDr.Jitをインストールします。
!pip install drjit --quiet
次にPOT(Python Optimal Transport)をインストールします。これは最適輸送(Optimal Transport, OT))を解くためのライブラリです。詳細は後述します。
!pip install POT --quiet
無事インストールできました。
次に必要なライブラリをimportしていきます。
Dr.JitとPOT以外にもnumpyやmatploblibなどもimportしていきます。
# 定番ライブラリ
import math
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
from scipy import stats
# Dr.JIT
import drjit as dr
from drjit.cuda import Array2f, Array3f, Texture2f, Texture3f, Float, TensorXf
#from drjit.llvm import Array2f, Array3f, Texture2f, Texture3f, Float, TensorXf
# POT
import ot
import ot.plot
各種ライブラリをimportできました。
次に便利な機能を先に作っておきます。画像処理をShaderToyのようにフラグメントシェーダー的に書くための便利な関数も作ります。
# ---------------------------
# 便利な定数と関数の準備
# ---------------------------
# ドライブをマウントし、Notebookがあるフォルダまで移動する
import google.colab as colab
colab.drive.mount('/content/drive/')
%cd "/content/drive/MyDrive/colab/infinite patterns"
# テクスチャのサンプル
def sample(tex,uv):
return Array3f(tex.eval_cubic(uv))
#return Array3f(tex.eval(uv))
Texture2f.sample = sample
Texture3f.sample = sample
# テクスチャの値をRGBの配列にする
def texToRGB(tex: Texture2f):
tv = tex.value().numpy()
tv = tv.reshape(tex.shape[0]*tex.shape[1], 3)
return tv
# ファイル名からテクスチャを生成
def loadTexture(filename : str) -> Texture2f:
img = np.asarray(Image.open(filename))
# 4チャンネル画像も3チャンネルに変換
img = img[:,:,:3]
img = img / 256.0
tex = Texture2f(TensorXf(img))
return tex
# DrJitにfmodがないので追加
def fmod(a,b):
return a - dr.floor(a/b) * b
dr.fmod = fmod
# shadertoyのように画像処理を書けるようにする関数
def shadertoy(f, width):
# UVの生成
x = dr.linspace(Float, 0.0, 1.0, width)
y = dr.linspace(Float, 0.0, 1.0, width)
u, v = dr.meshgrid(x, y)
uv = Array2f(u, v)
# レンダリングする関数を呼ぶ
img = f(uv)
# テンソル型にして表示
imgt = TensorXf(dr.ravel(img), shape=(width, width,3))
dpi = 80
height_ratio = 4
figsize = width / float(dpi), (width*(1+height_ratio)/height_ratio) / float(dpi)
fig, ax = plt.subplots(2, 1, figsize=figsize,tight_layout=True, gridspec_kw={'height_ratios': [height_ratio, 1]})
ax[0].axis("off")
ax[0].imshow(imgt)
# ヒストグラムをカラーチャンネル毎に表示する
ax[1].tick_params(left=False)
ax[1].axes.yaxis.set_ticklabels([])
ax[1].set_xlim([0.0,1.0])
img_tmp = img.numpy()
ax[1].hist(img_tmp[:,0].flatten(), bins=256, histtype='step', color='red')
ax[1].hist(img_tmp[:,1].flatten(), bins=256, histtype='step', color='green')
ax[1].hist(img_tmp[:,2].flatten(), bins=256, histtype='step', color='blue')
plt.tight_layout()
plt.show()
元画像の表示¶
これで下準備が整いました。
ちゃんと想定通り動くか試すために画像を単純に貼り付けてみます。
# イメージファイルからテクスチャをロードして、そのまま貼り付ける
# 画像はPolyで生成しました
# https://withpoly.com/texture/create?asset_id=Q9NrGKoKuj
tex_example = loadTexture(f'moss_128.png')
def simpleDisplay(uv : Array2f) -> Array3f:
return tex_example.sample(uv)
shadertoy(simpleDisplay, 512)
ちゃんと想定通り動きました。
簡単にコードの説明をします。
- shadertoy関数は、画像を生成し、それを表示します。第1引数に画像全体で実行したいフラグメントシェーダーのような関数、第2引数にその画像のサイズを指定します。今回は512x512の画像を最終的に出力することになります。
- simpleDisplay関数は、UVを引数にとります。これはArray2f型で、UV座標が全ピクセル数分だけ入ってます。今回の場合は512x512=262144個のUVが格納されています。
- texは今回貼り付けるテクスチャが入っており、sample関数でそのテクスチャの指定した座標のピクセル値をサンプルして返します。
- 下に表示されているのはチャンネル毎のヒストグラムです。
変数uvや、関数tex.sampleの戻り値は全ピクセルの値が入った(Dr.Jitの)配列であることに注意してください。そのため、フラグメントシェーダーのようなピクセル単位の分岐は通常のif文ではできないことに注意してください。もう一つの例題として、試しにGL_REPEATのような絵を作ってみます。
# 単純にリピートして表示
def simpleRepeat(uv : Array2f) -> Array3f:
uv = dr.fmod(uv*4.0, 1.0)
return tex_example.sample(uv)
shadertoy(simpleRepeat, 512)
繰り返して貼り付けることに成功しました。当たり前ですが、ヒストグラムは元の画像と完全に同じになっています。
さて次は、入力画像から似たような画像を作る話の前に、ガウシアンノイズ画像から別のガウシアンノイズ画像を作る話をします。
ガウシアンノイズ¶
「入力画像と似たようなパターンを無限に作る」という問題に取り組む前に、シンプルな問題を考えてみます。具体的には、ガウシアンノイズ画像から別のガウシアンノイズ画像を作ってみます。真っ先に思い浮かぶ方法は、ガウシアンノイズ画像を足し上げていって、それを単に平均してしまう方法です。"ランダムな"画像を足しても、"ランダムな"画像になる気がします。
まずはガウシアンノイズ画像を作って表示する部分を作ります。
# ガウスノイズの画像を作成する
def createGaussianNoiseTexture(size):
# Truncated normal distribution
lower, upper = 0, 1.0
mu, sigma = 0.5, 0.25
X = stats.truncnorm((lower - mu) / sigma, (upper - mu) / sigma, loc=mu, scale=sigma)
# 1chで生成して、それを3chに複製する
img = X.rvs([size[0],size[1],1])
img = np.tile(img,(1,3))
tex = Texture2f(TensorXf(img))
return tex
# 生成したガウスノイズテクスチャを表示する
tex_gaussian_noise = createGaussianNoiseTexture((512,512))
def showGaussianNoiseTexture(uv : Array2f) -> Array3f:
return tex_gaussian_noise.sample(uv)
shadertoy(showGaussianNoiseTexture, 512)
ガウシアンノイズが生成されました。下のヒストグラムは、画像の中に登場するピクセルの値をヒストグラムにしたものです。ガウシアンなので、単峰性で、中央から離れれば離れるほど、急激にその値のピクセル数が少なくなっているのがわかります。
次にこのガウシアンを複数枚作って、平均をとります。具体的には4つのガウシアンを加算して加算した分だけ除算してみます。また同じ確率分布に従うノイズテクスチャが生成されるような気がします。
# 生成したガウシアンノイズテクスチャを表示する
texSize = (512,512)
tex_gaussian_noise0 = createGaussianNoiseTexture(texSize)
tex_gaussian_noise1 = createGaussianNoiseTexture(texSize)
tex_gaussian_noise2 = createGaussianNoiseTexture(texSize)
tex_gaussian_noise3 = createGaussianNoiseTexture(texSize)
def simpleCompositeGaussianNoiseTexture(uv : Array2f) -> Array3f:
p0 = tex_gaussian_noise0.sample(uv)
p1 = tex_gaussian_noise1.sample(uv)
p2 = tex_gaussian_noise2.sample(uv)
p3 = tex_gaussian_noise3.sample(uv)
return (p0+p1+p2+p3)/4.0
shadertoy(simpleCompositeGaussianNoiseTexture, 512)
残念ながら想定通りのノイズテクスチャになりませんでした。
どうもぼやけたような画像になってしまいました。これは目の錯覚ではありません。ヒストグラムを見ると、実際に分布のバラツキが減っています。結果見た目のコントラストも減ってしまっていたわけです。
確率分布に従う確率変数の和の確率分布の分散は、元の確率分布の分散の加算#%E6%80%A7%E8%B3%AA)になります。今回は同じ確率分布に従うので、元のガウス分布の4つ分、分散は4倍になっています。標準偏差は分散の平方根なので$\sqrt{4}=2$倍となります。今回は足した確率変数をさらにその数の4で割るので、$2/4=1/2$倍が最終的な分布の標準偏差になります。つまり元の標準偏差の半分になってしまっているわけです。確かにヒストグラムを見比べるとすそ野の幅が半分になっているように見えます。
これをどうにかしてコントラストを復活させることはできないでしょうか。
コントラストを復活させる¶
ばらつきが半分になったというのであれば、分布の中心から2倍に引き伸ばしてしてしまえば、同じ分布になってくれないかなと考えます。ヒストグラムをみると確かに確かに、分布の形状は同じだけれども、横に縮んだだけに見えるので、単に横に引き延ばすだけでうまくいきそうです。
やってみます。
# ガウスノイズの合成の分散を復活させる
def simpleCompositeGaussianNoiseTexture(uv : Array2f) -> Array3f:
p0 = tex_gaussian_noise0.sample(uv)
p1 = tex_gaussian_noise1.sample(uv)
p2 = tex_gaussian_noise2.sample(uv)
p3 = tex_gaussian_noise3.sample(uv)
sE = (p0+p1+p2+p3)/4.0
# ガウスノイズの期待値
Ex = 0.5
# 各項のウェイト
wi = 1.0/4.0
# 補正用の値
sW = math.sqrt(4*(wi**2))
# 最終的な値
ret = (sE - Ex)/sW + Ex
return ret
shadertoy(simpleCompositeGaussianNoiseTexture, 512)
うまくいってしまいました！
どうやらガウシアンノイズを足したものは、ガウシアンノイズのバラツキだけがが変化したものになるようです。しかしこれは一般的な画像でも同じことが言えるのでしょうか？
# 画像のオフセット
offsets = np.random.random((16,2))
# 画像の平均値を出す
def texMu(tex):
rgb = texToRGB(tex)
r,g,b = rgb[:,0], rgb[:,1], rgb[:,2]
return np.array([np.mean(r),np.mean(g),np.mean(b)])
mu = texMu(tex_example)
# UVを変えて数回合成する
def composeMultipleView(uv : Array2f) -> Array3f:
sum = 0.0
for offset in offsets:
nuv = dr.fmod(uv+offset, 1.0)
sum += tex_example.sample(nuv)
sum /= len(offsets)
# 各項のウェイト
wi = 1.0/len(offsets)
# 補正用の値
deno = math.sqrt(4*(wi**2))
# 最終的な値
ret = (sum - mu)/deno + mu
return ret
shadertoy(composeMultipleView, 512)
元画像と似ても似つかないものになってしまいました。さらにヒストグラムも元の画像と全く違い、ガウス分布のようになってしまっています。(「分布を足していけばどんどんガウス分布になる」という中心極限定理の効果がまさに出てしまっています)
これはテスト画像の色の分布が、ガウス分布のように再生成が成立する分布ではなかったためです。もし画像を何らかの方法で、合成の瞬間だけガウス分布のようにすることができればヒストグラムを変形させることがなくなりそうです。
実はこれが可能です。そのような変換を行うために、最適輸送という道具を使います。
# 砂山を移動する例
# 公式サンプルの"Optimal Transport for 1D distributions"を使わせていただきました。
# https://pythonot.github.io/auto_examples/plot_OT_1D.html#sphx-glr-auto-examples-plot-ot-1d-py
from ot.datasets import make_1D_gauss
# 砂の数は100個
n = 100
x = np.arange(n, dtype=np.float64)
# 一つ目の砂山(青)
a = ot.datasets.make_1D_gauss(n, m=20, s=5)
# 二つ目の砂山(赤)
b = ot.datasets.make_1D_gauss(n, m=60, s=10)
plt.figure(1, figsize=(6.4, 3))
plt.plot(x, a, 'b', label='source')
plt.plot(x, b, 'r', label='target')
plt.legend()
plt.show()
このグラフは、高低で砂の量を表した一次元の砂山を表しています。青い山を赤の山にするのが目標です。青い山を構成する砂を右の赤い山に移動していくわけです。そして、どこの砂をどこに移動させるのが、最もトータルの輸送コストが低いかということを解きます。
# コスト行列を作る
# 今回は全ての砂は単純に距離に比例するコスト
M = ot.dist(x.reshape((n, 1)), x.reshape((n, 1)))
M /= M.max()
# POTのソルバーを使って最適な輸送順番を導き表示する
G0 = ot.emd(a, b, M)
plt.figure(3, figsize=(5, 5))
ot.plot.plot1D_mat(a, b, G0, 'OT matrix G0')
plt.show()
解けました。
グラフの読み方としては、左の青のプロットと、上の赤のプロットの交差するところが黄色になっている部分が、移動前後で同じ砂が存在する位置です。この例は単純な1Dでかつ、全ての砂の輸送コストが全て同じ条件なので、単に右から順番に運んでいけばいいということが示されただけですが、次元が上がったり、輸送コストが砂によって違う場合などは問題の難易度が格段に上がります。これから解こうとしている問題はRGBの三チャンネルの問題なので、三次元の最適輸送の問題であり、難しい問題です。しかしPOTを使えば、簡単に解けます。
画像のガウス分布化¶
さっそくPOTを使って入力画像を、再生性があるガウス分布に変換してみましょう！
# テクスチャのピクセルをランダムにサンプルする
# 全テクセルに対して最適化をしたくなるが、
# 計算量的にできないので、いくつかのテクセルをサンプルする
rgb = texToRGB(tex_example)
nb = 2000
rnd = np.random.RandomState(0x1234)
idx = rnd.randint(rgb.shape[0], size=(nb,))
Xs = rgb[idx, :]
# 三次元ガウシアンからサンプルする
# 各次元は無相関なので、共分散行列は対角行列
# 分散の大きさは1/36は、論文中の数字をそのまま持ってきた
# Xtが範囲の外に出たら単にクリップ
mu = np.array([0.5, 0.5, 0.5])
cov = np.diag(np.ones(3)/36)
Xt = np.random.multivariate_normal(mean=mu, cov=cov, size=nb)
Xt = np.clip(Xt, 0.0, 1.0)
# EMDで最適化
transport = ot.da.EMDLaplaceTransport() # ot.da.EMDTransport()
transport.fit(Xs=Xs, Xt=Xt)
最適化が終わりました。最適化の結果を用いて、サンプルが実際にどこに輸送されたのかプロットします。先ほどの例は一次元でしたが、今回は三次元の問題です。三次元をプロットすると見づらいので、最初の二次元だけをプロットします。
# どう転送される行列が作られたかプロット
fig = plt.figure(figsize = (8,8))
ot.plot.plot2D_samples_mat(Xs[:,:2], Xt[:,:2], transport.coupling_, c=[.5, .5, 1])
plt.xlim(0.0,1.0)
plt.ylim(0.0,1.0)
plt.plot(Xs[:, 0], Xs[:, 1], '+b', label='source')
plt.plot(Xt[:, 0], Xt[:, 1], 'xr', label='target')
plt.legend()
plt.show()
青の点が元の画像の上のピクセルの値で、赤の点が多次元のガウス分布です。その間にある青の線が、どの点とどの点が対応するかを示しています。注意点としては、同じ赤の点には一つの赤の点しか割り当てできないので、必ずしも青い点に一番近い赤い点には割り当てられないことです。お互いが譲り合って、全体として移動コストが最も小さくなるように配分されます。プロットをみるとうまく行っていそうですが、かなり複雑な解であることも伝わってきます。
ではこの解を用いて入力画像をガウス分布にしてみます。
# ガウス分布化画像のテクスチャを作成
flatten_gaussianized_img = transport.transform(Xs=rgb)
gaussianized_tex = Texture2f(TensorXf(flatten_gaussianized_img.reshape(tex_example.shape)))
# ガウス化画像をテクスチャにして表示する
def showGaussianized(uv : Array2f) -> Array3f:
return gaussianized_tex.sample(uv)
shadertoy(showGaussianized, 512)
ガウス化された画像をみてもうまくいっているのかよくわかりませんが、ヒストグラムをみると、ガウス分布っぽくなっているのでおおよそ成功しているようです。
※ RGBの全チャンネルがぴったり一致した綺麗なヒストグラムになるのが理想なのですが、最適輸送のアルゴリズムの計算量がO(n^3)であるため、画像のピクセルをランダムサンプルするしかなく、その結果綺麗なガウス分布になっていません。改良手法[2]では3次元の最適輸送をせず、全ピクセルをみても十分高速な方法になったので、この問題が緩和されました。今回は原理の説明なので「大体ガウス分布になった」と目をつぶってください。
逆方向の変換も行って、入力画像にちゃんと戻るかを確認します。
# ガウス分布化画像を元に戻したテクスチャ(=元に戻って欲しいテクスチャ)を作成
degaussianized_img = transport.inverse_transform(Xt=flatten_gaussianized_img)
degaussianized_tex = Texture2f(TensorXf(degaussianized_img.reshape(tex_example.shape)))
# ガウス化を元に戻したテクスチャを表示する
def showDegaussianized(uv : Array2f) -> Array3f:
return degaussianized_tex.sample(uv)
shadertoy(showDegaussianized, 512)
ちゃんと元の画像に戻っています。
では、いよいよ、ガウス化した画像を使って合成して、別のパターンを作れるかを試していきます。例題なので、画像全体に対してUVのオフセットをいくつか作り、それらのUVでガウス化画像を平均して、新しい画像を作ります。
# 画像のオフセット
offsets = np.random.random((4,2))
# 画像の平均値。理想的には0.5になっているはずなので本来は不要。
tex_mu = texMu(gaussianized_tex)
# UVを変えて数回合成する
def composeMultipleView(uv : Array2f) -> Array3f:
sum = 0.0
for offset in offsets:
nuv = dr.fmod(uv+offset, 1.0)
sum += gaussianized_tex.sample(nuv)
sum /= len(offsets)
ns = len(offsets)
# 各項のウェイト
wi = 1.0/ns
# 補正用の値
deno = math.sqrt(ns*(wi**2))
# ガウス分布を合成した最終結果
gaussian_result = (sum - mu)/deno + mu
# 元画像のヒストグラムにする
return TensorXf(transport.inverse_transform(Xt=gaussian_result.numpy()))
shadertoy(composeMultipleView, 512)
ついにできました！ 入力画像と似たような見た目でかつ、別のパターンにちゃんとなっています。
話していないことリスト¶
核になるアイデアだけを話しましたが、これだけでは実用にはなりません。主に次のトピックについて話していません。
- 無限にパターンを作っただけでは一枚の画像にしかなりません。無限に連続していくテクスチャを作るためには、それをスムースに繋ぐための方法が必要です。テクスチャ空間で領域を三角形に分割して、ブレンドする率を重心座標に従って滑らかに変えることで、無限に連続していくテクスチャを実現します。
- ガウス分布の画像を入力画像のヒストグラムに戻す時に、POTを直接使いましたが、ゲームやレンダリングでは、ランタイムで行う部分なので高速に動作しなければなりません。[1]では$32^3$キューブテクスチャをLUTとして使うことで高速化を図っています。
- 最適輸送はとても重いアルゴリズムなので、全てのピクセルではなく間引いたピクセルの上でだけ計算しました。その結果、ガウス分布化した画像のヒストグラムも少し歪んでしまいます。[2]では完全な最適輸送をせずに、主成分分析のように、RGBの固有ベクトルを求めたのちに、その上で各次元の1次元での最適輸送をしています。1次元なので全ピクセルを対象にしても高速に求まります。デモページもあります。
おわりに¶
昨年のアドベントカレンダーに続き、Heitz氏の論文の実装を行なってみました。簡単な計算でこの結果が出るのは面白いですね。ヒストグラムが保存されるとうれしいのはわかるのですが、それだけでここまで見た目がいい感じになるのが意外でした。
ぜひ自前のレイトレーシングレンダラーにも組み込んでみましょう！
参照¶
[1]E. Heitz and F. Neyret, “High-Performance By-Example Noise using a Histogram-Preserving Blending Operator,” Proceedings of the ACM on Computer Graphics and Interactive Techniques, vol. 1, no. 2, p. Article No. 31:1-25, Aug. 2018, doi: 10.1145/3233304.
[2]T. Deliot and E. Heitz, “Procedural Stochastic Textures by Tiling and Blending,” GPU Zen 2
[3]M. S. Mikkelsen, “Practical Real-Time Hex-Tiling,” Journal of Computer Graphics Techniques (JCGT), vol. 11, no. 3, pp. 77–94, Aug. 2022, [Online]. Available: http://jcgt.org/published/0011/03/05/
[4]B. Burley, “On Histogram-Preserving Blending for Randomized Texture Tiling,” Journal of Computer Graphics Techniques (JCGT), vol. 8, no. 4, pp. 31–53, Nov. 2019, [Online]. Available: http://jcgt.org/published/0008/04/02/
# htmlに変換
!jupyter nbconvert --to html "infinite patterns.ipynb"