koturnの日記

普通の人です.ブログ上のコードはコピペ自由です.

レイマーチングにおける正規化法線ベクトルの計算と実装

たまには趣向を変えて,ですます調で記事を書いてみようと思います.

レイマーチングではライティングを行うために正規化法線ベクトルの計算が必要となります. この記事ではいくつかの正規化法線ベクトルの計算手法について述べることにします.

なお,僕は普段Unityに触れている機会が多いため,本記事ではshaderlab(Cg/HLSL)をシェーダー言語として取り上げることにします. GLSLについては,記事中ではわずかに言及するのみですが,ShadertoyにGLSLでの実装コードを置いています

陰関数の法線

$\boldsymbol{p} = \begin{pmatrix} x & y & z \end{pmatrix}^T$ とします. これは3次元空間における座標の表現となるわけですね.

よく知られた事実として,陰関数 $f(\boldsymbol{p}) = 0$ の正規化法線ベクトル $\boldsymbol{n}(\boldsymbol{p})$ は以下のようになります.

\begin{eqnarray} \boldsymbol{n}(\boldsymbol{p}) & = & normalize(\nabla f(\boldsymbol{p})) \nonumber \\ & = & normalize \left( \begin{pmatrix} \dfrac{\partial f}{\partial x}(\boldsymbol{p}) & \dfrac{\partial f}{\partial y}(\boldsymbol{p}) & \dfrac{\partial f}{\partial z}(\boldsymbol{p}) \end{pmatrix}^T \right) \label{NormalizedNormal} \end{eqnarray}

レイマーチングの符号付き距離関数はまさしく陰関数であるため,式\eqref{NormalizedNormal}を基に正規化法線ベクトルを求めることができるわけです.

正規化法線ベクトルの計算と実装

本章でいくつかの正規化法線ベクトルの計算と実装を紹介します.

目次兼リンクを作ってみたのですが,大きくわけて4個,細かくわけて14個というなかなかの数となりました.

中心差分による法線計算

最初に中心差分による計算手法を紹介します. 中心差分による正規化法線ベクトルの計算は以下に基づいています. $h$ は微小な正の実数です.

\begin{eqnarray} \boldsymbol{n}(\boldsymbol{p}) & = normalize \left( \begin{pmatrix} \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right)}{2h} \\ \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right)}{2h} \\ \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right)}{2h} \end{pmatrix} \right) \nonumber \\ & = normalize \left( \begin{pmatrix} f \left(\boldsymbol{p} + \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) - f \left(\boldsymbol{p} - \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) \end{pmatrix} \right) \label{CentralDifference} \end{eqnarray}

正規化を行うので定数倍の項である $\dfrac{1}{2h}$ は除外できます. これは明らかな事実で,任意のベクトル $\boldsymbol{v}$ と $a$ ($a \neq 0$)の乗算を考えると,正規化について下記のようになるからですね.

\begin{equation} normalize(a \boldsymbol{v}) = \dfrac{a \boldsymbol{v}}{\| a \boldsymbol{v} \|} = \dfrac{a \boldsymbol{v}}{a \| \boldsymbol{v} \|} = \dfrac{\boldsymbol{v}}{\| \boldsymbol{v} \|} = normalize(\boldsymbol{v}) \end{equation}

中心差分

式\eqref{CentralDifference}をそのままシェーダーのレイマーチングの正規化法線ベクトルの計算関数に落とし込むと以下のようになります. 数式をそのままコードにしましたといえるものになっていますね.

