Compute area of a sphere through a Dirac delta I've been having issues with integrating with a Dirac delta. To compute the area of a sphere centered at $(0,0,0)$ it seems to work just fine:
$$\int_{0}^{2\pi}{\int_{0}^{\pi}{\int_{0}^{\infty}{\delta(r-\rho)r^2\sin\theta\, dr}\,d\theta}\,d\phi} = 4\pi\rho^2$$
Now I will take the same sphere but offset by $(0,0,\rho)$, that is: $x^2 + y^2 + (z-\rho)^2 = \rho^2$. Going to spherical coordinates yields: $r^2\cos^2\phi\sin^2\theta + r^2\sin^2\phi\sin^2\theta + (r\cos\theta-\rho)^2 = \rho^2$, which yields: $r(r-2\rho\cos\theta)=0$, and we can express the sphere in spherical coordinates as: $r(\theta) = 2\rho\cos\theta, \theta \in [0,\pi/2], \phi\in[0,2\pi]$.
Integrating yields:
$$\int_{0}^{2\pi}{\int_{0}^{\frac{\pi}{2}}{\int_{0}^{\infty}{\delta(r-2\rho\cos\theta)r^2\sin\theta\, dr}\,d\theta}\,d\phi} = \frac{8\pi\rho^2}{3}$$
Now this is not right clearly. The only reason I can think of has to do something with properties of the Dirac  delta I am unaware of. Note that I have not studied measure theory. I need the Dirac delta and not a surface integral  because I will be using this to compute transformations of probability density functions which I will need to write through a Dirac delta.
Edit: References covering this for engineers/computer science students are welcome.
Edit 2:
Taking into account David Holden's answer I came up with the following fact which must hold (I hope it's correct):
$$\int_{V}{\delta(f(x)) \,dx} = \int_{S = \{x|f(x) =0\}}{\,dA}$$
Edit 3:
I found some more information on the subject:
Impulse functions over curves and surfaces
Properties_in_n_dimensions
Surface area from indicator function
Property of Dirac delta function in $\mathbb{R}^n$
Does the coarea formula hold for delta-function?
I believe the issue was that whenever I offset the sphere the Dirac delta changed such that $\delta(f(r)) \rightarrow \delta(g(r,\theta))$ and $g$ was then a non-trivial mapping (so it's not the one-dimensional dirac delta I am used to anymore). Based on the first article I believe that I can rewrite it as a surface Dirac delta $\delta(g(r,\theta)) = \delta_S(r,\theta)$ which yields the surface integral giving a correct result. The other threads and wikipedia state that I should have a normalization by the magnitude of the gradient. I think I am missing an important piece since for the result to be correct this normalization factor should cancel out with something. More precisely:
$$\int_V{\delta(r-2\rho\cos\theta)r^2\sin\theta \,dr\,d\theta\,d\phi} = \int_S{\frac{\,d\sigma}{\sqrt{r^2+\rho^2-2r\rho\cos\theta}}}$$
The only idea I have is that somehow the normalization factor will pop out of the $\,d\sigma$. No idea though since it's supposed to be a 'Minkowski content measure' which is way over my head as a computer science student.
To add to this I would also like to be able to solve the same problem with a heaviside function (for integrating the volume of the offset ball). I am unsure whether similar considerations apply there, however if I integrate it, the result seems correct. I still want to make sure this is valid for other volumes also (maybe it's just a coincidence like the sphere centered at $(0,0,0)$). So I would be grateful if somebody with more knowledge on geometric measure theory could clarify all of the points.
 A: To summarize the discussion in the comments, the definition of $\delta(f)$ is derived from postulating two basic properties: the substitution rule
$$\int_{\mathbb R^n} \delta(f(\mathbf x)) \,\phi(\mathbf x) \,d\mathbf x =
\int_U \delta(f(\mathbf x(\mathbf u))) \,\phi(\mathbf x(\mathbf u))
 \left| \det D \mathbf x(\mathbf u) \right| d\mathbf u$$
and
$$\int_{\mathbb R^n} \delta(x_1) \,\phi(\mathbf x) \,d\mathbf x =
\int_{\mathbb R^{n - 1}} \phi(\mathbf x) \rvert_{x_1 = 0} \,dx_2 \cdots dx_n.$$
If you try to set
$$\small \int \delta(f(\mathbf x)) d\mathbf x =
\int_{f(\mathbf x) = 0} dS =
\int_{2 f(\mathbf x) = 0} dS =
\int \delta(2 f(\mathbf x)) d\mathbf x,$$
you violate the first rule. If you try to set
$$\small \iint \delta(r - f(\theta)) \phi(r, \theta) dr d\theta \neq
\int \phi(f(\theta), \theta) d\theta,$$
you violate the second rule. If you adopt the definition that the rest of the world is using derived from the two stated properties, you get the identity
$$\int_{\mathbb R^n} \delta(f(\mathbf x))
 \left| \nabla f(\mathbf x) \right| \phi(\mathbf x) \,d\mathbf x =
\int_{f(\mathbf x) = 0} \phi(\mathbf x) \,dS(\mathbf x),$$
which is formally the same as the coarea formula because both are essentially the same change of variables formula. The first two formulas in your question will in fact be correct, while the last two will be incorrect.
A: when you write integrating yields... you may be making an unwarranted assumption about what is the "element of area". a simple geometric approach suggests:
$$
dA = \rho d(2\theta) \rho \sin 2\theta d\phi = 2 \rho^2 \sin 2\theta d\theta d\phi
$$
note that with $dA$ thus defined:
$$
\int_0^{2\pi}\int_0^{\frac{\pi}2} dA = 4\pi\rho^2 
$$
A: I finally figured out why I am getting a 'wrong' result. As expected I cannot substitute with the delta directly since it's a composition with a submersion. However the following equality holds from the coarea formula:
$$\int_{R^n}{f(x)\delta(g(x))\,dx} = \int_{g^{-1}(0)}{\frac{f(x)}{|\nabla g(x)|}\,d\sigma(x)}$$
Where $g:R^n\rightarrow R$, $|\nabla g(x)|\ne 0$, and $d\sigma$ is the measure on the surface $g^{-1}(0)$.
Let us consider the non-normalized uniform probability density function on the sphere with center $(0,0,0)$ and radius $\rho$ in spherical coordinates: $p_A(x,y,z) = \delta(r-\rho)r^2\sin\theta$. Unsurprisingly integrating it yields $4\pi\rho^2$:
$$\int_{0}^{2\pi}{\int_{0}^{\pi}{\int_{0}^{\infty}{\delta(r-\rho)r^2\sin\theta\,dr}\,d\theta}\,d\phi} = \int_{0}^{2\pi}{\int_{0}^{\pi}{\frac{\rho^2\sin\theta}{1}\,d\theta}\,d\phi} = 4\pi\rho^2$$
Note that the division by $1$ is to emphasize that $|\nabla g| = 1$. That is, I have used the coarea formula above even if it may seem unnecessary (but as we'll see later it is actually important for other mappings, and this is simply a special case where we have the standard delta function). Now let us compute the area of the translated sphere. Going to Cartesian coordinates gives us: $p_B(x,y,z) = \frac{p_A(r,\theta)}{|r^2\sin\theta|} = \delta(\sqrt{x^2+y^2+z^2}-\rho)$ (we have used the invertible pdf transformation theorem). Translating by $(0,0,\rho)$ yields: $p_C(x,y,z) = p_B(x,y,z-\rho)$, where the Jacobian of this transformation is $1$. Finally going back to spherical coordinates we have: 
$$p_D(r,\theta) = \delta(\sqrt{r^2\cos^2\phi\sin^2\theta + r^2\sin^2\phi\sin^2\theta + r^2\cos\theta^2 + \rho^2 - 2r\rho\cos\theta}-\rho)r^2\sin\theta = \\ 
= \delta(\sqrt{r^2+\rho^2-2r\rho\cos\theta}-\rho)r^2\sin\theta$$
We compute the gradient of the mapping as: $\nabla g(r,\theta) = \frac{1}{2\sqrt{r^2+\rho^2-2\rho\cos\theta}}(2r-2\rho\cos\theta, 2\frac{r}{r}\sin\theta,0)$. Finally $|\nabla g(r,\theta)| = 1$. We may compute $g^{-1}(0) = \{(2\rho\cos\theta, \theta, \phi),\theta \in [0,\frac{\pi}{2}], \phi \in [0,2\pi]\}$. The surface area element is $dA = 2\rho^2\sin2\theta\,d\theta\,d\phi$. Then finally:
$$\int_{0}^{2\pi}{\int_{0}^{\frac{\pi}{2}}{\int_{0}^{\infty}{p_D(r,\theta)\,dr}\,d\theta}\,d\phi}
= \int_{0}^{2\pi}{\int_{0}^{\frac{\pi}{2}}{\rho^2\sin2\theta\,d2\theta}\,d\phi} = 4\pi\rho^2$$
Now let us consider a slightly different variant: $p_A(r) = \delta(r^2 - \rho^2)r^2\sin\theta$, $|\nabla g(r)| = 2r$
$$\int_{0}^{2\pi}{\int_{0}^{\pi}{\int_{0}^{\infty}{\delta(r^2-\rho^2)r^2\sin\theta\,dr}\,d\theta}\,d\phi} = \int_{0}^{2\pi}{\int_{0}^{\pi}{\frac{\rho^2\sin\theta}{2\rho}\,d\theta}\,d\phi} = 2\pi\rho$$
Rather surprisingly (at least for me) we get a different result, which however for the $\delta$ defined as is, is supposedly correct (I believe that the result being $2\pi\rho$ is just a lucky coincidence). So one has to be careful about the mapping.
After transforming to cartesian coordinates, translating and returning to spherical coordinates we get $p_D(r, \theta) = \delta(r^2-2\rho\cos\theta)r^2\sin\theta$, $|\nabla g(r,\theta)| = 2\sqrt{r^2+\rho^2-2r\rho\cos\theta}$. Using the coarea formula once again:
$$\int_{0}^{2\pi}{\int_{0}^{\frac{\pi}{2}}{\int_{0}^{\infty}{p_D(r,\theta)\,dr}\,d\theta}\,d\phi}
= \int_{0}^{2\pi}{\int_{0}^{\frac{\pi}{2}}{\frac{\rho^2\sin2\theta}{2\sqrt{\rho^2}}\,d2\theta}\,d\phi} = 2\pi\rho$$
In conclusion, it seems that it is not correct to substitute directly when the delta function is composed with a function different than the identity (or $\pm const$). In that specific case the coarea formula has to be used. Additionally we seem to have the relationship $\delta_S(x) = \delta(g(x))|\nabla g(x)|$, where $S=g^{-1}(0)$:
$$\int_{R^n}{f(x)\delta(g(x))|\nabla g(x)|\,dx} = \int_{R^n}{f(x)\delta_S(x)\,dx} = \int_{S}{f(x)\,d\sigma(x)}$$
I very much appreciate the input from Maxim and David Holden, which ultimately helped me figure this out.
Edit: A very interesting read I found later: https://www.mathpages.com/home/kmath663/kmath663.htm
It certainly helps to understand the problem from an intuitive point of view also.
