Painting With Math
From Shaders to Scenes
Raymarching SDFs Friendly Optimizations
Who am I?
-
Name: Diego Fonseca
-
Full-Time Student at NVCC
-
Visiting Student at Stanford
-
4+ years of experience developing custom rendering engines
Big Boi Stuff
I Work(ed) On
Main MacOS Programmer
(Since 2024)
Lead Engine Programmer
(Since 2021)
Let's Begin!
Inspired By Inigo Quilez
Let's First Discuss Raycasting
What Is It?
-
Rays are shot from a camera's viewpoint through each pixel of the screen
-
Determines what object, if any, the ray hits first
-
The color of that object is used for that specific pixel
Variations of Raycast Algorithms
}
-
Raytracing
-
Raymarching
How Does Raymarching
Work?
Let's see the simplest example
Sphere Equation
(x-c_x)^2+(y-c_y)^2+(z-c_z)^2=R^2
\text{Let $c_x$, $c_y$, $c_z$, and $R$ be constants.}
What Does This Have to do With Rays?
\text{Let }r(t) =
\begin{pmatrix}
x(t) & y(t) & z(t)
\end{pmatrix}^T
\text{Let }c =
\begin{pmatrix}
c_x & c_y & c_z
\end{pmatrix}^T
(x-c_x)^2+(y-c_y)^2+(z-c_z)^2=R^2
\sqrt{(x-c_x)^2+(y-c_y)^2+(z-c_z)^2}=R
What Does This Have to do With Rays?
||r(t)-c||-R=0
Continue..
\text{Let } f(r(t))=||r(t)-c||-R
\text{So what is $f(r(t))$?}
I got lazy so I'mma freestyle the explanation
That's why it's called
"Signed Distance Function"
What if we Didn't Intersect The Sphere's Surface?
Well
We reiterate!
\text{For now on $p = r(t)$}
Ok Then?
What's the Math?
I'm Glad You Asked!
\frac{dt}{ds}=f(p)
dt=f(p)ds
Where
\int_{0}^{t_f}{dt}=\int_{0}^{N}{f(p)ds}
t_f\approx\sum_{k=0}^{N-1}{f(p_k)\Delta{s_k}}, \Delta{s} = 1 \ \text{(For this case)}
#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}#define N 80
// Sphere constants
const vec3 c = vec3(0., 0., 0.)
const float R = 1.0;
float scene(vec3 p) {
return length(p - c) - R; // function `f(p)`
}
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 uv = (fragCoord * 2. - iResolution.xy) / iResolution.y;
vec3 ro = vec3(0, 0, -3); // Ray Origin
vec3 rd = normalize(vec3(uv, -1)); // Ray Direction
float t = 0.; // Total distance travelled
for(int i = 0; i < N; i++){
vec3 p = ro + rd * t;
float dt = scene(p); // ds = 1
t += dt;
if(dt < 0.001) { // Exit Epsilon
fragColor = vec4(1., 0., 0., 1.); // If hit, then red
return;
}
}
fragColor = vec4(0., 0., 0., 1.);
}How do we Express Adding Multiple Objects in a Scene?
\text{$scene(p) = f_1(p) \oplus f_2(p) \oplus ... \oplus f_n(p)$}
\text{$scene(p) = f_1(p) \oplus (f_2(p) \oplus ...\oplus f_n(p))$}
\text{$scene(p) = (f_1(p) \oplus f_2(p)) \oplus ... \oplus f_n(p)$}
X \otimes Y = X + Y
Understanding Multiplication
\log xy = \log x + \log y
\text{Let } l_b = log_b
\text{such that } l_b \in \mathbb{R_{>0}}
\mathbb{R_{>0}}
\mathbb{R}
l_b
l_b^{-1}
(b > 1)
\mathbb{R_{>0}}
\mathbb{R}
l_b
l_b^{-1}
l_b \big(l_b^{-1}(X) \cdot l_b^{-1}(Y) \big) = X \otimes Y
\log_b (b^X \cdot b^Y) = X \otimes Y
\log_b (b^{X+Y}) = X \otimes Y
X+Y = X \otimes Y
X \oplus_b Y = l_b\big(l_b^{-1}(X) + l_b^{-1}(Y)\big)
X \oplus_b Y = \log_b\big(b^X + b^Y\big)
X \oplus_b Y = \lim_{b \to 0^+} \log_b\big(b^X + b^Y\big) = \min(X,Y)
OR
X \oplus_b Y = \lim_{b \to \infin^+} \log_b\big(b^X + b^Y\big) = \max(X,Y)
Squeeze Theorem!
\text{$scene(p) = f_1(p) \oplus f_2(p) \oplus f_3(p)$}
Visualization
Visualization
What was the point of Tropical Algebra
The Point (Partially)
-\frac{1}{k}\big((-kX) \oplus_e (-kY)\big)
-\frac{1}{k}\big((-kX) \oplus_e (-kY)\big)
-\frac{1}{k}\big((-kX) \oplus_e (-kY)\big)
-\frac{1}{k}\ln(e^{-kX}+e^{-kY})
\text{Do you see it?}
LogSumExp Function
-\frac{1}{k}\ln(e^{-kX}+e^{-kY})
LogSumExp Function
My Demo
Creating Custom Formulas
(For Art Purposes)
Creating Custom Formulas
(For Art Purposes)
Creating Custom Formulas
(For Art Purposes)
Character Modeling
Character Modeling
Thanks For Listening!
(End of Part 1)
Any Questions?
Any More Questions?
C'mon. Don't Be Shy.
Understanding
the Series
t_f\approx\sum_{k=0}^{N-1}{f(p_k)\Delta{s_k}}
For Optimization Purposes
t_f=\sum_{k=0}^{N-1}{f(p_k)\Delta{s_k}}=t_\alpha \quad (N \rightarrow \infin)
\text{Assuming $t_f$ will intersect a surface}
\text{where } t_\alpha \in \mathbb{R}
\text{such that } \sum^{\infin}_{k=0}a_k = t_\alpha
\text{Riemann Rearrangement Theorem}
\sum_{k=0}^{\infin} a_k: \text{ cond. conv.,} \quad t_\beta \in \mathbb{R}
\sum_{k=0}^{\infin} a_k = t_\alpha \xrightarrow{\text{rearrange}} \sum_{k=0}^{\infin} b_k = t_\beta
\text{Riemann Rearrangement Theorem}
\sum_{k=0}^{\infin} a_k: \text{ cond. conv.,} \quad t_\beta \in \mathbb{R}
\sum_{k=0}^{\infin} a_k = t_\alpha \xrightarrow{\text{rearrange}} \sum_{k=0}^{\infin} b_k = t_\beta
\text{Riemann Rearrangement Theorem}
\sum_{k=0}^{\infin} a_k: \text{ cond. conv.,} \quad t_\beta \in \mathbb{R}
\sum_{k=0}^{\infin} a_k = t_\alpha \xrightarrow{\text{rearrange}} \sum_{k=0}^{\infin} b_k = t_\beta
\text{Riemann Rearrangement Theorem}
a_0+a_1+a_2+a_3+a_4+...
a_2+a_0+a_4+a_1+a_3+...
\text{Riemann Rearrangement Theorem}
a_0+a_1+a_2+\textcolor{red}{a_3}+a_4+...
a_2+a_0+a_4+a_1+\textcolor{red}{a_3}+...
\text{Appear Once}
\text{Tying it all together}
t_f=\sum_{k=0}^{N-1}{f(p_k)\Delta{s_k}}=t_\alpha
\text{Let } K = [0,N-1] \quad \therefore \sum_{k\in K}{f(p_k)\Delta{s_k}}=t_\beta
\text{Tying it all together}
\text{Let } I,J \subset K
\textcircled{1} \quad I \cup J = K
\textcircled{2} \quad I \cap J = \emptyset
\textcircled{3} \quad card(I) = card(J)
\text{Tying it all together}
\text{Let } I,J \subset K
\textcircled{1} \quad I \cup J = K
\textcircled{2} \quad I \cap J = \emptyset
\textcolor{red}{\textcircled{3} \quad card(I) = card(J)}
\text{Over simplification}
t_\beta = \gamma(N)f(x_0) + \sum_{k\in K} f(p_k)\Delta{s}_k
\text{Goal: } \quad card(K) \equiv 2
\text{such that now } K = [1-\gamma(N), N-1]
\gamma(N) = \begin{cases}
1, \text{if}\ N\ \text{is odd} \\
0, \text{else}
\end{cases}
\text{Over simplification}
\text{Tying it all together}
\text{Let } I,J \subset K
\textcircled{1} \quad I \cup J = K
\textcircled{2} \quad I \cap J = \emptyset
\textcircled{3} \quad card(I) = card(J)
\text{Over simplification}
t_\beta = \gamma(N)f(x_0) + \sum_{k\in K} f(p_k)\Delta{s}_k
t_\beta = \gamma(N)f(x_0) + \sum_{i\in I} f(p_i)\Delta{s}_i + \sum_{j\in J} f(p_j)\Delta{s}_j
\text{(The real equation involves \textcolor{red}{tuples} but this is simpler and holds true)}
#define MAX_STEPS 100
float raymarch(vec3 ray_origin, vec3 ray_direction) {
float offset_size = 10.0;
float t = 0.0;
int iter = 0;
float t_j = cam_block.far;
float min_dist = cam_block.far;
int i = int((MAX_STEPS & 1) == 0);
#if (MAX_STEPS & 1) == 0
{
t = scene(ray_origin + t * ray_direction, 2.0, offset_size);
if(t < cam_block.near) {
return t;
}
}
#endif
while(i <= (MAX_STEPS >> 1)) {
vec3 p_i = ray_origin + t * ray_direction;
float dist_i = scene(p_i, 2.0, offset_size);
if(dist_i < cam_block.near) {
return t;
}
if(t > cam_block.far) {
return -1.0;
}
if(dist_i - min_dist > MIN_GROWTH) {
vec3 p_j = ray_origin + t_j * ray_direction;
float dist_j = scene(p_j, 2.0, offset_size);
if(abs(dist_i - dist_j) <= abs(t_j - t)
&& (dist_i + dist_j) >= abs(t_j - t)) {
return -1.0;
}
t_j -= dist_j;
}else {
t += dist_i;
dist_i = scene(ray_origin + t * ray_direction, 2.0, offset_size);
if(dist_i < cam_block.near) {
return t;
}
if(t > cam_block.far) {
return -1.0;
}
}
min_dist = min(min_dist, dist_i);
t += dist_i;
i++;
}
return t;
}How Does it Work?
Thanks For Listening!