執筆者：廣瀬
映像業界を経てモノリスソフトへ入社。 以来、テクニカルアーティストとして主にエフェクト関連の業務を担当。 好きな食べものはソフトクリーム。
TECH BLOG
Height MapとNormal Mapの相互変換
はじめに
こんにちは。モノリスソフト テクニカルアーティストの廣瀬です。今回はCGでよく利用される、Height Map と Normal Map について、それらを相互変換する手法について解説していきたいと思います。
初めに、Height Map から Normal Map を作成する方法を解説します。Normal Map は、単純な計算を行うだけで作成することができます。
次に、Normal Map から Height Map を作成する方法を解説します。この Normal Map からの Height Map の作成は、難しい問題なため、いくつかの手法が考案されています。
今回は、
- 線積分
- ポアソン方程式
- 最小化問題
による作成手法ついて解説します。
上記の各手法に関しては、NumPy と SciPy による実装コードも記載しておきます。
Height Map → Normal Map
検証にはこちらの Height Map を使用しています。
まず、Height Map から Normal Map を作成する方法を解説します。
Normal Map は座標の取り方で正負の方向が逆転することに注意してください。（よく言われる OpenGL形式 と DirectX形式）
Height Map を z(x,y)、Normal Map を 3次元の単位ベクトル n(x,y) とします。
Height Map z(x,y) に対して x軸、y軸方向に偏微分を行います。ここでは、z(x,y) の中心差分を取ることとします。
\(
\displaystyle
\begin{eqnarray}
\frac{\partial z}{\partial x} &=& \frac{z(x+1,y) - z(x-1,y)}{2\Delta x}, \\
\frac{\partial z}{\partial y} &=& \frac{z(x,y+1) - z(x,y-1)}{2\Delta y}. \hspace{60pt} (1)
\end{eqnarray}
\)
式 (1) の偏微分の値から x-z軸方向のベクトル x(x,y)、y-z軸方向のベクトル y(x,y) を作成します。
\(
\displaystyle
\begin{eqnarray}
\vec{x}(x,y) &=& (\Delta x, 0, \frac{\partial z}{\partial x}), \\
\vec{y}(x,y) &=& (0, \Delta y, \frac{\partial z}{\partial y}). \hspace{60pt} (2)
\end{eqnarray}
\)
式 (2) のベクトル同士で外積を取り、それを正規化したのが Normal Map のベクトル n(x,y) です。
\(
\displaystyle
\vec{n}(x,y) = \frac{\vec{x}(x,y) \times \vec{y}(x,y)}{\| \vec{x}(x,y) \times \vec{y}(x,y) \|}. \hspace{60pt} (3)
\)
単位ベクトル n(x,y) の各成分は -1 ~ 1 の範囲になるため、一般的にテクスチャとして保存する際は、各成分に * 0.5 + 0.5 を行い 0 ~ 1 の範囲にリマップした値が使用されます。
【左】Height Map【右】作成されたNormal Map
ソースコード
import os import numpy as np import matplotlib.image as mpimg import matplotlib.pyplot as plt filename = 'C:/Users/%USER%/Downloads/Approximate_Earth_Heigh_Map.png' reverse = True # テクスチャカラーGの反転 OpenGL or DirectX img = mpimg.imread(filename, format='png') img = img[:,:,:3] # アルファチャンネルを削除 ny, nx, _ = img.shape dy, dx = np.gradient(img[:,:,0], 2.0/ny, 2.0/nx) if reverse: dy *= -1 norm = np.sqrt(dx*dx + dy*dy + 1.0*1.0) img[:,:,0] = -dx / norm img[:,:,1] = -dy / norm img[:,:,2] = 1.0 / norm img = img * 0.5 + 0.5 # -1 - 1 から 0 - 1 の範囲へ変更 mpimg.imsave(os.path.splitext(filename)[0] + '_output.png', np.clip(img, 0, 1)) plt.imshow(img, interpolation='nearest') plt.show()
Normal Map → Height Map
検証にはこちらの Normal Map を使用しています。
File:Normal map example - Map.png - Wikimedia Commons
Image by Julian Herzog, licensed under CC BY 4.0
NOTE:
Normal Mapを作成した際の差分法は正確なHeight Mapを作成する上では重要な要素ではありますが、多くの場合はそういった詳細な情報はないため、数値誤差として無視してしまいます。
次に Normal Map から Height Map を作成する方法です。
Height Map を z(x,y)、Normal Map を -1から1の範囲の3次元ベクトル n(x,y) とします。
Normal Map n(x,y) から x軸方向の偏微分の値 p(x,y)、y軸方向の偏微分の値 q(x,y) を復元します。
\(
\displaystyle
\begin{eqnarray}
p(x,y) &=& -\frac{\vec{n}_{x}}{\vec{n}_{z}} \simeq \frac{\partial z}{\partial x}, \\
q(x,y) &=& -\frac{\vec{n}_{y}}{\vec{n}_{z}} \simeq \frac{\partial z}{\partial y}. \hspace{60pt} (4)
\end{eqnarray}
\)
ここで復元した p(x,y) と q(x,y) と Normal Map の作成元の Height Map の偏微分の値の間には、テクスチャの圧縮などにより誤差が含まれていることに注意してください。
この p(x,y)、q(x,y) を元に Height Mapを作成します。
線積分
最も単純な方法は、各軸に沿って線積分を行う方法です。
\(
\displaystyle
z(u,v) = \int_{0}^{v} q(0,y) dy + \int_{0}^{u} p(x,v) dx + C. \hspace{60pt} (5)
\)
C は偏微分を行った際に消えてしまった積分定数で、Height Map 全体のオフセット量になります。今回、このオフセット量は重要ではないため、無視してしまいます。
この手法は単純明快で実装も容易ですが、欠点として数値誤差による影響が大きく、作成された Height Map に縞状のノイズが乗ってしまうことがあります。
【左】式 (5) の計算結果【右】左の計算結果にガンマをかけてノイズを分かりやすく可視化した状態
画像を3D化した結果
ソースコード
import os import numpy as np import matplotlib.image as mpimg import matplotlib.pyplot as plt filename = 'C:/Users/%USER%/Downloads/Normal_map_example_-_Map.png' reverse = True # テクスチャカラーGの反転 OpenGL or DirectX img = mpimg.imread(filename, format='png') img = img[:,:,:3] # アルファチャンネルを削除 img = img * 2.0 - 1.0 # 0 - 1 から -1 - 1 の範囲へ変更 ny, nx, _ = img.shape # 式(4) p = -img[:,:,0] / img[:,:,2] q = -img[:,:,1] / img[:,:,2] if reverse: q *= -1 p *= 2.0 / nx q *= 2.0 / ny # 式(5) C = 0.0 zy = 0.0 for y in range(ny): zy += q[y,0] zx = 0.0 for x in range(nx): zx += p[y,x] img[y,x,0:3] = zy + zx + C mpimg.imsave(os.path.splitext(filename)[0] + '_output.png', np.clip(img, 0, 1)) plt.imshow(img, interpolation='nearest') plt.show()
ポアソン方程式
別の手法として、ポアソン方程式を解くことで Height Map を求める方法があります。
2次元のベクトル場 r を以下の様に定義します。
\(
\displaystyle
\textbf{r}(x,y) = (p(x,y),q(x,y)). \hspace{60pt} (6)
\)
r の発散が解となる z を求めることで Height Map を作成することができます。
\(
\displaystyle
\Delta z = - (\mathrm{div} \; \textbf{r}). \hspace{60pt} (7)
\)
このポアソン方程式は、様々な方法で解くことができます。
ポアソン方程式の解法の例については、詳しくは「Houdiniでぽあそん」の記事を参照してください。
今回は SciPy の疎行列を用いた直接法で解いています。
式 (7) の計算結果
画像を3D化した結果
ソースコード
import os import numpy as np import matplotlib.image as mpimg import matplotlib.pyplot as plt from scipy.sparse import diags, kronsum from scipy.sparse.linalg import spsolve filename = 'C:/Users/%USER%/Downloads/Normal_map_example_-_Map.png' reverse = True # テクスチャカラーGの反転 OpenGL or DirectX img = mpimg.imread(filename, format='png') img = img[:,:,:3] # アルファチャンネルを削除 img = img * 2.0 - 1.0 # 0 - 1 から -1 - 1 の範囲へ変更 ny, nx, _ = img.shape # 式(4) p = -img[:,:,0] / img[:,:,2] q = -img[:,:,1] / img[:,:,2] if reverse: q *= -1 p *= 2.0 / nx q *= 2.0 / ny # 式(7) # ノイマン境界条件の2次元ラプラシアン tx = diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(nx, nx), format='lil') ty = diags([-1.0, 2.0, -1.0], [-1, 0, 1], shape=(ny, ny), format='lil') tx[0,0] = tx[-1,-1] = 1.0 ty[0,0] = ty[-1,-1] = 1.0 A = kronsum(tx.tocsr(), ty.tocsr(), format='csr') # -(div r) divr = -(np.gradient(p, axis=1) + np.gradient(q, axis=0)) divr = divr.reshape(-1) # 解を一意に決定するため、境界上の1点にディリクレ条件を付与 A = A.tolil() A[0,:] = 0.0 A[0,0] = 1.0 A = A.tocsr() divr[0] = 0.0 # Height Map 全体のオフセット量 # ポアソン方程式を解く z = spsolve(A, divr).reshape(ny, nx) img[:,:,0] = z img[:,:,1] = img[:,:,0] img[:,:,2] = img[:,:,0] mpimg.imsave(os.path.splitext(filename)[0] + '_output.png', np.clip(img, 0, 1)) plt.imshow(img, interpolation='nearest') plt.show()
最小化問題
最後に、最小化問題を解く手法を紹介します。
以下の目的関数の W を最小化する z を求めます。
\(
\displaystyle
W = \iint_{\Omega} (|\frac{\partial}{\partial x} z - p|^{2} + |\frac{\partial}{\partial y} z - q|^{2})dxdy. \hspace{60pt} (8)
\)
上記の方程式を解くためにフーリエ変換を使用します。
2次元 フーリエ変換 の定義は以下です。ここで j は虚数単位です。
\(
\displaystyle
\hat{z}(u,v) = \iint_{\Omega} z(x,y)e^{-j(ux + vy)} dxdy \hspace{60pt} (9)
\)
2次元 逆フーリエ変換 の定義は以下です。
\(
\displaystyle
z(x,y) = \frac{1}{2 \pi} \iint_{\Omega} \hat{z}(u,v)e^{j(ux + vy)} dudv \hspace{60pt} (10)
\)
フーリエ変換の導関数の性質は以下です。
\(
\displaystyle
\begin{eqnarray}
\mathscr{F}[\frac{\partial}{\partial x} z(x,y)] &=& ju \mathscr{F}[z(x,y)], \\
\mathscr{F}[\frac{\partial}{\partial y} z(x,y)] &=& jv \mathscr{F}[z(x,y)]. \hspace{60pt} (11)
\end{eqnarray}
\)
パーセバルの定理より以下の式が成り立ちます。
\(
\displaystyle
\iint_{\Omega} |z(x,y)|^{2} dxdy = \frac{1}{2 \pi} \iint_{\Omega} |\hat{z}(u,v)|^{2} dudv. \hspace{60pt} (12)
\)
式 (8) (11) (12) より 次の式が得られます。
\(
\displaystyle
\frac{1}{2 \pi} \iint_{\Omega} (|ju \hat{z}(u,v) - \hat{p}(u,v)|^{2} + |jv \hat{z}(u,v) - \hat{q}(u,v)|^{2}) dudv → minimum. \hspace{60pt} (13)
\)
上記の式を展開すると次の式が得られます。ここで * は共役の記号です。
\(
\displaystyle
\begin{eqnarray}
&& \frac{1}{2\pi} \iint_{\Omega} (u^{2} \hat{z} \hat{z}^{*} - ju \hat{z} \hat{p}^{*} + ju \hat{z}^{*} \hat{p} + \hat{p} \hat{p}^{*} \\
&& \hspace{30pt} + v^{2} \hat{z} \hat{z}^{*} - jv \hat{z} \hat{q}^{*} + jv \hat{z}^{*} \hat{q} + \hat{q} \hat{q}^{*}) dudv → minimum. \hspace{60pt} (14)
\end{eqnarray}
\)
上記の式を ẑ と ẑ* で微分すると、式 (8) の最小化条件が導き出されます。
\(
\displaystyle
\begin{eqnarray}
&& (u^{2} + v^{2}) \hat{z} + ju \hat{p} + jv \hat{q} = 0, \\
&& (u^{2} + v^{2}) \hat{z}^{*} - ju \hat{p}^{*} - jv \hat{q}^{*} = 0. \hspace{60pt} (15)
\end{eqnarray}
\)
式 (15) の2つの式の和と差をそれぞれ取ると次の式が得られます。
\(
\displaystyle
\begin{eqnarray}
&& (u^{2} + v^{2})(\hat{z} + \hat{z}^{*}) + ju(\hat{p} - \hat{p}^{*}) + jv(\hat{q} - \hat{q}^{*}) = 0, \\
&& (u^{2} + v^{2})(\hat{z} - \hat{z}^{*}) + ju(\hat{p} + \hat{p}^{*}) + jv(\hat{q} + \hat{q}^{*}) = 0. \hspace{60pt} (16)
\end{eqnarray}
\)
上記の式を u2 + v2 ≠ 0 について解くと
\(
\displaystyle
\hat{z}(u,v) = \frac{-ju \hat{p}(u,v) - jv \hat{q}(u,v)}{u^{2} + v^{2}}. \hspace{60pt} (17)
\)
が得られます。
u2 + v2 = 0 は周波数 0 の時の値であるため、画像全体のオフセット値に相当します。
つまり、 p(x,y) と q(x,y) をフーリエ変換し、式 (17) を解き、逆フーリエ変換を行うと Height Map を作成することができます。
画像に周期性がない場合はゼロパディングを行うことでより正確な Height Map を作成することができます。
最小化問題の計算結果
画像を3D化した結果
ソースコード
import os import numpy as np import matplotlib.image as mpimg import matplotlib.pyplot as plt filename = 'C:/Users/%USER%/Downloads/Normal_map_example_-_Map.png' reverse = True # テクスチャカラーGの反転 OpenGL or DirectX zeropadding = 3 # ゼロパディングを行う画像サイズの倍数 1のときゼロパディングなし zeropadding = max(1, zeropadding) img = mpimg.imread(filename, format='png') img = img[:,:,:3] # アルファチャンネルを削除 img = img * 2.0 - 1.0 # 0 - 1 から -1 - 1 の範囲へ変更 ny, nx, _ = img.shape # 式(4) p = -img[:,:,0] / img[:,:,2] q = -img[:,:,1] / img[:,:,2] if reverse: q *= -1 p *= 2.0 / nx q *= 2.0 / ny # ゼロパディング if zeropadding >= 2: pady = ny * (zeropadding - 1) padx = nx * (zeropadding - 1) p = np.pad(p, ((0, pady), (0, padx)), mode='constant') q = np.pad(q, ((0, pady), (0, padx)), mode='constant') # u,vを定義 u = (2.0 * np.pi * np.fft.fftfreq(nx * zeropadding)).reshape(1, nx * zeropadding) v = (2.0 * np.pi * np.fft.fftfreq(ny * zeropadding)).reshape(ny * zeropadding, 1) # 順方向2次元高速フーリエ変換 P = np.fft.fft2(p) Q = np.fft.fft2(q) # 式(17) denom = u*u + v*v mask = denom <= 1e-12 denom[mask] = 1.0 # ゼロ除算回避 Z = ( u * P.imag + v * Q.imag) / denom \ + (-u * P.real - v * Q.real) / denom * 1.0j Z[mask] = 0.0 + 0.0j # 逆方向2次元高速フーリエ変換 z = np.fft.ifft2(Z) img[:,:,0] = z.real[:ny,:nx] img[:,:,1] = img[:,:,0] img[:,:,2] = img[:,:,0] mpimg.imsave(os.path.splitext(filename)[0] + '_output.png', np.clip(img, 0, 1)) plt.imshow(img, interpolation='nearest') plt.show()
まとめ
今回は Height Map と Normal Map の相互変換についていくつかの手法を解説してみました。この記事が Height Map と Normal Map のより深い理解へとつながれば幸いです。
参考
Reinhard Klette, Karsten Schluns, "Height data from gradient fields", Jan 1996