/*!
 * @brief Calculate normal of the objects using central differences.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal01(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float2 d = float2(h, 0.0);

    return normalize(
        float3(
            map(p + d.xyy) - map(p - d.xyy),
            map(p + d.yxy) - map(p - d.yxy),
            map(p + d.yyx) - map(p - d.yyx)));
}

上記から生成されるDirect3D11用のアセンブリコードのうち,上記の部分のみを抜粋すると下記の通りです. map() 関数は後述する半径0.5の球を描画するだけのコードとしています. 同じような処理が6回並んでいるのが確認できると思います. ループで書いていないので勝手にループ処理にされるわけはないですが,ここが map() の呼び出し6回分に相当しています.

  30: add r1.xyz, r1.xyzx, cb3[3].xyzx
  31: add r2.xyz, r0.xyzx, l(0.000100, 0.000000, 0.000000, 0.000000)
  32: dp3 r0.w, r2.xyzx, r2.xyzx
  33: sqrt r0.w, r0.w
  34: add r0.w, r0.w, l(-0.500000)
  35: add r2.xyz, r0.xyzx, l(-0.000100, -0.000000, -0.000000, 0.000000)
  36: dp3 r1.w, r2.xyzx, r2.xyzx
  37: sqrt r1.w, r1.w
  38: add r1.w, r1.w, l(-0.500000)
  39: add r2.x, r0.w, -r1.w
  40: add r3.xyz, r0.xyzx, l(0.000000, 0.000100, 0.000000, 0.000000)
  41: dp3 r0.w, r3.xyzx, r3.xyzx
  42: sqrt r0.w, r0.w
  43: add r0.w, r0.w, l(-0.500000)
  44: add r3.xyz, r0.xyzx, l(-0.000000, -0.000100, -0.000000, 0.000000)
  45: dp3 r1.w, r3.xyzx, r3.xyzx
  46: sqrt r1.w, r1.w
  47: add r1.w, r1.w, l(-0.500000)
  48: add r2.y, r0.w, -r1.w
  49: add r3.xyz, r0.xyzx, l(0.000000, 0.000000, 0.000100, 0.000000)
  50: dp3 r0.w, r3.xyzx, r3.xyzx
  51: sqrt r0.w, r0.w
  52: add r0.xyzw, r0.xyzw, l(-0.000000, -0.000000, -0.000100, -0.500000)
  53: dp3 r0.x, r0.xyzx, r0.xyzx
  54: sqrt r0.x, r0.x
  55: add r0.x, r0.x, l(-0.500000)
  56: add r2.z, -r0.x, r0.w
  57: dp3 r0.x, r2.xyzx, r2.xyzx
  58: rsq r0.x, r0.x
  59: mul r0.xyz, r0.xxxx, r2.xyzx

このコードにおける map() は符号付き距離関数を統合する関数です. 符号付き距離関数そのものと考えてもよいです.

例えば単純な球を描画するなら以下のようになります. 先程の出力アセンブリコード例では map() を下記の実装にしていました.

/*!
 * @brief SDF of Sphere.
 * @param [in] p  Position of the tip of the ray.
 * @param [in] r  Radius of the sphere.
 * @return Signed Distance to the Sphere.
 */
float3 sdSphere(float3 p, float r)
{
    return length(p) - r;
}

/*!
 * @brief SDF integration function.
 * @param [in] p  Position of the tip of the ray.
 * @return Signed Distance to the one of object.
 */
float3 map(float3 p)
{
    reutrn sdSphere(p, 0.5);
}

発展させて,2つの球を合体させるなら下記のようになります.

/*!
 * @brief SDF of Sphere.
 * @param [in] r  Radius of the sphere.
 * @return Signed Distance to the Sphere.
 */
float3 sdSphere(float3 p, float r)
{
    return length(p) - r;
}

/*!
 * @brief SDF integration function.
 * @param [in] p  Position of the tip of the ray.
 * @return Signed Distance to the one of object.
 */
float3 map(float3 p)
{
    reutrn min(
        sdSphere(p - float3(0.25, 0.0, 0.0), 0.5),
        sdSphere(p + float3(0.25, 0.0, 0.0), 0.5));
}

中心差分(ループ)

以下のように6つのベクトルを用意します.

\begin{equation} \begin{cases} \boldsymbol{k}_0 & = & \begin{pmatrix}1 & 0 & 0\end{pmatrix}^T \\ \boldsymbol{k}_1 & = & \begin{pmatrix}-1 & 0 & 0\end{pmatrix}^T \\ \boldsymbol{k}_2 & = & \begin{pmatrix}0 & 1 & 0\end{pmatrix}^T \\ \boldsymbol{k}_3 & = & \begin{pmatrix}0 & -1 & 0\end{pmatrix}^T \\ \boldsymbol{k}_4 & = & \begin{pmatrix}0 & 0 & 1\end{pmatrix}^T \\ \boldsymbol{k}_5 & = & \begin{pmatrix}0 & 0 & -1\end{pmatrix}^T \end{cases} \label{SelectorVector01} \end{equation}

これを用いると,式\eqref{CentralDifference}は以下のように総和の形で記述し直すことができます.

\begin{equation} \boldsymbol{n}(\boldsymbol{p}) = normalize \left( \sum_{i} \boldsymbol{k}_i f(\boldsymbol{p} + h \boldsymbol{k}_i) \right) \label{CentralDifferenceSum} \end{equation}

