$\newcommand{\d}{\mathrm{d}}$I came up with this myself, but my proof will appear fairly brief, so I suspect it is not fully rigorous - although I do not see any mistakes, hence the question.
First I would like to show that:
$$\tag{1}\int_{\Omega}(g\circ f)(x)\,\d\mu(x)=\int_{f(\Omega)}g(x)\,\d(\mu\circ f^{-1})(x)$$
Where $\Omega,f(\Omega),(g\circ f)(\Omega)$ are measurable spaces, with $\mu$ the measure on $\Omega$ and $g,f$ measurable functions w.r.t the relevant sigma algebras. I think this is called "push-forward", but I have not seen a proof of it. My attempt:
Assuming $g$ and $f$ measurable, both integrals in the LHS and RHS make sense, where $f^{-1}$ is understood as preimage rather than a functional inverse, as $\mu\circ f^{-1}$ is a measure as $f^{-1}(A)$ is measurable in $\Omega$ if $A$ is measurable in $f(\Omega)$, and $f^{-1}(A\cup B)=f^{-1}(A)\cup f^{-1}(B)$ so the additivity requirements on a measure are satisfied, etc.
Use the definition of Lebesgue integral by simple functions: if $\phi\in S^+({g^+})$ defined on $(\mu\circ f^{-1})$-measurable sets $A_k$ then: $$\begin{align}\int_{f(\Omega)}g^+\,\d(\mu\circ f^{-1})&\ge\int_{f(\Omega)}\phi\,\d(\mu\circ f^{-1})\\&=\sum_{k}c_k\cdot(\mu\circ f^{-1})(A_k)\\&=\sum_kc_k\cdot\mu(E_k)\\&\le\int_{\Omega}(g\circ f)^+\,\d\mu\end{align}$$Where $E_k=f^{-1}(A_k)$. We find that a simple function in $f(\Omega)$ is a lower bound to $g^+$ iff and only if it is a lower bound to $(g\circ f)^+$ in $\Omega$, and it has the same integral in both spaces. Since the two spaces $S^+(g^+,\mu\circ f^{-1})=S^+((g\circ f)^+,\mu)$ are equal in integration, the limit-from-below definition of Lebesgue integral finds that: $$\int_{f(\Omega)}g^+\,\d(\mu\circ f^{-1})=\int_{\Omega}(g\circ f)^+\,\d\mu$$From which it follows that the integrals in $(1)$ are indeed the same.
Note - that was very, very brief. The equality of the two integrals felt obvious to me, but perhaps a true proof is more difficult than this.
I'll now move on to substitution:
Let $\Omega,f(\Omega),(g\circ f)(\Omega)$ all be measurable spaces, and suppose that $f(\Omega)$ forms a $\sigma$-finite measure space with the measure $\mu$, where $g,f$ are measurable functions. Now suppose that $f$ is bijective, with functional inverse $\psi$, and that $\psi:f(\Omega)\to\Omega$ is measurable. It follows from the measurability of $f$ that $\Omega$ is $\sigma$-finite w.r.t the measure $\mu\circ f$, but suppose also that $\Omega$ forms a sigma-finite measure space with the measure $\nu$. Then: $$\int_{f(\Omega)}g(x)\,\d\mu=\int_{f(\Omega)}(g\circ f)(\psi(x))\,\d\mu=\int_{\Omega}(g\circ f)(x)\,\d(\mu\circ \psi^{-1})=\int_{\Omega}(g\circ f)(x)\,\d(\mu\circ f)$$Where $(1)$ was used. By the $\sigma$-finiteness, we can take the Radon-Nikodym derivative of $(\mu\circ f)$ w.r.t $\nu$, and find: $$\tag{2}\int_{f(\Omega)}g(x)\,\d\mu=\int_{\Omega}(g\circ f)(x)\cdot\frac{\d(\mu\circ f)}{\d\nu}\,\d\nu$$
Let $\Omega\subset\Bbb R^n,f(\Omega)\subset\Bbb R^n$, $\mu=\mu_n$ the Lebesgue measure on $\Bbb R^n$, and $\nu=\mu_n$ also. Suppose $g$ is a measurable function, and suppose $f$ is $C^1$ in $\Bbb R^n$ and $\det J_f\neq 0$ everywhere in $\Omega$. Then the inverse function theorem gives that $\psi$ exists as a $C^1$ map everywhere in $f(\Omega)$, and is therefore measurable. As the reals are sigma finite w.r.t these measures, we can employ $(2)$, so must find the measurable function $\varphi$ such that: $$\mu_m(f(A))=\int_A\varphi(x)\,\d\mu_n(x)$$For all Lebesgue measurable $A$ in $\Bbb R^n$. Assume $A$ to not be a null set, which is fine since if it were a null set both LHS and RHS would be zero regardless of $\varphi$; indeed, take $A$ w.l.o.g (as $\varphi$ is guaranteed to be unique up to a.e. equivalence) to be a ball $B(r)$ centred about any point $x$ in $\Omega$. Then the Lebesgue differentiation theorem gives that, a.e.: $$\lim_{r\to0^+}\frac{(\mu_n\circ f)(B_r(x))}{\mu_n(B_r(x))}=\varphi(x)$$Let $\epsilon\gt0$ be fixed. As $f$ is $C^1$, we can write $f(x+h)-f(x)=J_f(x)\cdot h+o(\epsilon)$ when $\|h\|$ is small, say $0\lt\|h\|\lt\delta$. Suppose $r\lt\delta$; then, by translational invariance of Lebesgue's measure and subadditivity: $$\begin{align}\mu_n(f(B_r(x)))&=\mu_n\{f(x+h):h\in\Bbb R^n,\,\|h\|\lt r\}\\&=\mu_n\{f(x+h)-f(x):h\in\Bbb R^n,\,\|h\|\lt r\}\\&\le\mu_n\{J_f(x)\cdot h:h\in\Bbb R^n,\,\|h\|\lt r\}+\mu_n(B_{r+|o(\epsilon)|}(x)\setminus B_r(x))\end{align}$$Lebesgue's measure has the nice property that it scales by $|\det T|$ for any $T$ a linear map, so this becomes: $$\begin{align}0&\le(\mu_n\circ f)(B_r(x))\\&\le|\det J_f(x)|\mu_n(B_r(x))+\mu_n(B_{r+|o(\epsilon)|}(x)\setminus B_r(x))\\&=|\det J_f(x)|\mu_n(B_r(x))+o(\epsilon^n)\end{align}$$In division by $\mu_n(B_r(x))$, the limit becomes: $$\varphi(x)=\lim_{r\to0^+}|\det J_f(x)|+\frac{o(\epsilon^n)}{\mu_n(B_r(x))}=|\det J_f(x)|$$Therefore $|\det J_f(x)|$ is the desired R-N derivative, and using $2$ we finally have: $$\int_{f(\Omega)}g(x)\,\d\mu_n(x)=\int_{\Omega}(g\circ f)(x)|\det J_f(x)|\,\d\mu_n(x)$$
Is this a complete proof? I'm worried it's too simple...
