陰関数 $f$ （ $f: \mathbb{R}^3 \mapsto \mathbb{R}$ ）の正規化法線ベクトル $\boldsymbol{n}$ （ $\boldsymbol{n} \in \mathbb{R^3}$ ） は下記のように求められる．
$$
\begin{align}
\boldsymbol{n} &= normalize(\nabla f( \boldsymbol{p} )) \\
&= normalize \begin{pmatrix}
\dfrac{\partial f(\boldsymbol{p})}{\partial x}
\\
\dfrac{\partial f(\boldsymbol{p})}{\partial y}
\\
\dfrac{\partial f(\boldsymbol{p})}{\partial z}
\end{pmatrix} \\
&= normalize \begin{pmatrix}
\lim_{h \to 0} \dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}h & 0 & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\\
\lim_{h \to 0} \dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}h & 0 & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\\
\lim_{h \to 0} \dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}h & 0 & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\end{pmatrix}
\end{align}
$$
$h$ を極小の値とし，数値微分を考えると，正規化により定数倍は無視できるので下記のようになる．
$$
\begin{align}
\boldsymbol{n} &= normalize(\nabla f( \boldsymbol{p} )) \\
&= normalize \begin{pmatrix}
\dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}h & 0 & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\\
\dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}0 & h & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\\
\dfrac{f\left(\boldsymbol{p} + \begin{pmatrix}0 & 0 & h \end{pmatrix}^T \right) - f(\boldsymbol{p})}{h}
\end{pmatrix} \\
&= normalize \begin{pmatrix}
f\left(\boldsymbol{p} + \begin{pmatrix}h & 0 & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})
\\
f\left(\boldsymbol{p} + \begin{pmatrix}0 & h & 0 \end{pmatrix}^T \right) - f(\boldsymbol{p})
\\
f\left(\boldsymbol{p} + \begin{pmatrix}0 & 0 & h \end{pmatrix}^T \right) - f(\boldsymbol{p})
\end{pmatrix}
\end{align}
$$
誤差が小さくなるように中央差分法を用いる場合は下記のようになる．
$$
\boldsymbol{n} = normalize \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}
$$
これをGLSLで実装すると，下記のようになる．
map() が数式の $f$ に対応する関数である．
vec3 calcNormal(vec3 p)
{
const float h = 0.0001 ; // replace by an appropriate valueconst vec2 k = vec2 (h, 0.0 );
return normalize (
vec3 (
map(p + k.xyy) - map(p - k.xyy),
map(p + k.yxy) - map(p - k.yxy),
map(p + k.yyx) - map(p - k.yyx)));
}ループを用いる場合は下記のように書ける．
レイマーチングの map() は大きくなりがちな関数であるため，ループ形式にするということはコードサイズ，コンパイル時間の面で利がある．
ただし，GLSLは積極的にインライン展開とループアンローリングを行うため，真にループを用いるコードを生成するためには細工が必要であるが割愛する．
vec3 calcNormal(vec3 p)
{
const float h = 0.0001 ; // replace by an appropriate valueconst vec3 s = vec3 (1.0 , - 1.0 , 0.0 ); // used only for generating k.const vec3 k[6 ] = vec3 [](s.xzz, s.yzz, s.zxz, s.zyz, s.zzx, s.zzy);
vec3 normal = vec3 (0.0 , 0.0 , 0.0 );
for (int i = 0 ; i < 6 ; i++ ) {
normal += k[i] * map(p + h * k[i]);
}
return normalize (normal);
}四面体法 (Tetrahedron Teqnique) 下記の4つのベクトルを設ける．
$$
\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}
$$
この4つのベクトルを用いると下記のように打ち消しが発生し，ゼロベクトルとなる．
$$
\sum_i \boldsymbol{k}_i f(\boldsymbol{p}) =
\begin{pmatrix}
f(\boldsymbol{p})_x \\
-f(\boldsymbol{p})_y \\
-f(\boldsymbol{p})_z
\end{pmatrix} +
\begin{pmatrix}
-f(\boldsymbol{p})_x \\
f(\boldsymbol{p})_y \\
-f(\boldsymbol{p})_z
\end{pmatrix} +
\begin{pmatrix}
-f(\boldsymbol{p})_x \\
-f(\boldsymbol{p})_y \\
f(\boldsymbol{p})_z
\end{pmatrix} +
\begin{pmatrix}
f(\boldsymbol{p})_x \\
f(\boldsymbol{p})_y \\
f(\boldsymbol{p})_z
\end{pmatrix}
= \boldsymbol{0}
$$
従って， $h$ を微小な値とすると，
$$
\begin{align}
\boldsymbol{m} &= \sum_i \boldsymbol{k}_i \dfrac{f(\boldsymbol{p} + h \boldsymbol{k}_i)}{h} \\
&= \sum_i \boldsymbol{k}_i \dfrac{f(\boldsymbol{p} + h \boldsymbol{k}_i)}{h} - \boldsymbol{0} \\
&= \sum_i \boldsymbol{k}_i \dfrac{f(\boldsymbol{p} + h \boldsymbol{k}_i)}{h} - \sum_i \boldsymbol{k}_i \dfrac{f(\boldsymbol{p})}{h} \\
&= \sum_i \left( \boldsymbol{k}_i \dfrac{f(\boldsymbol{p} + h \boldsymbol{k}_i)}{h} - \boldsymbol{k}_i \dfrac{f(\boldsymbol{p})}{h} \right) \\
&= \sum_i \boldsymbol{k}_i \dfrac{f(\boldsymbol{p} + h \boldsymbol{k}_i) - f(\boldsymbol{p})}{h} \\
&= \sum_i \boldsymbol{k}_i \nabla_{\boldsymbol{k}_i} f(\boldsymbol{p}) \\
&= \sum_i \boldsymbol{k}_i (\boldsymbol{k}_i \cdot \nabla f(\boldsymbol{p}))
\end{align}
$$
となる．
ここで， $\boldsymbol{m}$ の成分 $x$ だけに注目すると，以下を得る．
$$
\begin{align}
m_x &= \sum_i k_{ix} (\boldsymbol{k}_i \cdot \nabla f(\boldsymbol{p})) \\
&= \sum_i \left( (k_{ix} \boldsymbol{k}_i) \cdot \nabla f(\boldsymbol{p}) \right) \\
&= \nabla f(\boldsymbol{p}) \cdot \sum_i k_{ix} \boldsymbol{k}_i \\
&= \nabla f(\boldsymbol{p}) \cdot \begin{pmatrix} 4 & 0 & 0 \end{pmatrix}^T \\
&= 4 \dfrac{\partial f(\boldsymbol{p})}{\partial x}
\end{align}
$$
$y$ , $z$ に関しても同様であるため，最終的に $\boldsymbol{m}$ は下記のようになる．
$$
\begin{align}
\boldsymbol{m} &= \begin{pmatrix}
m_x & m_y & m_z
\end{pmatrix}^T \\
&=
\begin{pmatrix}
4 \dfrac{\partial f(\boldsymbol{p})}{\partial x} & 4 \dfrac{\partial f(\boldsymbol{p})}{\partial y} & 4 \dfrac{\partial f(\boldsymbol{p})}{\partial z}
\end{pmatrix}^T \\
&=
4 \begin{pmatrix}
\dfrac{\partial f(\boldsymbol{p})}{\partial x} & \dfrac{\partial f(\boldsymbol{p})}{\partial y} & \dfrac{\partial f(\boldsymbol{p})}{\partial z}
\end{pmatrix}^T \\
&=
4 \nabla f(\boldsymbol{p})
\end{align}
$$
正規化を行うなら定数倍は無視できるので，正規化法線ベクトル $\boldsymbol{n}$ は
$$
\begin{align}
\boldsymbol{n}
&=
normalize \left( \nabla f(\boldsymbol{p}) \right) \\
&=
normalize \left( \dfrac{1}{4} \boldsymbol{m} \right) \\
&=
normalize(\boldsymbol{m}) \\
&=
normalize \left( \sum_i \boldsymbol{k}_i \cdot f(\boldsymbol{p} + h \boldsymbol{k}_i) \right)
\end{align}
$$
となる．
これをGLSLで実装すると下記のようになる．
vec3 calcNormal(vec3 p)
{
const float h = 0.0001 ; // replace by an appropriate valueconst vec2 s = vec2 (1.0 , - 1.0 );
const vec2 hs = h * s;
return normalize (
s.xyy * map(p + hs.xyy)
+ s.yxy * map(p + hs.yxy)
+ s.yyx * map(p + hs.yyx)
+ s.xxx * map(p + hs.xxx));
}ループを用いる場合は下記のようになる．
vec3 calcNormal(vec3 p)
{
const float h = 0.0001 ; // replace by an appropriate valueconst vec2 s = vec2 (1.0 , - 1.0 ); // used only for generating k.const vec3 k[4 ] = vec3 [](s.xyy, s.yxy, s.yyx, s.xxx);
vec3 normal = vec3 (0.0 , 0.0 , 0.0 );
for (int i = 0 ; i < 4 ; i++ ) {
normal += k[i] * map(p + h * k[i]);
}
return normalize (normal);
}$f(x + h)$ と $f(x - h)$ のTaylor展開はそれぞれ下記の通り．
$$
f(x + h) = f(x) + h f'(x) + \dfrac{h^2}{2!} f''(x) + \dfrac{h^3}{3!} f'''(x) + O(h^4)
$$
$$
f(x - h) = f(x) - h f'(x) + \dfrac{h^2}{2!} f''(x) - \dfrac{h^3}{3!} f'''(x) + O(h^4)
$$
従って，
$$
\dfrac{f(x + h) - f(x)}{h} = f'(x) + \dfrac{h}{2!} f''(x) + \dfrac{h^2}{3!} f'''(x) + O(h^3) = f'(x) + O(h)
$$
$$
\dfrac{f(x) - f(x - h)}{h} = f'(x) - \dfrac{h}{2!} f''(x) + \dfrac{h^2}{3!} f'''(x) + O(h^3) = f'(x) + O(h)
$$
$$
\dfrac{f(x + h) - f(x - h)}{2h} = \dfrac{1}{2h} \left( 2 hf'(x) + 2 \dfrac{h^3}{3!} f'''(x) + O(h^5) \right) = f'(x) + O(h^2)
$$
$h$ は十分に小さな正の値であるため， $|O(h^2)| < |O(h)|$ ．
すなわち，中心差分は前方差分，後方差分より誤差が小さい．
$$
\nabla_{\boldsymbol{v}} f(\boldsymbol{x}) = \lim_{h \to 0} \dfrac{f(\boldsymbol{x} + h\boldsymbol{v}) - f(\boldsymbol{x})}{h}
$$
$$
\nabla_{\boldsymbol{v}} f(\boldsymbol{x}) = \nabla f(\boldsymbol{x}) \cdot \boldsymbol{v}
$$