$\boldsymbol{k}_i$ は各次元ごとに独立しており,正負の符号を与えるものになっています. そのため,このベクトル,定数倍のベクトルの和はどれか1つの次元のみの加算になるわけなので,式\eqref{CentralDifference}を総和の形に変形できたわけですね.

さて,式\eqref{CentralDifferenceSum}はシェーダープログラムとしてはループを用いて記述でき,以下のようになります.

/*!
 * @brief Calculate normal of the objects using central differences with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal01Loop(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);  // used only for generating k.
    static const float3 k[6] = {s.xzz, s.yzz, s.zxz, s.zyz, s.zzx, s.zzy};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 6; i++) {
        normal += k[i] * map(p + h * k[i]);
    }

    return normalize(normal);
}

変数 s というのはsignの意味で s という名前にしています. 都合上, 0も必要なので,第3要素を 0.0 にしています. どの位置にどの値を入れるかは迷いどころでしたが,0を第3要素としているとswizzle演算子を用いたときに, s.z と書くことができ,zからZeroを連想することができるからこの順にしようと考えた次第です.

UNITY_LOOP を指定していますが,これは環境差吸収用のマクロであり,例えば [loop] に展開されます. この指示があるとコンパイラは積極的にループ命令を生成してくれるようになります.

  31: mov r2.xyz, l(0,0,0,0)
  32: mov r0.w, l(0)
  33: loop 
  34:   ige r1.w, r0.w, l(6)
  35:   breakc_nz r1.w
  36:   mad r3.xyz, icb[r0.w + 0].xyzx, l(0.000100, 0.000100, 0.000100, 0.000000), r0.xyzx
  37:   dp3 r1.w, r3.xyzx, r3.xyzx
  38:   sqrt r1.w, r1.w
  39:   add r1.w, r1.w, l(-0.500000)
  40:   mad r2.xyz, icb[r0.w + 0].xyzx, r1.wwww, r2.xyzx
  41:   iadd r0.w, r0.w, l(1)
  42: endloop 
  43: dp3 r0.x, r2.xyzx, r2.xyzx
  44: rsq r0.x, r0.x
  45: mul r0.xyz, r0.xxxx, r2.xyzx

逆にUNITY_UNROLL[unroll]) を指定するとコンパイラは積極的にループ展開を行ってくれるようになります. この場合に生成されるコードは 中心差分 のものとほぼ同一になります.

  30: add r1.xyz, r1.xyzx, cb3[3].xyzx
  31: add r2.xyz, r0.xyzx, l(0.000100, 0.000000, 0.000000, 0.000000)
  32: dp3 r0.w, r2.xyzx, r2.xyzx
  33: sqrt r0.w, r0.w
  34: add r0.w, r0.w, l(-0.500000)
  35: add r2.xyz, r0.xyzx, l(-0.000100, 0.000000, 0.000000, 0.000000)
  36: dp3 r1.w, r2.xyzx, r2.xyzx
  37: sqrt r1.w, r1.w
  38: add r1.w, r1.w, l(-0.500000)
  39: mul r2.xyz, r1.wwww, l(-1.000000, 0.000000, 0.000000, 0.000000)
  40: mad r2.xyz, r0.wwww, l(1.000000, 0.000000, 0.000000, 0.000000), r2.xyzx
  41: add r3.xyz, r0.xyzx, l(0.000000, 0.000100, 0.000000, 0.000000)
  42: dp3 r0.w, r3.xyzx, r3.xyzx
  43: sqrt r0.w, r0.w
  44: add r0.w, r0.w, l(-0.500000)
  45: mad r2.xyz, r0.wwww, l(0.000000, 1.000000, 0.000000, 0.000000), r2.xyzx
  46: add r3.xyz, r0.xyzx, l(0.000000, -0.000100, 0.000000, 0.000000)
  47: dp3 r0.w, r3.xyzx, r3.xyzx
  48: sqrt r0.w, r0.w
  49: add r0.w, r0.w, l(-0.500000)
  50: mad r2.xyz, r0.wwww, l(0.000000, -1.000000, 0.000000, 0.000000), r2.xyzx
  51: add r3.xyz, r0.xyzx, l(0.000000, 0.000000, 0.000100, 0.000000)
  52: dp3 r0.w, r3.xyzx, r3.xyzx
  53: sqrt r0.w, r0.w
  54: add r0.xyzw, r0.xyzw, l(0.000000, 0.000000, -0.000100, -0.500000)
  55: mad r2.xyz, r0.wwww, l(0.000000, 0.000000, 1.000000, 0.000000), r2.xyzx
  56: dp3 r0.x, r0.xyzx, r0.xyzx
  57: sqrt r0.x, r0.x
  58: add r0.x, r0.x, l(-0.500000)
  59: mad r0.xyz, r0.xxxx, l(0.000000, 0.000000, -1.000000, 0.000000), r2.xyzx
  60: dp3 r0.w, r0.xyzx, r0.xyzx
  61: rsq r0.w, r0.w
  62: mul r0.xyz, r0.wwww, r0.xyzx

何も指示しなかった場合, map() のコードサイズに応じてループ展開を行うかどうかをコンパイラが決定してくれます. 単純なスフィアトレーシングであればループ展開され,処理の多い複雑な map() だとループ命令が生成されることは確認しました.

わざわざループにするのはプログラムの高速化の考えに反しますが,一般的に map() は長大な関数になりがちであるため,ループにすることで出力コードサイズの削減になります. また,どの部分が法線計算であるか明確になるため,出力コードが読みやすくなるというメリットもあります. そもそもマーチングループの中で map() を呼び出しているため,法線計算のための数回の呼び出しをループにすることはそこまでの抵抗を感じるものではないかなと思います.

ところで,下記のように hk の積のルックアップテーブルを作成する実装も考えることができます.

/*!
 * @brief Calculate normal of the objects using central differences with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal01LoopEx(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);  // used only for generating hs, k and hk.
    static const float3 hs = h * s;  // used only for generating hk.
    static const float3 k[6] = {s.xzz, s.yzz, s.zxz, s.zyz, s.zzx, s.zzy};
    static const float3 hk[6] = {hs.xzz, hs.yzz, hs.zxz, hs.zyz, hs.zzx, hs.zzy};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 6; i++) {
        normal += k[i] * map(p + hk[i]);
    }

    return normalize(normal);
}

しかし,p + h * k[i] のような積和の計算は単一のmad命令となるため,p + hk[i] のような単一の加算命令と比較しても速度的な差はないでしょう.

中心差分(LUT無しループ)

ルックアップテーブルを用いた形のコードを提示しましたが,ルックアップテーブル用のConstant Bufferが作成されることになります. その領域を削りたい場合はループインデックス $i$ から係数ベクトル $\boldsymbol{k}_i$ を作るとよいでしょう.

/*!
 * @brief Calculate normal of the objects using central differences with loop without Look Up Table.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal01LoopNoLut(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 6; i++) {
        const int j = i >> 1;
        const float4 v = float4(int4((int3(j + 3, i, j) >> 1), i) & 1);
        const float3 k = v.xyz * (v.w * 2.0 - 1.0);
        normal += k * map(p + h * k);
    }

    return normalize(normal);
}

float3 6つ分のConstant Bufferを削りたいモチベーションはそこまでなく,このコードは単に興味の一環で作成したというものになります. もちろん k の計算を行う分,ループ内の命令数がいくつか増えることになります. DirectX11向けのアセンブリコードの当該箇所を抜粋すると下記の通りで,8命令になっていることがわかります(あらかじめ i >> 1 の一時変数を用意したり,整数ベクトルにまとめてビットシフトとビット論理積を行ったり,mad命令が生成されるように浮動小数点ベクトルに変換してから * 2.0 - 1.0 を行って命令数が少なくなるようにしています).

  36:   ishr r3.y, r0.w, l(1)
  37:   iadd r3.x, r3.y, l(3)
  38:   mov r3.w, r0.w
  39:   ishr r3.xyz, r3.xwyx, l(1)
  40:   and r3.xyzw, r3.xyzw, l(1, 1, 1, 1)
  41:   itof r3.xyzw, r3.xyzw
  42:   mad r1.w, r3.w, l(2.000000), l(-1.000000)
  43:   mul r3.xyz, r1.wwww, r3.xyzx

やや複雑なビット演算を行っているので,テクニカルなコード書ける俺スゲー感を得たい場合にはよいかもしれませんね.

中心差分(2重ループ)

最初に言っておくと,これは採用する理由のない計算手法になります.

式\eqref{SelectorVector01}を眺めてみると,

\begin{equation} \begin{cases} \boldsymbol{k}_0 & = & -\boldsymbol{k}_1 \\ \boldsymbol{k}_2 & = & -\boldsymbol{k}_3 \\ \boldsymbol{k}_4 & = & -\boldsymbol{k}_5 \end{cases} \end{equation}

となっていることがわかります. そこで,

\begin{equation} \begin{cases} \boldsymbol{k}'_0 & = & \begin{pmatrix}1 & 0 & 0\end{pmatrix}^T \\ \boldsymbol{k}'_1 & = & \begin{pmatrix}0 & 1 & 0\end{pmatrix}^T \\ \boldsymbol{k}'_2 & = & \begin{pmatrix}0 & 0 & 1\end{pmatrix}^T \end{cases} \label{SelectorVector02} \end{equation}

\begin{equation} \begin{cases} s_0 & = & 1 \\ s_1 & = & -1 \end{cases} \end{equation}

とおけば,式\eqref{CentralDifferenceSum}は以下のように2重の総和の形になります.

\begin{equation} \boldsymbol{n}(\boldsymbol{p}) = normalize \left( \sum_{i} \sum_{j} s_j \boldsymbol{k}'_i f(\boldsymbol{p} + h s_j \boldsymbol{k}'_i) \right) \end{equation}

プログラムとしては,以下のような2重ループの形になります.

/*!
 * @brief Calculate normal of the objects using central differences with double loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal01LoopW(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);
    static const float3 k[3] = {s.xzz, s.zxz, s.zzx};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 3; i++) {
        UNITY_LOOP
        for (int j = 0; j < 2; j++) {
            const float3 sk = s[j] * k[i];
            normal += sk * map(p + h * sk);
        }
    }

    return normalize(normal);
}

しかし,これは $s_j$ と $\boldsymbol{k}'_i$ の積を計算する必要があったり,二重ループのためループ指示を2箇所に書く必要があるなど,採用する理由がない計算手法となります.

前方差分による法線計算

中心差分による関数 $f$ すなわち,map() の評価回数は6回でした. レイマーチングにおいて map() は大部分を占める関数であり,基本的に高コストな関数であるため,極力呼び出し回数は減らしたいものです.

そこで,計算誤差の増加と引き換えに,前方差分により法線計算を行うことを考えてみます.

\begin{eqnarray} \boldsymbol{n}(\boldsymbol{p}) & = normalize \left( \begin{pmatrix} \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right)}{2h} \\ \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right)}{2h} \\ \dfrac{f \left(\boldsymbol{p} + \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right)}{2h} \end{pmatrix} \right) \nonumber \\ & = normalize \left( \begin{pmatrix} f \left(\boldsymbol{p} + \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) - f \left(\boldsymbol{p} \right) \end{pmatrix} \right) \label{ForwardDifference} \end{eqnarray}

$f(p)$ の値を保持すれば,関数 $f$ の評価回数は4回となりますね.

前方差分

式\eqref{ForwardDifference}に従い,愚直に実装すると以下のようになります. 数式通りのコードですね.

/*!
 * @brief Calculate normal of the objects using forward differences.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal02(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float2 d = float2(h, 0.0);

    const float mp = map(p);

    return normalize(
        float3(
            map(p + d.xyy) - mp,
            map(p + d.yxy) - mp,
            map(p + d.yyx) - mp));
}

前方差分(ループ)

式\eqref{ForwardDifference}に対し,式\eqref{SelectorVector02}の $\boldsymbol{k}'_i$ を用いると,以下のようになります.

\begin{eqnarray} \boldsymbol{n}(\boldsymbol{p}) & = & normalize \left( \sum_{i} \boldsymbol{k}'_i \left\{ f(\boldsymbol{p} + h \boldsymbol{k}'_i) - f(\boldsymbol{p}) \right\} \right) \nonumber \\ & = & normalize \left( \sum_{i} \boldsymbol{k}'_i f(\boldsymbol{p} + h \boldsymbol{k}'_i) - \sum_{i} \boldsymbol{k}'_i f(\boldsymbol{p}) \right) \nonumber \\ & = & normalize \left( \sum_{i} \boldsymbol{k}'_i f(\boldsymbol{p} + h \boldsymbol{k}'_i) - \begin{pmatrix} f(\boldsymbol{p}) & f(\boldsymbol{p}) & f(\boldsymbol{p}) \end{pmatrix}^T \right) \end{eqnarray}

/*!
 * @brief Calculate normal of the objects using forward differences with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal02Loop(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);  // used only for generating k.
    static const float3 k[3] = {s.xzz, s.zxz, s.zzx};

    float3 normal = (-map(p)).xxx;

    UNITY_LOOP
    for (int i = 0; i < 3; i++) {
        normal += k[i] * map(p + h * k[i]);
    }

    return normalize(normal);
}

shaderlab(Cg/HLSL)ではスカラーに対してswizzle演算子を用いることができるため,

const float nmp = -map(p);
float3 normal = float3(nmp, nmp, nmp);

と書かなくても

float3 normal = (-map(p)).xxx;

のように1行で記述できるのが楽でよいですね.

GLSLではスカラーに対してswizzle演算子を用いることができませんが,ベクトルの生成に1要素のみを指定することで,全要素が同じ値のベクトルを作ることができるので,これも便利な文法ですね.

vec3 normal = vec3(-map(p));

前方差分(LUT無しループ)

これもループインデックス $i$ から係数ベクトル $\boldsymbol{k}'_i$ を作ることができます.

/*!
 * @brief Calculate normal of the objects using forward differences with loop without Look Up Table.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal02LoopNoLut(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value

    float3 normal = (-map(p)).xxx;

    UNITY_LOOP
    for (int i = 0; i < 3; i++) {
        const float3 k = float3(int3((i + 3) >> 1, i, i >> 1) & 1);
        normal += k * map(p + h * k);
    }

    return normalize(normal);
}

前方差分(NGループ)

この計算手法も採用する理由がないものとなります.

前方差分(ループ)ではループの形にしたものの, map(p) の計算をループ外に置かざるをえず,コンパイラがループを生成したとしても,map(p) の呼び出し箇所がループの外と中の2箇所となり,コードサイズ面で中心差分(ループ)より不利になります.

ここで思考実験として, map(p) の計算を無理矢理ループ内に納めることを考えてみましょう. 以下の6つの係数を設けると,

\begin{equation} \begin{cases} z_0 & = & 1 \\ z_1 & = & 0 \\ z_2 & = & 1 \\ z_3 & = & 0 \\ z_4 & = & 1 \\ z_5 & = & 0 \end{cases} \end{equation}

式\eqref{SelectorVector01}の $\boldsymbol{k}_i$ と併用し,以下のように正規化法線ベクトルを計算することができます.

\begin{equation} \boldsymbol{n}(\boldsymbol{p}) = normalize \left( \sum_{i} \boldsymbol{k}_i f(\boldsymbol{p} + h z_i \boldsymbol{k}_i) \right) \end{equation}

これをシェーダープログラムにすると以下のようになります. しかし,これでは map() の呼び出し回数が6回となり,わざわざ精度を犠牲にしてまで前方差分を採用した意味がなくなってしまいますね.

/*!
 * @brief Calculate normal of the objects using forward differences with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal02LoopNG(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);  // used only for generating k.
    static const float3 k[6] = {s.xzz, s.yzz, s.zxz, s.zyz, s.zzx, s.zzy};
    static const float z[6] = {1.0, 0.0, 1.0, 0.0, 1.0, 0.0};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 6; i++) {
        normal += k[i] * map(p + h * z[i] * k[i]);
    }

    return normalize(normal);
}

前方差分による法線計算(f(p)無視)

レイマーチングにおいて, $f(\boldsymbol{p}) \approx 0$ となるように $\boldsymbol{p}$ を求めていました. そこで大胆にも式\eqref{ForwardDifference}において,$f(\boldsymbol{p}) = 0$ とおき,法線を計算するという考えがあるようです.

\begin{equation} \boldsymbol{n}(\boldsymbol{p}) = normalize \left( \begin{pmatrix} f \left(\boldsymbol{p} + \begin{pmatrix} h & 0 & 0 \end{pmatrix}^T \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & h & 0 \end{pmatrix}^T \right) \\ f \left(\boldsymbol{p} + \begin{pmatrix} 0 & 0 & h \end{pmatrix}^T \right) \end{pmatrix} \right) \end{equation}

実際に実装してみるとわかるのですが,かなり精度が悪く,描画結果が残念なものになりました. レイマーチングの打ち切り条件をかなり厳しくしないと良い描画結果が得られないのではないでしょうか?

前方差分(f(p)無視)

map(p) の減算が無いため,コードはかなりシンプルになります.

/*!
 * @brief Calculate normal of the objects ignoring f(p).
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal03(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float2 d = float2(h, 0.0);

    return normalize(
        float3(
            map(p + d.xyy),
            map(p + d.yxy),
            map(p + d.yyx)));
}

前方差分(f(p)無視)(ループ)

前方差分(ループ)と同様です.

/*!
 * @brief Calculate normal of the objects ignoring f(p)
 * and using forward differences with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal03Loop(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float3 s = float3(1.0, -1.0, 0.0);  // used only for generating k.
    static const float3 k[3] = {s.xzz, s.zxz, s.zzx};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 3; i++) {
        normal += k[i] * map(p + h * k[i]);
    }

    return normalize(normal);
}

前方差分(f(p)無視)(LUT無しループ)

前方差分(LUT無しループ)と同様です.

/*!
 * @brief Calculate normal of the objects ignoring f(p)
 * and using forward differences with loop without Look Up Table.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal03LoopNoLut(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 3; i++) {
        const float3 k = float3(int3((i + 3) >> 1, i, i >> 1) & 1);
        normal += k * map(p + h * k);
    }

    return normalize(normal);
}

Tetrahedron Technique による法線計算

Tetrahedron TechniqueとはInigo Quilez氏の記事で紹介されている正規化法線ベクトルの計算手法の1つです. 割と有名な手法かもしれませんね.

まず,以下のように4つのベクトルを定めます.

\begin{equation} \begin{cases} \boldsymbol{k}''_0 & = & \begin{pmatrix}1 & -1 & -1\end{pmatrix}^T \\ \boldsymbol{k}''_1 & = & \begin{pmatrix}-1 & -1 & 1\end{pmatrix}^T \\ \boldsymbol{k}''_2 & = & \begin{pmatrix}-1 & 1 & -1\end{pmatrix}^T \\ \boldsymbol{k}''_3 & = & \begin{pmatrix}1 & 1 & 1\end{pmatrix}^T \end{cases} \label{SelectorVector03} \end{equation}

すると,正規化法線ベクトルは以下のように計算できます.

\begin{equation} \boldsymbol{n}(\boldsymbol{p}) = normalize \left( \sum_i \boldsymbol{k}''_i f(\boldsymbol{p} + h \boldsymbol{k}''_i) \right) \label{TetrahedronTechnique} \end{equation}

導出過程については,Inigo Quilez氏の記事を参照してください. $\sum_i \boldsymbol{k}''_i f(\boldsymbol{p}) = \boldsymbol{0}$であることと,方向微分の公式を用いてうまく式変形しているので,感銘を受けるものがあります.

式\eqref{CentralDifferenceSum}と式\eqref{TetrahedronTechnique}が同じ形であること,また,式\eqref{SelectorVector01}の$\boldsymbol{k}_i$について

\begin{equation} \sum_i \boldsymbol{k}_i = \boldsymbol{0} \end{equation}

そして,式\eqref{SelectorVector03}の$\boldsymbol{k}''_i$についても

\begin{equation} \sum_i \boldsymbol{k}''_i = \boldsymbol{0} \end{equation}

となることから,式\eqref{TetrahedronTechnique}のように正規化法線ベクトルが計算できることは自然なことであるように直感的に感じることができます.

Tetrahedron Technique

数式の通りにコードを書くと以下のようになります. ライブコーディング等でゼロから記述する際に,提案した手法の中では最も楽に記述できるコードかつ美しいのではないだろうかと個人的には思います.

/*!
 * @brief Calculate normal of the objects using tetrahedron techniue.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal04(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float2 s = float2(1.0, -1.0);
    static const float2 hs = h * s;

    return normalize(
        s.xyy * map(p + hs.xyy)
            + s.yxy * map(p + hs.yxy)
            + s.yyx * map(p + hs.yyx)
            + map(p + hs.xxx).xxx);
}

数式通りなら,map(p + hs.xxx).xxx の部分は s.xxx * map(p + hs.xxx) ですが,s.xxx は全要素 1.0 のベクトルなので,この2つは同一の結果になるわけです. まぁ,ベクトルに対するスカラーの加算は,スカラー側が全要素同じ値のベクトルとして扱われるため,単に map(p + hs.xxx) と書いてもよかったのですが,数学的にはベクトルとスカラーの乗算はあってもベクトルとスカラーの加算はなく,直感的ではないということで前述の表現を用いることにしました.

Tetrahedron Technique(ループ)

式\eqref{TetrahedronTechnique}が総和の形なので,それをそのまま for を用いたプログラムにしてみます.

/*!
 * @brief Calculate normal of the objects using tetrahedron techniue with loop.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal04Loop(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value
    static const float2 s = float2(1.0, -1.0);  // used only for generating k.
    static const float3 k[4] = {s.xyy, s.yxy, s.yyx, s.xxx};

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 4; i++) {
        normal += k[i] * map(p + h * k[i]);
    }

    return normalize(normal);
}

Tetrahedron Technique(LUT無しループ)

これも $i$ から $\boldsymbol{k}''_i$ を生成することができ,ルックアップテーブルが不要となります.

/*!
 * @brief Calculate normal of the objects using tetrahedron techniue with loop without Look Up Table.
 * @param [in] p  Position of the tip of the ray.
 * @return Normal of the objects.
 */
