Is there a simple proof for the (co)area formula in the case of equal dimensions? Let $\Omega \subseteq \mathbb{R}^n$ be a nice domain (say a ball), and let
$f:\Omega \to \mathbb{R}^n$ be smooth. Then
$$ \int_{\Omega} |\det df|=\int_{\mathbb{R}^n} |f^{-1}(y)|dy, \tag{1}$$
where $|f^{-1}(y)|$ is the cardinality of the set $\{f^{-1}(y)\}$.
Formula $(1)$ is a special case of the area and co-area formulas, when the dimensions of the source and target are equal. 
In general, the area\co-area formulas:


*

*Deal with domains of different dimensions.

*Allow the domain to be any measurable set.

*Allow $f$ to be less regular (require only Lipschitzity).
The proofs for the general claims are not very short or easy.

Question: Is there a simpler proof for the specific case mentioned above? (Equal dimensions, $f$ smooth, nice smooth domain).

A reference would suffice (I couldn't find one).
 A: I think I have it:
First assume $f$ is orientation-preserving, i.e. $df$ is invertible and has positive determinant everywhere.
Then by degree theory $$ \int_{\Omega} |\det df|=\int_{\Omega} \det df=\int_{\mathbb{R}^n} | f^{-1}(y)|. \tag{1}$$
If $f$ is orientation-reversing, a similar computation works.
Now, for the general case: By Sard's theorem, for almost every $y \in \mathbb{R}^n$, $y$ is a regular value of $f$. Let $Z \subseteq \mathbb{R}^n$ be the set of non-regular values.
Now partition $\Omega=\Omega_0 \cup \Omega_+ \cup \Omega_-,$ where


*

*$ \Omega_0=\{ x \in \Omega \, | \, \det df_x =0 \}$

*$ \Omega_+=\{ x \in \Omega \, | \, \det df_x >0 \}$

*$ \Omega_-=\{ x \in \Omega \, | \, \det df_x <0 \}$


Note that $\Omega_+,\Omega_-$ are open sets. Then, using the above mentioned result for the orientation preserving (reversing) maps $f_{\Omega_+}$  ($f_{\Omega_-}$), we obtain
$$ \int_{\Omega} |\det df|=\int_{\Omega_+} |\det df|+\int_{\Omega_-} |\det df|=\int_{\mathbb{R}^n} | f_{\Omega_+}^{-1}(y)|+\int_{\mathbb{R}^n} | f_{\Omega_-}^{-1}(y)|=$$
$$ \int_{\mathbb{R}^n} | f_{\Omega_+}^{-1}(y)|+|f_{\Omega_-}^{-1}(y)|=\int_{\mathbb{R}^n \setminus Z}| f_{\Omega_+}^{-1}(y)|+|f_{\Omega_-}^{-1}(y)|=\int_{\mathbb{R}^n \setminus Z} | f^{-1}(y)|$$
where we used the fact $Z$ is of measure zero, and that on $\mathbb{R}^n \setminus Z$, every pre-image belongs to $\Omega_+$ or $\Omega_-$.
