The geometric approach works.

Let’s compute the volume of the $2n$ dimensional ball, $D^{2n}$, in two ways. One way is extremely clever but has been known for centuries and provides interesting insights: it’s based on Liouville’s trick. Specifically, we will compute two integrals in polar coordinates, one of which is the volume of the ball and the other of which reduces to a product of one-dimensional integrals. Both integrands will depend (at most) on the radial coordinate $r$, which lets us separate out the surface area of the boundary of the ball as a common factor. Write this surface area as $S_{2n-1}$.

There’s essentially just one way to do this trick: integrate $\exp(-r^2)$. Its integral over $\mathbb{R}^{2n}$ equals

$$S_{2n-1} \int_0^\infty {\exp\left(- r^2 \right) r^{2n-1} dr}.$$

However, because $r^2 = x_1^2 + x_2^2 + \ldots + x_{2n}^2$, the integrand (in *Cartesian* coordinates $\left( x_1, x_2, \ldots, x_{2n} \right)$) factors as $\exp\left(-r^2 \right) = \exp\left(-x_1^2 \right) \cdots \exp\left(-x_{2n}^2 \right)$, each of which must be integrated from $-\infty$ to $+\infty$. Whence

$$S_{2n-1} \int_0^\infty {\exp \left(- r^2 \right) r^{2n-1} dr} = \left( \int_{- \infty}^ \infty {\exp \left( -x^2 \right) dx} \right) ^{2n}.$$

I will call the left hand integral $\tfrac{1}{2} \Gamma \left(n \right)$, because that is what it is (as a simple change of variables shows). In the same notation, $\Gamma \left(1/2 \right) = \int_{-\infty}^\infty {\exp\left(-x^2 \right) dx}$. Algebraic re-arrangement of the foregoing yields the volume of $D^{2n}$ as

$$|D^{2n} | = S_{2n - 1} \int_0^1 {r^{2n - 1} dr} = \frac{{S_{2n - 1} }}
{{2n}} = \frac {\Gamma \left(1/2 \right)^{2n}} { n \Gamma \left(n \right) }.$$

That was the first method: the result is a familiar one, but has been left expressed in a way that better reveals its origins in polar and Cartesian integration.

The next way to compute the ball's volume is, I believe, new. It is inspired by Cavalieri’s Principle: the idea that you can shift slices of a solid around without changing the volume of the solid. The generalization is to move two-dimensional slices around *and to change their shapes while you do so,* but in a way that does not change their areas. It follows that the new solid has the same (hypervolume) as the original, although it might have a completely different shape.

We will compute the volume of a region $Q_n$ in $\mathbb{R}^{2n}$. It is conveniently described by identifying $\mathbb{R}^{2n}$ with $\mathbb{C}^{n}$, using coordinates $z_i = \left( x_{2i - 1}, x_{2i} \right)$, in terms of which

$$Q_n = \{ \mathbf{z} \in \mathbb{C}^n :0 \leq \left| {z_1 } \right| \leq \left| {z_2 } \right| \leq \cdots \leq \left| {z_n } \right| \leq 1 \}.$$

If these were *real* variables, we could make the volume-preserving transformation $w_1 = z_1, w_2 = z_2 – z_1, \ldots , w_i = z_i – z_{i-1}, \ldots, w_n = z_n – z_{n-1}$, with the sole restriction that the sum of the $w_i$ (all of which are nonnegative) not exceed 1. Because they are *complex* variables, though, we have to consider the area of an annulus bounded by $z_{i-1}$ and $z_i$: it is proportional to $z_i^2 – z_{i – 1}^2$. The circle of the same area has radius $w_i$ for which $w_i^2 = z_i^2 – z_{i – 1}^2$. Therefore, *if we define new variables* $w_i$ *according to this formula*, we obtain a new region- - one of substantially different shape- - having the same volume. This region is defined by $\left| {w_1 }^2 \right| + \cdots + \left| {w_n }^2 \right| \le 1$: that is, it’s our old friend $D^{2n}$. Therefore, **the volume of** $Q_n$ **equals the volume of** $D^{2n}$ .

Now for the punch line: $Q_n$ is a fundamental domain for the action of $S[n]$, the symmetric group, on the product of $n$ disks $T^{2n} = \left( D^2 \right) ^n$; $S[n]$ acts by permuting the Complex coordinates $z_1, \ldots, z_n$. The volume of $T^{2n}$ equals $|D^2|^n = \pi ^n$. Writing $|S[n]|$ for the number of permutations and equating our two completely different calculations of the volume of the $2n$ ball gives

$$\pi ^ n / |S[n]| = \frac {\Gamma \left(1/2 \right)^{2n}} { n \Gamma \left(n \right) },$$

whence

$$|S[n]| = \frac{{\pi ^n n\Gamma \left( n \right)}}{{\Gamma \left( {1/2} \right)^{2n} }}.$$

This simplifies: the volume formula for $n = 2$ must give the area of the unit circle, equal to $\pi$, whence $\Gamma \left( 1/2 \right)^2 = \pi$. Finally, then,

$$|S[n]| = n\Gamma \left( n \right).$$

I will finish by remarking that Liouville's method is a perfectly natural thing to encounter when working with the multivariate Normal distribution, so it's not really an isolated trick, but is rather a pretty basic result expressing a defining property of Normal (Gaussian) variates. There are, of course, many other ways to compute the volume of $D^{2n}$, but this one gives us the Gamma function directly.