float3 getNormal04LoopNoLut(float3 p)
{
    static const float h = 0.0001;  // replace by an appropriate value

    float3 normal = float3(0.0, 0.0, 0.0);

    UNITY_LOOP
    for (int i = 0; i < 4; i++) {
        const float3 k = float3(int3((i + 3) >> 1, i, i >> 1) & 1) * 2.0 - 1.0;
        normal += k * map(p + h * k);
    }

    return normalize(normal);
}

コード中の int3 は要素が 01 のベクトルです. これに2を掛け,1を減ずることで, -11 のベクトルになります.

各手法の長所と短所

各手法の長所と短所をまとめると以下の通りとなります.

手法 SDF評価回数 ループ利用時SDF呼び出し箇所数 誤差
中心差分 6 1
前方差分 4 2
前方差分(f(p)無視) 3 1
Tetrahedron Technique 4 1 おそらく前方差分より小

総合的に見て,一番バランスが良いのがTetrahedron Techniqueだと思います.

まとめ

本記事では下記3つの手法それぞれの正規化法線ベクトルの計算手法を示しました.

  • 中心差分
  • 前方差分
  • Tetrahedron Technique

またそれぞれの手法についてループを用いた実装も示しました. ループを用いた実装にする利点としては以下の通りです.

  • 生成コードサイズが小さくなる
    • 生成コードの可読性が向上する
  • UNITY_UNROLL[unroll]) を指定することで容易にループ展開可能

