How to treat Dirac delta function of two variable? We can treat one variable delta function as 
$$\delta(f(x)) = \sum_i\frac{1}{|\frac{df}{dx}|_{x=x_i}} \delta(x-x_i).$$   
Then how do we treat two variable delta function, such as $\delta(f(x,y))$?
for example, how calculate $\int\int \delta(x-y)$ ?
I first thought using  $\int f(x)\delta(x-y)dx = f(y)$
$$\int\int \delta(x-y)dxdy = \int\int1 * \delta(x-y)dx dy = \int dy.$$ 
but this is nonsense, since we can also think it as $\int\int \delta(x-y)dxdy = \int dx$

$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\mathrm{d}y\, f(y)g(y),$$ or
$$\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\mathrm{d}x\, f(x)g(x),$$
which is the same since the integration limits are the same, and the integration variable is just a dummy variable, whose actual name doesn't matter. Now what if the limits are not the same? For example:
$$\int_{a}^{b}\mathrm{d}x\,\int_{c}^{d} \mathrm{d}y\, f(x)g(y)\delta(x,y).$$
Then we can use the Heaviside step function to write:
$$h(x)=f(x)(\theta(x-a)-\theta(x-b)),$$ and
$$l(y)=g(y)(\theta(y-c)-\theta(y-d)),$$ and so:
$$\int_{a}^{b}\mathrm{d}x\,\int_{c}^{d} \mathrm{d}y\, f(x)g(y)\delta(x,y) = \int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty} \mathrm{d}x\,\mathrm{d}y\, h(x)l(y)\delta(x,y),$$
and we can use the previous result.
Note that I have used one dimensional Cartesian variables. When dealing with non-Cartesian coordinates (for example, spherical coordinates), the Dirac delta $\delta(x,y)$ can (depending on your notation), become a scalar in the first variable $x$ and a scalar density in the second variable $y$, and you may need to insert a factor of the square root of the determinant of the metric to make everything work. If you need more details let me know - coincidentally a couple of people had asked me about this recently, so I have a lot of it LaTeXed, just not handy at the moment.
