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$
 A: The simplest way to start is to think about how the delta function is defined in 1D,
$$\int_{-\infty}^{\infty} f(x)\delta(x-a)dx = f(a)$$
which you correctly identified above.  This is saying that some function $f(x)$ multiplied by the a shifted delta function $\delta(x-a)$ and integrated everywhere returns the same function evaluated at the value given by the shift, $f(a)$.  We want the same sort of definition to apply when we integrate everywhere in 2D, so we have to have a 2-variable function $f(x,y)$ and it has to have a 2-variable delta function $\delta(x,y)$.  (Notice that $\delta(x-y)\neq \delta(x,y)$ in general.) Applying it to the 2D case in this way we should give
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = f(a,b)$$
which should be true whether you integrate either $x$ or $y$ first as shown below.
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = \int_{-\infty}^{\infty} f(a,y)\delta(y-b) \ dy  = f(a,b)$$
$$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} f(x,y)\delta(x-a,y-b) \ dx dy = \int_{-\infty}^{\infty} f(x,b)\delta(x-a) \ dx  = f(a,b)$$
These are easiest to solve when $f(x,y)$ is a separable function like $f(x,y)=g(x) h(y)$, but some harder functions are still quite doable.
A: $$\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.
A: I think something is still wrong. Because according to your answer in the integration limits $(-\infty,+\infty)$ it doesn't matter which is the function inside the $\delta$, the answer is always the same.
