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.
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.