My description of the Mahalanobis distance at Bottom to top explanation of the Mahalanobis distance? includes two key results:
By definition, it does not change when the regressors are uniformly shifted.
The squared Mahalanobis distance between vectors $x$ and $y$ is given by $$D^2(x,y) = (x-y)^\prime \Sigma^{-1}(x-y)$$ where $\Sigma$ is the covariance of the data.
(1) allows us to assume the means of the regressors are all zero. It remains to compute $h_i$. However, for the claim to be true, we need to add one more assumption:
The model must include an intercept.
Allowing for this, let there be $k \ge 0$ regressors and $n$ data, writing the value of regressor $j$ for observation $i$ as $x_{ij}$. Let the column vector of these $n$ values for regressor $j$ be written $\mathbf{x}_{,j}$ and the row vector of these $k$ values for observation $i$ be written $\mathbf{x}_i$. Then the model matrix is
$$X = \pmatrix{ 1 &x_{11} &\cdots &x_{1k} \\ 1 &x_{21} &\cdots &x_{2k} \\
\vdots &\vdots &\vdots &\vdots \\ 1 &x_{n1} &\cdots &x_{nk}}$$
and, by definition, the hat matrix is
$$H = X(X^\prime X)^{-1} X^\prime,$$
whence entry $i$ along the diagonal is
$$h_i = h_{ii} = (1; \mathbf{x}_i) (X^\prime X)^{-1} (1; \mathbf{x}_i)^\prime.\tag{1}$$
There's nothing for it but to work out that central matrix inverse--but by virtue of the first key result, it's easy, especially when we write it in block-matrix form:
$$X^\prime X = n\pmatrix{1 & \mathbf{0}^\prime \\ \mathbf{0} & C}$$
where $\mathbf{0} = (0,0,\ldots,0)^\prime$ and
$$C_{jk} = \frac{1}{n} \sum_{i=1}^n x_{ij} x_{ik} = \frac{n-1}{n}\operatorname{Cov}(\mathbf{x}_j, \mathbf{x}_k) = \frac{n-1}{n}\Sigma_{jk}.$$
(I have written $\Sigma$ for the sample covariance matrix of the regressors.)
Because this is block diagonal, its inverse can be found simply by inverting the blocks:
$$(X^\prime X)^{-1} = \frac{1}{n}\pmatrix{1 & \mathbf{0}^\prime \\ \mathbf{0} & C^{-1}} = \pmatrix{\frac{1}{n} & \mathbf{0}^\prime \\ \mathbf{0} & \frac{1}{n-1}\Sigma^{-1}}.$$
From the definition $(1)$ we obtain
$$\eqalign{
h_i &= (1; \mathbf{x}_i) \pmatrix{\frac{1}{n} & \mathbf{0}^\prime \\ \mathbf{0} & \frac{1}{n-1}\Sigma^{-1}}(1; \mathbf{x}_i)^\prime \\
&=\frac{1}{n} + \frac{1}{n-1}\mathbf{x}_i \Sigma^{-1}\mathbf{x}_i^\prime \\
&=\frac{1}{n} + \frac{1}{n-1} D^2(\mathbf{x}_i, \mathbf{0}).
}$$
Solving for the squared Mahalanobis length $D_i^2 = D^2(\mathbf{x}_i, \mathbf{0})$ yields
$$D_i^2 = (n-1)\left(h_i - \frac{1}{n}\right),$$
QED.
Looking back, we may trace the additive term $1/n$ to the presence of an intercept, which introduced the column of ones into the model matrix $X$. The multiplicative term $n-1$ appeared after assuming the Mahalanobis distance would be computed using the sample covariance estimate (which divides the sums of squares and products by $n-1$) rather than the covariance matrix of the data (which divides the sum of squares and products by $n$).
The chief value of this analysis is to impart a geometric interpretation to the leverage, which measures how much a unit change in the response at observation $i$ will change the fitted value at that observation: high-leverage observations are at large Mahalanobis distances from the centroid of the regressors, exactly as a mechanically efficient lever operates at a large distance from its fulcrum.
R code to show that the relation indeed holds:
x <- mtcars
# Compute Mahalanobis distances
h <- hat(x, intercept = TRUE); names(h) <- rownames(mtcars)
M <- mahalanobis(x, colMeans(x), cov(x))
# Compute D^2 of the question
n <- nrow(x); D2 <- (n-1)*(h - 1/n)
# Compare.
all.equal(M, D2) # TRUE
print(signif(cbind(M, D2), 3))