おまけ

本記事ではshaderlab(Cg/HLSL)での実装例を示しましたが,ShadertoyにGLSLでの実装例も置いてあります. GLSLの実装コードが欲しい場合は見てください.

ただし,ルックアップテーブルを作成する際の配列の宣言と同時に初期化を行っている部分と,ルックアップテーブルを使用しない場合の係数ベクトル作成における $i$ のビット演算はGLES 3.0以降でなければ利用できないため注意が必要です.

配列については宣言と同時に初期化しなければ利用可能なので,例えば

const vec3 k[4] = vec3[](s.xyy, s.yxy, s.yyx, s.xxx);

vec3 k[4];
k[0] = s.xyy;
k[1] = s.yxy;
k[2] = s.yyx;
k[3] = s.xxx;

とすれば,古いGLSLでもコンパイルが通るコードとなります.

また,Inigo Quilez氏の記事の最後では,GLSLでループ命令を生成するためのHackが紹介されています. GLSLでは [loop] のようなコンパイラへの指示がないため,あえてループ命令を生成しようとした場合,uniform変数に依存した定数ではないゼロを作る等の工夫が必要になるようです.

記事中のコードを本記事に合わせた関数名・変数名に変更すると以下のようになります. 謎に $0.5773$ を掛けている部分がありますが, $\dfrac{1}{\sqrt{3}} \approx 0.5773$ なので, k の正規化を図っているのでしょう.

vec3 getNormal(vec3 p)
{
    const float h = 0.0001;  // replace by an appropriate value
#define ZERO (min(iFrame, 0))  // non-constant zero

    vec3 normal = vec3(0.0);
    for (int i = ZERO; i < 4; i++) {
        vec3 k = 0.5773 * (2.0 * vec3(((i + 3) >> 1) & 1, (i >> 1) & 1, i & 1) - 1.0);
        normal += k * map(pos + h * k);
    }

    return normalize(normal);
}

このようなHackが必要になるのは,特定のプラットフォーム(特にWebGL)におけるシェーダーコンパイラmap() が長大である場合,利用可能な命令数を超過し,クラッシュすることがあるためのようです.

参考文献