先日のポストでは DVR [Niemeyer et al., CVPR'20] を参考にDNNでSigned Distance Function (SDF)をモデル化し，学習させてみた． forwardではスフィアトレーシングによるレンダリングを行い，backwardではオブジェクトとの交点の微分を陰関数微分により求めることで，微分可能レンダリングの枠組みでDNNの最適化を試みた．しかし，現状全く学習がうまく進んでいない．
DVR [Niemeyer et al., CVPR'20] ではSigned Distance Functionではなく，Occupancy Networkをモデル化している．なぜこれをSigned Distance Functionに置き換えたのかというと，GLSLでそのままスフィアトレーシングによりレンダリングしたかったからだ． しかしOccupancy Networkの方が，例えば物体の投影された領域が分かっている場合などは，教師信号を与えることができる点で最適化しやすいように思う．また，Signed Distance Functionをモデル化している SDFDiff [Jiang et al., CVPR'20] などは3D deconvolutionによるボクセルベースのSigned Distance Functionを用いており，ビジョンにおける画像生成が2D deconvolutionによって成功していることを考えると，convolutionの力を借りないと難しいのかもしれない．
でもOccupancy Networkなどで表現されたオブジェクトなどのように，Signed Distance Functionが直接分からない場合でもスフィアトレーシングできるケースがあることを知った．
話が変わるが，フラクタルの有名な例としてJulia集合，Mandelbrot集合などがある． 以下のような2次関数が与えられた際に，
$$f_{c}(z)=z^{2}+c$$
Julia集合$\mathcal{J_{c}}$，Mandelbrot集合$\mathcal{M}$は以下で定義される．
$$\mathcal{J_{c}}=\{z: \lim_{n \to \infty}f_{c}^{n}(z) \not\to \infty\}$$ $$\mathcal{M}=\{c: \lim_{n \to \infty}f_{c}^{n}(z) \not\to \infty, f'_{c}(z)=0\}$$
$z,c$は通常複素数であり，$\mathcal{J_{c}},\mathcal{M}$を複素平面上に表すと，フラクタルが現れる．
2D Julia set & Mandelbrot set
ここでこれら$\mathcal{J_{c}},\mathcal{M}$を3次元に拡張したいと思うのは自然なことであり，$z,c$を四元数（ クォータニオン）として表して，最後に適当な超平面に射影すれば，実際に3次元のフラクタルを得ることができる． ここでなぜ三元数を考えないのかと思ったが，$f_{c}$を適用することを考えると，$z,c$が定義される代数系は加法と乗法について閉じている必要があるが，三元数は乗法について閉じていないので，三元数を通して$\mathcal{J_{c}},\mathcal{M}$を得ることはできない，ということらしい． 実際に，三元数$a+bi+cj$を考え，$$ij=\alpha+\beta i+\gamma j$$を満たす実数$\alpha,\beta,\gamma$が存在するか調べると， $$(\alpha\gamma-\beta)+(\alpha+\beta\gamma)i+(\gamma^{2}+1)j=0$$ より，$$\gamma^{2}+1=0$$を満たす実数$\gamma$は存在しない．よって三元数は乗法について閉じていない． ということでわざわざ四元数を持ち出しているということだ．
これでクォータニオンを持ち込めば，3次元のフラクタルが得られることが分かったが，それをどうレンダリングすればいいのだろう． 2次元の場合は各ピクセルについて単に$\lim_{n \to \infty}f_{c}(z)$が収束するか否かを調べれば良かったが，3次元の場合はピクセルに向けて飛ばしたレイ上の全ての点をしらみつぶしに調べるわけにもいかない． フラクタルまでの距離がわかれば，スフィアトレーシングできるが，結論から言うと$\mathcal{J_{c}},\mathcal{M}$に対するDistance Functionは存在しない．
が，Distance Functionの近似なら求めることができるようだ．
The Orsay Notes [Douady and Hubbard, 1984] によると，点$z$におけるJulia集合$\mathcal{J_{c}}$に対するポテンシャル$G(z)$は以下で与えられる．
$$G(z)=\lim_{n \to \infty}\cfrac{1}{2^{n}}log|f_{c}^{n}(z)|$$
$\mathcal{J_{c}}$の定義から，$\mathcal{J_{c}}$に含まれる点$z$で$G(z)=0$となる． このポテンシャル$G(z)$を用いると，$\mathcal{J_{c}}$の外側のある点$z$から$\mathcal{J_{c}}$の境界までの距離$d(z)$のboundは以下のように導かれる．
$$\cfrac{\sinh{G(z)}}{2e^{G(z)}|G'(z)|}<d(z)<\cfrac{2\sinh{G(z)}}{|G'(z)|}$$
The Science of Fractal Images [Fisher, 1988] のAppendix Dにここの導出が詳しく書いてある． 実用上は$z$が$\mathcal{J_{c}}$に十分近いとして，
$$d(z)>\cfrac{|{f_{c}}^{n}(z)|}{2|{f_{c}^{n}}'(z)|}\log|f_{c}^{n}(z)|$$
とすることが多いようだ．これはポテンシャル$G(z)$をテイラー展開して1次近似してるのかと思ったが，厳密には違うっぽい． いずれにせよ下限がわかれば，その分だけレイを進めることができそうだ．
ここで，${f_{c}^{n}}'(z)$はどうやって求めるのだろう． 考えられるのは数値微分か自動微分．今回は双対数（Dual Number）を用いたフォワード型の自動微分を行った． 今回微分したい$f_{c}(z)$は入力，出力ともにクォータニオンなので，クォータニオンの双対数（Dual Quaternion）を考えることになり，結構ややこしそうだ．
まずはクォータニオン用の演算を定義しておく．以降のソースは全てGLSLフラグメントシェーダ．
// ---------------- quaternion ---------------- //
vec4 qAdd(vec4 q1, vec4 q2)
{
// return vec4(q1.x + q2.x, q1.yzw + q2.yzw);
return q1 + q2;
}
vec4 qSub(vec4 q1, vec4 q2)
{
// return vec4(q1.x - q2.x, q1.yzw - q2.yzw);
return q1 - q2;
}
vec4 qMul(vec4 q1, vec4 q2)
{
return vec4(q1.x * q2.x - dot(q1.yzw, q2.yzw), q2.x * q1.yzw + q1.x * q2.yzw + cross(q1.yzw, q2.yzw));
}
vec4 qConj(vec4 q)
{
return vec4(q.x, -q.yzw);
}
float qNorm(vec4 q)
{
// return sqrt(qMul(q, qConj(q)).x);
return length(q);
}
vec4 qInv(vec4 q)
{
return qConj(q) / pow(qNorm(q),2);
}
vec4 qDiv(vec4 q1, vec4 q2)
{
return qMul(q1, qInv(q2));
}
vec4 qPow(vec4 q, int n)
{
vec4 p = vec4(1.0, vec3(0.0));
for (int i = 0; i < n; ++i)
{
p = qMul(p, q);
}
return p;
}
次にDual Quaternion用の演算．ここら辺は定義に従って実装するだけ． 演算子のオーバーロードができないのが辛い．
// ---------------- dual ---------------- //
struct DualQ
{
vec4 q;
vec4 d;
};
DualQ dqAdd(DualQ dq1, DualQ dq2)
{
return DualQ(qAdd(dq1.q, dq2.q), qAdd(dq1.d, dq2.d));
}
DualQ dqSub(DualQ dq1, DualQ dq2)
{
return DualQ(qSub(dq1.q, dq2.q), qSub(dq1.d, dq2.d));
}
DualQ dqMul(DualQ dq1, DualQ dq2)
{
return DualQ(qMul(dq1.q, dq2.q), qAdd(qMul(dq1.d, dq2.q), qMul(dq1.q, dq2.d)));
}
DualQ dqDiv(DualQ dq1, DualQ dq2)
{
return DualQ(qDiv(dq1.q, dq2.q), qDiv(qSub(qMul(dq1.d, dq2.q), qMul(dq1.q, dq2.d)), qMul(dq2.q, dq2.q)));
}
DualQ dqPow(DualQ dq, int n)
{
DualQ dp = DualQ(vec4(1.0, vec3(0.0)), vec4(0.0, vec3(0.0)));
for (int i = 0; i < n; ++i)
{
dp = dqMul(dp, dq);
}
return dp;
}
一旦ここら辺を丁寧に実装すれば，低レベルのややこしい演算は抽象化できる． 双対数上で定義された関数$f$に双対数$x+\epsilon$をぶちこむと，$f(x+\epsilon)=f(x)+f'(x)\epsilon$となり，その虚部は常に導関数となっている． 双対数の虚部を，微分したい変数に関しては$1$，それ以外の変数に関しては$0$にしておき，あとは普通に関数に通すだけで良い． 普通にすごいな双対数．
今回で言うと，Julia集合$\mathcal{J_{c}}$に関しては$z$の虚部を$1$，Mandelbrot集合$\mathcal{M}$に関しては$c$の虚部を$1$にしておく．
つまづいたところは，${f_{c}^{n}}(z)$は入力，出力ともにクォータニオンなので，${f_{c}^{n}}'(z)$はヤコビ行列となりそうだが，$|{f_{c}^{n}}'(z)|$は何を意味するのか，ということだった． いろいろな文献をみて回ったが，結局納得のできる答えは見つからなかった．
ただ，僕のMacBookでレンダリングすることを考えると，毎イテレーションでヤコビ行列を計算したくないのもあり， Ray Tracing Quaternion Julia Sets on the GPU [Crane, 2015] に従って，クォータニオンの実部に関する微分のみ計算することにした．
この場合はDual Quaternionの虚部のクォータニオンの実部を$1$，虚部を$0$にすれば良い．（ややこしい…）
これでJulia集合における${f_{c}^{n}}(z)$，${f_{c}^{n}}'(z)$を同時に求めることができそうだ． ${f_{c}^{n}}(z)$，${f_{c}^{n}}'(z)$がわかれば，距離推定により，スフィアトレーシングできる． 距離推定のソースは以下．
// ---------------- fractals ---------------- //
float sdfJulia(vec4 z, vec4 c, int power, int numIterations, float escapeCriteria)
{
DualQ dzx = DualQ(z, vec4(1.0, 0.0, 0.0, 0.0));
DualQ dcx = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
for (int i = 0; i < numIterations; ++i)
{
// forward-mode automatic differentiation
dzx = dqAdd(dqPow(dzx, power), dcx);
if (qNorm(dzx.q) > escapeCriteria) break;
}
return (qNorm(dzx.q) * log(qNorm(dzx.q))) / (2 * qNorm(dzx.d));
}
float sdfMandelbrot(vec4 c, vec4 z, int power, int numIterations, float escapeCriteria)
{
DualQ dzx = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dcx = DualQ(c, vec4(1.0, 0.0, 0.0, 0.0));
for (int i = 0; i < numIterations; ++i)
{
// forward-mode automatic differentiation
dzx = dqAdd(dqPow(dzx, power), dcx);
if (qNorm(dzx.q) > escapeCriteria) break;
}
return (qNorm(dzx.q) * log(qNorm(dzx.q))) / (2 * qNorm(dzx.d));
}
あと必要なのは法線である．ポテンシャル$G(z)$または推定距離$d(z)$の微分が計算できれば良さそうだ． $$G(z)=\lim_{n \to \infty}\cfrac{1}{2^{n}}log|f_{c}^{n}(z)|$$より， $$\nabla G(z)=\cfrac{1}{2^{n}|f_{c}^{n}(z)|}\nabla\sqrt{f_{c}^{n}(z)^{T}f_{c}^{n}(z)}=\cfrac{{f_{c}^{n}}'(z)^{T}f_{c}^{n}(z)}{2^{n}|f_{c}^{n}(z)|^{2}}$$
これで距離推定と同時に法線も推定できることがわかった．
上式より法線に関してはヤコビ行列 ${f_{c}^{n}}'(z)$を計算したいので，DualQuaternionの虚部がmat4になる．mat4についてDualQuaternionの演算を定義するのは大変なので，複数の勾配を同時に追跡するか，以下のようにシンプルにchain ruleでヤコビ行列を追跡しても良さそう．
$$f_{c}^{n}(z)=(f_{c}^{n-1}(z))^{2}+c$$
より，
$${f_{c}^{n}}’(z)=2\begin{pmatrix}f_{c}^{n-1}(z)_{x} & -f_{c}^{n-1}(z)_{y} & -f_{c}^{n-1}(z)_{z} & -f_{c}^{n-1}(z)_{w}\\ f_{c}^{n-1}(z)_{y} & f_{c}^{n-1}(z)_{x} & 0 & 0\\ f_{c}^{n-1}(z)_{z} & 0 & f_{c}^{n-1}(z)_{x} & 0\\ f_{c}^{n-1}(z)_{w} & 0 & 0 & f_{c}^{n-1}(z)_{x}\end{pmatrix}{f_{c}^{n-1}}'(z)$$
これで$f_{c}^{n}(z)$の計算と同時に${f_{c}^{n}}'(z)$も計算できる． 最初以下のように普通に計算しようとしたが，クォータニオンの場合にはそもそも計算できなかった（?）ので，上式は手で微分して突っ込んだ．
$${f_{c}^{n}}’(z)=2{f_{c}^{n-1}}(z){f_{c}^{n-1}}'(z)$$
chain ruleでこんなに簡単に追跡できるなら，双対数を導入する必要性はなかった気もするが，これを機にいろいろ知れたし．これはこれで良かったと思う． 法線推定のソースは以下．
// ---------------- fractals ---------------- //
vec4 normalJulia(vec4 z, vec4 c, int power, int numIterations, float escapeCriteria)
{
DualQ dzx = DualQ(z, vec4(1.0, 0.0, 0.0, 0.0));
DualQ dzy = DualQ(z, vec4(0.0, 1.0, 0.0, 0.0));
DualQ dzz = DualQ(z, vec4(0.0, 0.0, 1.0, 0.0));
DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 1.0));
DualQ dcx = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dcy = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dcz = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 0.0));
for (int i = 0; i < numIterations; ++i)
{
// forward-mode automatic differentiation
dzx = dqAdd(dqPow(dzx, power), dcx);
dzy = dqAdd(dqPow(dzy, power), dcy);
dzz = dqAdd(dqPow(dzz, power), dcz);
dzw = dqAdd(dqPow(dzw, power), dcw);
if (qNorm(dzx.q) > escapeCriteria) break;
}
mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
return dzx.q * J;
}
vec4 normalMandelbrot(vec4 c, vec4 z, int power, int numIterations, float escapeCriteria)
{
DualQ dzx = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dzy = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dzz = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dzw = DualQ(z, vec4(0.0, 0.0, 0.0, 0.0));
DualQ dcx = DualQ(c, vec4(1.0, 0.0, 0.0, 0.0));
DualQ dcy = DualQ(c, vec4(0.0, 1.0, 0.0, 0.0));
DualQ dcz = DualQ(c, vec4(0.0, 0.0, 1.0, 0.0));
DualQ dcw = DualQ(c, vec4(0.0, 0.0, 0.0, 1.0));
for (int i = 0; i < numIterations; ++i)
{
// forward-mode automatic differentiation
dzx = dqAdd(dqPow(dzx, power), dcx);
dzy = dqAdd(dqPow(dzy, power), dcy);
dzz = dqAdd(dqPow(dzz, power), dcz);
dzw = dqAdd(dqPow(dzw, power), dcw);
if (qNorm(dzx.q) > escapeCriteria) break;
}
mat4 J = mat4(dzx.d, dzy.d, dzz.d, dzw.d);
return dzx.q * J;
}
これで推定された法線を使ってシェーディングできる．今回はシンプルなPhong反射モデルを使ってシェーディングした． 残りのソースは以下．
#version 150
// ================ definitions ================ //
// ---------------- rendering components ---------------- //
struct App
{
float time;
vec2 resolution;
};
struct Camera
{
mat4 viewMatrix;
mat4 projectionMatrix;
vec3 position;
};
struct PointLight
{
vec3 position;
vec4 ambientColor;
vec4 diffuseColor;
vec4 specularColor;
};
struct DirectionalLight
{
vec3 direction;
vec4 ambientColor;
vec4 diffuseColor;
vec4 specularColor;
};
struct Material
{
vec4 ambientColor;
vec4 diffuseColor;
vec4 specularColor;
vec4 emissionColor;
float shininess;
float refractiveIndex;
};
struct March
{
int numIterations;
float convergenceCriteria;
float finiteDifferenceEpsilon;
};
// ---------------- primitives ---------------- //
struct Line
{
vec3 position;
vec3 direction;
};
struct Sphere
{
float radius;
vec3 position;
};
struct Box
{
vec3 size;
vec3 position;
// mat3 orientation;
};
struct Intersection
{
vec3 position;
bool intersected;
};
// ---------------- scene components ---------------- //
struct FractalParams
{
int power;
int numIterations;
float escapeCriteria;
};
struct Fractal
{
FractalParams params;
Material material;
};
struct Scene
{
vec4 backgroundColor;
DirectionalLight light;
Sphere boundingSphere;
Fractal fractal;
};
// ================ variables ================ //
// ---------------- rendering components ---------------- //
uniform App uApp;
uniform Camera uCamera;
uniform March uMarch = March(300, 1e-4, 1e-4);
// ---------------- scene components ---------------- //
uniform Scene uScene = Scene(
vec4(vec3(0.05), 1.0),
DirectionalLight(vec3(-1.0, -1.0, -1.0), vec4(vec3(0.05), 1.0), vec4(vec3(0.8), 1.0), vec4(vec3(1.0), 1.0)),
Sphere(1.7, vec3(0.0)),
Fractal(FractalParams(2, 10, 2.0), Material(vec4(vec3(0.1), 1.0), vec4(vec3(0.8), 1.0), vec4(vec3(1.0), 1.0), vec4(vec3(0.0), 1.0), 64.0, 1.5))
);
// ---------------- varying ---------------- //
out vec4 oFragColor;
// ================ functions ================ //
// ---------------- utilities ---------------- //
vec2 linmap(vec2 in_val, vec2 in_min, vec2 in_max, vec2 out_min, vec2 out_max)
{
return (in_val - in_min) / (in_max - in_min) * (out_max - out_min) + out_min;
}
// ---------------- primitives ---------------- //
Intersection intersectionSphereLine(Sphere sphere, Line line)
{
vec3 difference = line.position - sphere.position;
float a = dot(line.direction, line.direction);
float b = 2 * dot(difference, line.direction);
float c = dot(difference, difference) - pow(sphere.radius, 2);
float d = pow(b, 2) - 4 * a * c;
float t = (-b - sqrt(d)) / (2 * a);
return Intersection(line.position + t * line.direction, d >= 0);
}
float sdfSphere(Sphere sphere, vec3 position)
{
return length(position - sphere.position) - sphere.radius;
}
vec3 normalSphere(Sphere sphere, vec3 position)
{
return normalize(position - sphere.position);
}
float sdfBox(Box box, vec3 position)
{
position = (position - box.position); // * box.orientation;
vec3 difference = abs(position) - box.size;
return length(max(difference, 0.0)) + min(max(difference.x, max(difference.y, difference.z)), 0.0);
}
vec3 normalBox(Box box, vec3 position, float finiteDifferenceEpsilon)
{
return normalize(vec3(
sdfBox(box, vec3(position.x + finiteDifferenceEpsilon, position.y, position.z)) - sdfBox(box, vec3(position.x - finiteDifferenceEpsilon, position.y, position.z)),
sdfBox(box, vec3(position.x, position.y + finiteDifferenceEpsilon, position.z)) - sdfBox(box, vec3(position.x, position.y - finiteDifferenceEpsilon, position.z)),
sdfBox(box, vec3(position.x, position.y, position.z + finiteDifferenceEpsilon)) - sdfBox(box, vec3(position.x, position.y, position.z - finiteDifferenceEpsilon))
));
}
// ---------------- reflection ---------------- //
vec4 phongReflection(vec3 surfaceNormal, vec3 eyeDirection, DirectionalLight light, Material material)
{
surfaceNormal = normalize(surfaceNormal);
eyeDirection = normalize(eyeDirection);
light.direction = normalize(light.direction);
vec4 ambientColor = light.ambientColor * material.ambientColor;
vec4 diffuseColor = light.diffuseColor * material.diffuseColor;
vec4 specularColor = light.specularColor * material.specularColor;
diffuseColor *= max(dot(-light.direction, surfaceNormal), 0.0);
specularColor *= pow(max(dot(reflect(light.direction, surfaceNormal), -eyeDirection), 0.0), material.shininess);
vec4 color = clamp(ambientColor + diffuseColor + specularColor + material.emissionColor, 0.0, 1.0);
return color;
}
float fresnelReflectance(vec3 light, vec3 normal, float refractiveIndex)
{
float specularReflectance = pow((1.0 - refractiveIndex) / (1.0 + refractiveIndex), 2.0);
return specularReflectance + (1.0 - specularReflectance) * pow(1.0 - dot(-light, normal), 5.0);
}
// ---------------- ray marching ---------------- //
vec4 rayMarching(App app, Scene scene, March march, Line ray)
{
Intersection intersection = intersectionSphereLine(scene.boundingSphere, ray);
if (intersection.intersected)
{
ray.position = intersection.position;
vec4 juliaType = 0.45 * cos(vec4(0.5, 3.9, 1.4, 1.1) + app.time * 0.15 * vec4(1.2, 1.7, 1.3, 2.5)) - vec4(0.3, 0.0, 0.0, 0.0);
vec4 criticalPoint = vec4(0.0);
for (int i = 0; i < march.numIterations; ++i)
{
float sdf = sdfJulia(vec4(ray.position, 0.0), juliaType, scene.fractal.params.power, scene.fractal.params.numIterations, scene.fractal.params.escapeCriteria);
// ray marching
ray.position += sdf * ray.direction;
// collision detection
if (abs(sdf) < march.convergenceCriteria)
{
vec3 surfaceNormal = normalize(normalJulia(vec4(ray.position, 0.0), juliaType, scene.fractal.params.power, scene.fractal.params.numIterations, scene.fractal.params.escapeCriteria).xyz);
vec4 fragColor = phongReflection(surfaceNormal, ray.direction, scene.light, scene.fractal.material);
return fragColor;
}
if (sdfSphere(scene.boundingSphere, ray.position) > 0) break;
}
}
return scene.backgroundColor;
}
void main()
{
vec2 fragCoord = linmap(gl_FragCoord.xy, vec2(0, 0), uApp.resolution, vec2(-1, -1), vec2(1, 1));
// why this does not work?
// vec3 rayDirection = normalize((inverse(uViewMatrix) * inverse(uProjectionMatrix) * vec4(vec3(fragCoord, 1.0), 0.0)).xyz);
vec3 rayDirection = normalize((inverse(mat3(uCamera.viewMatrix)) * inverse(mat3(uCamera.projectionMatrix)) * vec3(fragCoord, 1.0)).xyz);
Line ray = Line(uCamera.position, rayDirection);
oFragColor = rayMarching(uApp, uScene, uMarch, ray);
}
実際にGLSLでレンダリングしてみた結果が以下． ここまで来るのは正直結構大変だったが，鏡面反射してくれてるのを見ると，法線が推定できている感じがして嬉しい もっと格好良くレンダリングしてみたかったが，もっと勉強しないと格好良くは描けない．
3D Julia set
3D Mandelbrot set