diff --git a/.gitignore b/.gitignore
index 3ef5ac2..181daa2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,5 @@ build/
 /output/
 /.vscode/
 /.pixi/
+*.aux
+*.log
diff --git a/docs/fmb-math.qmd b/docs/fmb-math.qmd
new file mode 100644
index 0000000..75ce728
--- /dev/null
+++ b/docs/fmb-math.qmd
@@ -0,0 +1,190 @@
+\newcommand{\L}{\mathcal{L}}
+\newcommand{\N}{\mathbb{N}}
+\newcommand{\R}{\mathbb{R}}
+\newcommand{\SE}{\mathrm{SE}}
+\newcommand{\argmax}{\operatorname*{argmax}}
+
+# Mathematical Description of Fuzzy Metaball Rendering
+
+## Notation
+
+$\N_n$ denotes the set $\{1, 2, \ldots, n\}$. $\SE(d)$ denotes the group of
+$d$-dimensional rigid transformations (hence $\SE(3)$ is the group of poses).
+$S^d$ denotes the $d$-dimensional sphere, i.e. the set of $d + 1$ dimensional
+unit vectors.
+
+The notation $\mathrm{symbol} := \mathrm{expression}$ offers the definition of "symbol" in terms of the "expression."
+
+Greek letters are unknowns, Latin letters are known
+
+Writing $(\tau, \rho) := \pi $ for $\pi \in \SE(3)$ means that $\pi$'s translational
+component is $\tau$ and it's rotational component is $\rho$. For $v \in \R^d$
+we write $\pi\cdot v$ to denote $v$ transformed by $\pi$.
+
+| Symbol | Type | Meaning |
+|:------:|:----:|:--------|
+| $N$    | $\N$ | number of fuzzy metaballs |
+| $T$    | $\N$ | number of time steps (frames) |
+| $H$    | $\N$ | height of each observed frame in pixels |
+| $W$    | $\N$ | width of each observed frame in pixels |
+| $i$    | $\N_H$ | index ranging over the rows in each observed frame |
+| $j$    | $\N_W$ | index ranging over the columns in each observed frame |
+| $k$    | $\N_N$ | index ranging over the fuzzy metaballs |
+| $t$    | $\N_T$ | index ranging over frames |
+| $v_{ij}$ | $S^2$ | direction of the ray corresponding to pixel $(i, j)$ in the camera frame |
+| $\pi^{(t)}$ | $\SE(3)$ | camera pose in world frame at frame $t$ |
+| $\mu_k$ | $\R^3$ | the mean vector of the $k$-the fuzzy metaball |
+| $\Sigma_k$ | $\R^{3\times 3}_{\succ 0}$ | the covariance matrix of the $k$-the fuzzy metaball |
+| $\lambda_k$ | $\R$ | log-weight of the $k$-the fuzzy metaball |
+| $d_{ijk}^{(t)}$    | $\R$ | intersection depth of ray corresponding to pixel $(i, j)$ at frame $t$ with the $k$-th fuzzy metaball |
+| $q_{ijk}^{(t)}$    | $\R$ | the value of the quadratic form of $k$-th fuzzy metaball at the intersection of ray $(i, j)$ at frame $t$ |
+| $w_{ijk}^{(t)}$    | $\R$ | depth blending weight of the $k$-th fuzzy metaball for pixel $(i, j)$ at frame $t$ |
+| $\bar{d}_{ij}^{(t)}$    | $\R$ | depth value of pixel $(i, j)$ at frame $t$ |
+| $c_{ij}^{(t)}$    | $[0, 1]$ | confidence value of pixel $(i, j)$ at frame $t$ |
+| $\L_c$    | $\R$ | contour loss  (a function of $c_{ij}^{(t)}$) |
+| $\L_d$    | $\R$ | depth loss  (a function of $d_{ij}^{(t)}$) |
+
+
+## Forward Pass
+
+\begin{align*}
+%\argmax_{d \in \R} (\mu_k - \Delta x^{(t)} - dv_{ij}^{(t)})^\top
+%\Sigma^{-1} (\mu_k - \Delta x^{(t)} - dv_{ij}^{(t)}) \\
+(\tau^{(t)}, \rho^{(t)}) &:= \pi^{(t)} \\
+v_{ij}^{(t)} &:= \rho^{(t)}v_{ij} \\
+d_{ijk}^{(t)} &:= \frac{(\mu_k - \tau^{(t)})^\top\Sigma_k^{-1}v_{ij}^{(t)}}
+    {{v_{ij}^{(t)}}^\top\Sigma_k^{-1}v_{ij}^{(t)}} \\
+q_{ijk}^{(t)} &:= (\mu_k - \tau^{(t)} - d_{ijk}^{(t)}v_{ij}^{(t)})^\top
+    \Sigma_k^{-1} (\mu_k - \tau^{(t)} - d_{ijk}^{(t)}v_{ij}^{(t)}) \\
+\tilde{w}_{ijk}^{(t)} &:=
+    \exp\left(\beta_1(-q_{ijk}^{(t)}/2 + \lambda_k) - \beta_2 d_{ijk}^{(t)}\right) \\
+w_{ijk}^{(t)} &:= \frac{\tilde{w}_{ijk}^{(t)}}{\sum_{k=1}^N \tilde{w}_{ijk}^{(t)}} \\
+\bar{d}_{ij}^{(t)} &:= \sum_{k=1}^N w_{ijk}^{(t)} d_{ijk}^{(t)} \\
+c_{ij}^{(t)} &:= 1 - \exp\left(\sum_{k=1} \exp(q_{ijk}^{(t)}/2 - \lambda_k)\right)
+\end{align*}
+
+![Dependency graph of the variables in the forward pass using plate
+notation](plate_diagram.png){width=540}
+
+## Depth Backward Pass
+
+\begin{align*}
+\frac{\partial \L_d}{\partial \tilde{w}_{ijk}^{(t)}} &= 
+    \frac{\partial \L_d}{\partial \bar{d}_{ij}^{(t)}}
+    \frac{\partial \bar{d}_{ij}^{(t)}}{\partial \tilde{w}_{ijk}^{(t)}}\\
+\frac{\partial \L_d}{\partial q_{ijk}^{(t)}} &= 
+    \frac{\partial \L_d}{\partial \tilde{w}_{ijk}^{(t)}}
+    \frac{\partial \tilde{w}_{ijk}^{(t)}}{\partial q_{ijk}^{(t)}}\\
+\frac{\partial \L_d}{\partial d_{ijk}^{(t)}} &= 
+    \frac{\partial \L_d}{\partial q_{ijk}^{(t)}}
+    \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+    + \frac{\partial \L_d}{\partial \tilde{w}_{ijk}^{(t)}}
+    \frac{\partial \tilde{w}_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+    + \frac{\partial \L_d}{\partial \bar{d}_{ij}^{(t)}}
+    \frac{\partial \bar{d}_{ij}^{(t)}}{\partial d_{ijk}^{(t)}} \\
+\frac{\partial \L_d}{\partial \tau^{(t)}} &=
+    \sum_{i,j}\sum_k \frac{\partial \L_d}{\partial q_{ijk}^{(t)}}
+    \frac{\partial q_{ijk}^{(t)}}{\partial \tau^{(t)}}
+    + \sum_{i,j}\sum_k \frac{\partial \L_d}{\partial d_{ijk}^{(t)}}
+    \frac{\partial d_{ijk}^{(t)}}{\partial \tau^{(t)}} \\
+\frac{\partial \L_d}{\partial v_{ij}^{(t)}} &=
+    \sum_k \frac{\partial \L_d}{\partial q_{ijk}^{(t)}}
+    \frac{\partial q_{ijk}^{(t)}}{\partial v_{ij}^{(t)}}
+    + \sum_k \frac{\partial \L_d}{\partial d_{ijk}^{(t)}}
+    \frac{\partial d_{ijk}^{(t)}}{\partial v_{ij}^{(t)}} \\
+\frac{\partial \L_d}{\partial \rho^{(t)}} &=
+    \sum_{i,j} \frac{\partial \L_d}{\partial v_{ij}^{(t)}}
+           \frac{\partial v_{ij}^{(t)}}{\partial \rho^{(t)}} \\
+\frac{\partial \L_d}{\partial \lambda_k} &= 
+    \sum_t \sum_{i, j}
+    \frac{\partial \L_d}{\partial \tilde{w}_{ijk}^{(t)}}
+    \frac{\partial \tilde{w}_{ijk}^{(t)}}{\partial \lambda_k} \\
+\frac{\partial \L_d}{\partial \mu_k} &= 
+    \sum_t \sum_{i, j} \frac{\partial \L_d}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial \mu_k}
+    + \sum_t \sum_{i, j} \frac{\partial \L_d}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial \mu_k} \\
+\frac{\partial \L_d}{\partial \Sigma_k} &= 
+    \sum_t \sum_{i, j} \frac{\partial \L_d}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial \Sigma_k}
+    + \sum_t \sum_{i, j} \frac{\partial \L_d}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial \Sigma_k} \\
+\end{align*}
+
+
+## Confidence Backward Pass
+
+\begin{align*}
+\frac{\partial \L_c}{\partial d_{ijk}^{(t)}} &= 
+    \frac{\partial \L_c}{\partial q_{ijk}^{(t)}}
+    \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}\\
+\frac{\partial \L_c}{\partial q_{ijk}^{(t)}} &= 
+    \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+    \frac{\partial c_{ij}^{(t)}}{\partial q_{ijk}^{(t)}}\\
+\frac{\partial \L_c}{\partial \tau^{(t)}} &= 
+    \sum_{i,j}\sum_k \frac{\partial \L_c}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial \tau^{(t)}}
+    + \sum_{i,j}\sum_k \frac{\partial \L_c}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial \tau^{(t)}}\\
+\frac{\partial \L_c}{\partial v_{ij}^{(t)}} &= 
+    \sum_{k} \frac{\partial \L_c}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial v_{ij}^{(t)}}
+    + \sum_{k} \frac{\partial \L_c}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial v_{ij}^{(t)}} \\
+\frac{\partial \L_c}{\partial \rho^{(t)}} &= 
+    \sum_{i,j} \frac{\partial \L_c}{\partial v_{ij}^{(t)}}
+        \frac{\partial v_{ij}^{(t)}}{\partial \rho^{(t)}} \\
+\frac{\partial \L_c}{\partial \lambda_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+        \frac{\partial c_{ij}^{(t)}}{\partial \lambda_k} \\
+\frac{\partial \L_c}{\partial \mu_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial \mu_k}
+    + \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial \mu_k} \\
+\frac{\partial \L_c}{\partial \Sigma_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial q_{ijk}^{(t)}}
+        \frac{\partial q_{ijk}^{(t)}}{\partial \Sigma_k}
+    + \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial \Sigma_k}
+\end{align*}
+
+We can simplify the backward pass by inlining those derivatives that depend on
+fixed $i$, $j$, and $t$. We get:
+
+\begin{align*}
+\frac{\partial \L_c}{\partial \tau^{(t)}} &= 
+    \sum_{i,j}\sum_k \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+        \frac{\partial c_{ij}^{(t)}}{\partial q_{ijk}^{(t)}} \left(
+            \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+            \frac{\partial d_{ijk}^{(t)}}{\partial \tau^{(t)}}
+            + \frac{\partial q_{ijk}^{(t)}}{\partial \tau^{(t)}}
+      \right)\\
+\frac{\partial \L_c}{\partial \rho^{(t)}} &= 
+    \sum_{i,j} \sum_{k} \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+    \frac{\partial c_{ij}^{(t)}}{\partial q_{ijk}^{(t)}} \left(
+        \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+        \frac{\partial d_{ijk}^{(t)}}{\partial v_{ij}^{(t)}}
+        + \frac{\partial q_{ijk}^{(t)}}{\partial v_{ij}^{(t)}}
+    \right) 
+        \frac{\partial v_{ij}^{(t)}}{\partial \rho^{(t)}} \\
+\frac{\partial \L_c}{\partial \lambda_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+        \frac{\partial c_{ij}^{(t)}}{\partial \lambda_k} \\
+\frac{\partial \L_c}{\partial \mu_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+        \frac{\partial c_{ij}^{(t)}}{\partial q_{ijk}^{(t)}} \left(
+        \frac{\partial q_{ijk}^{(t)}}{\partial \mu_k}
+        + \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+          \frac{\partial d_{ijk}^{(t)}}{\partial \mu_k}
+    \right) \\
+\frac{\partial \L_c}{\partial \Sigma_k} &= 
+    \sum_{t}\sum_{i,j} \frac{\partial \L_c}{\partial c_{ij}^{(t)}}
+        \frac{\partial c_{ij}^{(t)}}{\partial q_{ijk}^{(t)}} \left(
+        \frac{\partial q_{ijk}^{(t)}}{\partial \Sigma_k}
+        + \frac{\partial q_{ijk}^{(t)}}{\partial d_{ijk}^{(t)}}
+          \frac{\partial d_{ijk}^{(t)}}{\partial \Sigma_k}
+    \right) \\
+\end{align*}
+
+## 
diff --git a/docs/make_fig.sh b/docs/make_fig.sh
new file mode 100755
index 0000000..4ddc240
--- /dev/null
+++ b/docs/make_fig.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+pdflatex plate_diagram.tex
+convert -verbose -density 300 -trim plate_diagram.pdf -size 1080x1080 -quality 100 -flatten -sharpen 0x1.0 plate_diagram.png
diff --git a/docs/plate_diagram.png b/docs/plate_diagram.png
new file mode 100644
index 0000000..b621959
Binary files /dev/null and b/docs/plate_diagram.png differ
diff --git a/docs/plate_diagram.tex b/docs/plate_diagram.tex
new file mode 100644
index 0000000..2c0ea13
--- /dev/null
+++ b/docs/plate_diagram.tex
@@ -0,0 +1,77 @@
+\documentclass[border=10pt]{standalone}
+\usepackage{tikz}
+\usepackage{amsmath,amssymb}
+
+\usetikzlibrary{arrows.meta, positioning, fit, backgrounds}
+
+\begin{document}
+% if plate p binds index i and variable v appears in p, then v should vary over
+% i
+\begin{tikzpicture}[
+    node distance=1cm and 1cm,
+    >={Stealth[length=2mm]},
+    every node/.style={font=\small},
+    every edge/.style = {draw, -latex}
+]
+
+% Define coordinates for nodes
+\node (lambdak) at ( -1, 5) {$\lambda_k$};
+\node (muk)  at (  0, 5) {$\mu_k$};
+\node (Sigmak)     at (  1, 5) {$\Sigma_k$};
+
+\node (vij)     at (2.5, 4) {$v_{ij}^{(t)}$};
+\node (rho)     at (  4, 4) {$\rho^{(t)}$};
+
+\node (tau)     at (  4, 3.5) {$\tau^{(t)}$};
+
+\node (qijk)    at (  0, 2.5) {$q_{ijk}^{(t)}$};
+\node (dijk)    at (  1, 2.5) {$d_{ijk}^{(t)}$};
+
+\node (wijk)    at (  0, 1) {$\tilde{w}_{ijk}^{(t)}$};
+
+\node (cij)     at ( -1, -0.2) {$c_{ij}^{(t)}$};
+\node (dijbar)  at (  1, -0.2) {$\bar{d}_{ij}^{(t)}$};
+
+\node (Lc)      at ( -1, -1.4) {$\mathcal{L}_c$};
+\node (Ld)      at (  1, -1.4) {$\mathcal{L}_d$};
+
+
+\draw (Sigmak.south) edge (qijk.north) edge (dijk.north)
+      (muk.south) edge (qijk.north) edge (dijk.north)
+      (vij.south) edge (qijk.north) edge (dijk.north)
+      (rho.west) edge (vij.east)
+      (tau.west) edge (qijk.north) edge (dijk.north)
+      (lambdak.south) edge (wijk.north) edge (cij.north)
+      (qijk.south) edge (wijk.north) edge (cij.north)
+      (dijk.south) edge (wijk.north) edge (dijbar.north)
+      (dijk.west) edge (qijk.east)
+      (wijk.south) edge (dijbar.north)
+      (cij) edge (Lc)
+      (dijbar) edge (Ld);
+
+% Inner plate k
+\begin{scope}[on background layer]
+    \node[draw, rectangle, rounded corners=2pt, 
+          inner sep=4pt,
+          fit=(Sigmak)(muk)(lambdak)(wijk)(dijk),
+          label={[anchor=south east, font=\small]south east:$k$}] (plate_k) {};
+\end{scope}
+
+% Middle plate (i,j)
+\begin{scope}[on background layer]
+    \node[draw, rectangle, rounded corners=2pt,
+          inner sep=4pt,
+          fit=(wijk)(dijk)(vij)(cij)(dijbar),
+          label={[anchor=south east, font=\small]south east:$(i,j)$}] (plate_ij) {};
+\end{scope}
+
+% Outer plate t
+\begin{scope}[on background layer]
+    \node[draw, rectangle, rounded corners=2pt,
+          inner sep=4pt,
+          fit=(plate_ij)(rho)(tau),
+          label={[anchor=south east, font=\small]south east:$t$}] (plate_t) {};
+\end{scope}
+
+\end{tikzpicture}
+\end{document}