Tilings
For completeness, let's include Golomb's picture again.
(Appeared in Mathematics Magazine, March 1985.)
It is elegant to interpret this in light of this tiling-based proof of the Pythagorean theorem.
Take $\vec v_0 = (a,b)$ and $\vec v_1 = (c,d)$, let $\Lambda$ be the lattice $\{m \vec v_0 + n \vec v_1 : m, n \in \mathbb Z\}$, and let $P$ be the parallelogram spanned by $\vec v_0$ and $\vec v_1$. The parallelograms $\vec v + P$ with $\vec v \in \Lambda$ form a $\Lambda$-periodic tiling of the plane.
We can also try to produce a $\Lambda$-periodic tiling of the plane by taking the rectangles $\vec v+R_+$ for all $\vec v \in \Lambda$, where $R_+ = [0,a] \times [0,d]$ is the big rectangle in Golomb's picture. If $a$, $b$, $c$, and $d$ are as in the picture then this results in some overlap, which can be removed by cutting away the rectangles $\vec v+R_-$ for all $\vec v \in \Lambda$, where $R_- = [-c, 0] \times [-b, 0].$ (The small rectangle in Golomb's picture is $\vec v_0 + \vec v_1 + R_-$, which is one of the rectangles that we cut away.)
So we can cover each point in the plane exactly once by tiling the plane $\Lambda$-periodically with $a \times d$ rectangles $R_+$ and then fixing the overlap by cutting away a $c \times b$ rectangle $R_-$ from each tile. Since we can also just tile it $\Lambda$-periodically with parallelograms $P$, it follows that the area of $P$ is the area of $R_+$ minus the area of $R_-$, i.e. $|P| = ad - cb$.
Note that strictly speaking there are a lot of special cases depending on the signs of $a$, $b$, $c$, $d$, and $ad-bc$. For example, if $a, b, d > 0$ and $c < 0$ then the $a \times d$ rectangles don't cover the entire plane, but the $-c \times b$ rectangles exactly fill in the remaining gaps, so $|P| = ad + (-c)b$. Indeed, when $a > b > 0$ and $c = -b, d = a$, this gives us precisely the picture in the tiling proof of the Pythagorean theorem.
So we have a picture which provides an interpretation for each term for Leibniz's formula for $2 \times 2$ determinants.
Three dimensions
We'd like to have a similar picture in higher dimensions. The parallelepiped $P$ spanned by $\vec u = (a_{0,0}, a_{0,1}, a_{0,2})$, $\vec v = (a_{1,0}, a_{1,1}, a_{1,2})$, and $\vec w = (a_{2,0}, a_{2,1}, a_{2,2})$ tiles space $\Lambda$-periodically, where $\Lambda$ is the lattice $\Lambda = \{ n \vec u + m \vec v + l \vec w : n, m, l \in \mathbb Z \}$.
To get a picture that explains Leibniz's formula, we hope that we can also tile space $\Lambda$-periodically using a combination of axis-aligned $a_{0,0} \times a_{1,1} \times a_{2,2}$, $a_{0,1} \times a_{1,2} \times a_{2,0}$, and $a_{0,2} \times a_{1,0} \times a_{2,1}$ bricks, with $a_{0,0} \times a_{1,2} \times a_{2,1}$, $a_{0,1} \times a_{1,0} \times a_{2,2}$, and $a_{0,2} \times a_{1,1} \times a_{2,0}$ bricks to remove the overlap. However, it's challenging to place these bricks appropriately.
It is easier to start by visualizing Laplace's formula, using the usual algebraic arguments to motivate the picture. The idea always boils down to the fact that the contribution of $\vec u$ to the volume of $P$ can be broken down into the contributions of its axis-aligned parts $\vec u_0 = (a_{0,0}, 0, 0)$, $\vec u_1 = (0, a_{0,1}, 0)$, and $\vec u_2 = (0,0,a_{0,2})$. To explain Laplace's formula, we want a visual interpretation of the fact that the (signed) volume of $P$ is the sum of the (signed) volumes of $P_0$, $P_1$ and $P_2$, where $P_k$ is the parallelepiped spanned by $\vec u_k, \vec v, \vec w$.
Hopefully we can produce a tiling which covers every point in space exactly once (counted with sign) by $\Lambda$-periodic translates of $P_0$, $P_1$, and $P_2$. Taking $\Lambda$-translates of $P_0$, $P_1$, and $P_2$ themselves doesn't do the trick, but taking $\Lambda$-translates of $P_0, \vec u_0 + P_1, \vec u_0 + \vec u_1 + P_2$ does.
Visually this corresponds to the following: If $S$ is the line segment from $\vec 0$ to $\vec u$ and $F$ is the parallelogram spanned by $\vec v$ and $\vec w$, then
$$P = S + F = \{ \vec s+ \vec f : \vec s \in S, \vec p \in F \}.$$
Deform $S$ into a concatenation of three axis-aligned line segments: $S_0$ from $\vec 0$ to $\vec u_0$, $S_1$ from $\vec u_0$ to $\vec u_0 + \vec u_1$, and $S_2$ from $\vec u_0 + \vec u_1$ to $\vec u_0 + \vec u_1 + \vec u_2 = \vec u$. This deforms $P$ into the union of the three regions
\begin{align*}
P_0 &= S_0 + F, \\
\vec u_0 + P_1 &= S_1 + F, \text{ and } \\
\vec u_0 + \vec u_1 + P_2 &= S_2 + F.
\end{align*}
The deformation doesn't change the total signed volume. The $\Lambda$-translates of these three regions cover each generic point in space in space exactly once (where being covered by a negatively-oriented parallelogram counts negatively). So we have the visual interpretation of Laplace's formula that we were after.
To get an interpretation of Leibniz's formula, we can repeat this with each of our three tiling elements, with $\vec v$ and then $\vec w$ playing the role of $\vec u$. We triple the number of tiling elements each time. So this gives us a $\Lambda$-periodic tiling with $27$ axis-aligned bricks, only six of which are nondegenerate, corresponding to the six terms of Leibniz's formula.
However, visualizing the process step by step is tedious and unsatisfactory. It's more instructive to work out what the end result is in one stroke.
Leibniz formula, directly
We might as well work in $n$ dimensions. Given a positively-oriented basis $\vec v_0, \dots, \vec v_{n-1}$ of $\mathbb R^n$, the parametrization $\phi: \mathbb R^n \to \mathbb R^n$ given by $\phi(t_0, \dots, t_{n-1}) = t_0 \vec v_0 + \dots + t_{n-1} \vec v_{n-1}$ sends the standard $\mathbb Z^n$-periodic tiling of the domain by unit $n$-cubes to a $\Lambda$-periodic tiling of the codomain by parallelograms.
In the previous section, we replaced the segment from $\vec 0$ to $\vec v_0$ with some axis-aligned segments, so we will do the same thing here. Explicitly, if $\vec v_0 = (a_{0,0}, a_{0,1}, \dots, a_{0,n-1})$ let $\gamma_0:\mathbb [0, 1] \to \mathbb R^n$ be the piecewise-linear path through the points
\begin{align*}
\gamma_0(0) &= (0, 0, \dots, 0),\\
\gamma_0(1/n) &= (a_{0,0}, 0, \dots, 0), \\
\gamma_0(2/n) &= (a_{0,0}, a_{0,1}, 0, \dots, 0), \\
&\dots \\
\gamma_0(n/n) &= (a_{0,0}, a_{0,1}, \dots, a_{0,n-1}) = \vec v_0,
\end{align*}
and extend it to a path $\gamma_0:\mathbb R \to \mathbb R^n$ via $\gamma_0(t+1) = \gamma_0(t) + \vec v_0$. Similarly for other $\gamma_k$.
Let $\tilde \phi:\mathbb R^n \to \mathbb R^n$ be given by
$$\tilde \phi(t_0, \dots, t_{n-1}) = \gamma_0(t_0) + \dots + \gamma_{n-1}(t_{n-1}).$$
Then roughly speaking $\tilde \phi$ is a "crinkly" piecewise axis-aligned parametrization of $\mathbb R^n$. It hits some points several times, but generic points are hit a total of once if you count being hit by a negatively-oriented piece as being hit $-1$ times.
Each piece of the domain of $\tilde \phi$ is a set of the form $\vec m + \vec \tau/n + C$ where
\begin{align*}
\vec m &\in \mathbb Z^n,\\
\vec \tau &\in \{0, \dots, n-1\}^n, \text{ and }\\
C &= [0, 1/n]^n.
\end{align*}
The image of this piece is an axis-aligned brick, spanned by $a_{0,\tau_0}\vec e_{\tau_0},\, a_{1,\tau_1}\vec e_{\tau_1}, \dots$ and $a_{n,\tau_n}\vec e_{\tau_n}$. The signed volume of this brick is $$\operatorname{sign}(\tau) \prod_k a_{k,\tau_k},$$ where we think of $\tau$ as a function $\{0, \dots, n-1\} \to \{0, \dots, n-1\}$ and take $\operatorname{sign}(\tau)$ to be the sign of $\tau$ if $\tau$ is a permutation and $0$ otherwise.
This gives us a $\Lambda$-periodic tiling of $\mathbb R^n$ which hits each point a (signed) total of once, with each period made up of $n^n$ tiles, of which only $n!$ are nondegenerate. The total signed volume of the tiles in each period is
$$\sum_{\tau \in S_n} \operatorname{sign}(\tau) \prod_k a_{k,\tau_k}.$$
Since we can also $\Lambda$-periodically tile $\mathbb R^n$ with translates of the parallelepiped $P$, we have
$$|P| = \sum_{\tau \in S_n} \operatorname{sign}(\tau) \prod_k a_{k,\tau_k}.$$
Remarks
One can probably make this argument more formal by noting that $\phi, \tilde \phi:\mathbb R^n \to \mathbb R^n$ induce maps $\phi', \tilde \phi':\mathbb R^n / \mathbb Z^n \to \mathbb R^n / \Lambda$ between toruses, and these maps are homotopic.
Probably another way to make this argument more formal is to note that $T^{-1} \tilde \phi( T \vec t) \to \phi(\vec t)$ as $T \to \infty$, so the way these maps affect volume should be the same.
Strictly speaking, the picture this produces in two dimensions is not quite the one described in the first section. We get $\Lambda$-translates of the rectangles $R_+ = [c,a+c] \times [0,d]$ and $R_- = [a,a+c] \times [0,b]$, which is a translation of the original picture by $(c,0)$. A different choice of $\gamma_1$ would produce the original picture.
To visualize the crinkly map $\tilde \phi$ more clearly in a specific example, it helps to tweak it slightly so that the degenerate parallelepipeds are very thin but not quite degenerate.
In the two-dimensional case, it is possible to pick signs for the coefficients of our matrix such that all the nondegenerate bricks have positive volume, i.e. such that all terms in the Leibniz formula are positive. This gives a tiling with no overlap, which makes the argument less technical. Unfortunately, this is not possible generically for $3 \times 3$ matrices or above. The best one can do in the $3 \times 3$ case is a matrix with one zero entry, e.g. with signs $$\begin{pmatrix}+ & + & + \\ + & + & - \\ - & + & 0 \end{pmatrix}.$$ This gives four nondegenerate bricks rather than the full six. It would be interesting to produce a picture or a model of this situation.
