Just need confirmation that the following theorem is correct.
Let $X$ be a continuously distributed $n$-variate random variable with density $\varrho_X$ and $f : \mathrm{R}^n \rightarrow \mathrm{R}^n$ a continuously differentiable function. For the density $\varrho_{f(X)}$ of $f(X)$ it holds \begin{align*} \varrho_{f(X)}(y) = \int_{f^{-1}(y)} \frac{\varrho_X(x)}{\left| \det D f(x) \right|} \,V(\mathrm{d} x). \end{align*}
Edit: Was definitely not correct. Maybe
$$\varrho_{ f(\boldsymbol X)}(\boldsymbol y) = \int_{\mathbb R^n} \delta(\boldsymbol y - f(\boldsymbol x)) \varrho_{\boldsymbol X}(\boldsymbol x) d\boldsymbol x = \sum_{\boldsymbol x : f(\boldsymbol x)= \boldsymbol y} \frac{\varrho_{\boldsymbol X}(\boldsymbol x)}{ \left| \det D f(\boldsymbol x) \right|} \,\boldsymbol{1}_{\det D f(\boldsymbol x) \neq 0}(\boldsymbol x)\ ?$$ This reduces correctly to the univariate and diffeomorphic case.
What you have is a formula for $f: \mathbb R^n \to \mathbb R$. If $\nabla f \neq 0$ on the surface $f(\boldsymbol x) = y$, then $$\varrho_{f(\boldsymbol X)}(y) = \int_{\mathbb R^n} \delta(y - f(\boldsymbol x)) \varrho_{\boldsymbol X}(\boldsymbol x) d\boldsymbol x = \int_{f(\boldsymbol x) = y} \frac {\varrho_{\boldsymbol X}(\boldsymbol x)} {| \nabla f(\boldsymbol x)|} dS(\boldsymbol x).$$
If $\boldsymbol f$ is $\mathbb R^n \to \mathbb R^n$, then $\delta(\boldsymbol y - \boldsymbol f(\boldsymbol x))$ is a product of $n$ factors $\delta(y_i - f_i(\boldsymbol x))$. Assuming $\boldsymbol f$ has a smooth inverse in a neighborhood of $\boldsymbol x$ if $\boldsymbol f(\boldsymbol x) = \boldsymbol y$, $$\varrho_{\boldsymbol f(\boldsymbol X)}(\boldsymbol y) = \int_{\mathbb R^n} \delta(\boldsymbol y - \boldsymbol f(\boldsymbol x)) \varrho_{\boldsymbol X}(\boldsymbol x) d\boldsymbol x = \sum_{\boldsymbol x: \boldsymbol f(\boldsymbol x) = \boldsymbol y} \frac {\varrho_{\boldsymbol X}(\boldsymbol x)} {\left| \det D \boldsymbol f(\boldsymbol x) \right|}.$$
These can be viewed as special cases of the same general formula $$\varrho_{\boldsymbol f(\boldsymbol X)}(\boldsymbol y) = \int_{\mathbb R^n} \delta(\boldsymbol y - \boldsymbol f(\boldsymbol x)) \varrho_{\boldsymbol X}(\boldsymbol x) d\boldsymbol x = \int_{\boldsymbol f(\boldsymbol x) = \boldsymbol y} \frac {\varrho_{\boldsymbol X}(\boldsymbol x)} {\sqrt {\det (J J^t)}} dS(\boldsymbol x),$$ where $J = D \boldsymbol f(\boldsymbol x)$ is the $k \times n$ Jacobian matrix of $\boldsymbol f: \mathbb R^n \to \mathbb R^k$ and the integral is over an $(n - k)$-dimensional surface.
For $k > n$, we get a degenerate distribution: the pdf contains $k - n$ delta functions that haven't been integrated out.