How come that $\int \delta(H(p,q)-E)dpdq=\Omega(E)$ not infinity? In microcanonical ensembles we have (for one particle in 1 dimension)
$$\int \delta(H(p,q)-E)dpdq=\Omega(E)$$
I am not convinced and believe that this integral diverges. Take for example a harmonic oscillator in one dimension, with $H(p,q)=p^2/2+q^2/2=E$ ($m=k=1$ for simplicity). One can write the above integral as
$$\int f(p,q)dpdq$$
The integrand is $f(p,q)=\infty$ for all points $(p,q)$ with $H(p,q)=E$. Thus, one can think of $f(p,q)$ as a distribution with infinitely many (uncountably infinite) Dirac functions sprinkled around every point $(p,q)$ on the circle $H(p,q)=p^2/2+q^2/2=E$ (see figure). Each Dirac function will contribute 1 to the integral, but there are infinitely many of them on the circle so the integral will diverge. This argument is generalizable to higher dimensions.
I need help in clearing up my confusion as to why 1) this integral does not diverge and 2) it is proportional to $\Omega(E)$, i.e., the circumference of this circle.

 A: 
The integrand is $f(p,q)=\infty$ for all points $(p,q)$ with $H(p,q)=E$.

No, it's not.

Thus, one can think of $f(p,q)$ as a distribution with infinitely many (uncountably infinite) Dirac functions sprinkled around every point $(p,q)$ on the circle

No, you can't.
The hand-waving tale that the Dirac delta $\delta(x)$ equals "zero at $x\neq 0$ and infinity at $x=0$" is just that: a hand-waving tale. It is useful for building intuition in one dimension, but that's it.
What the Dirac delta actually represents is a distribution (also sometimes known as a "generalized function"), in the formal sense of distribution theory. A distribution is a function $\varphi: \mathscr F\to \mathbb R$ which acts on the space of well-behaved functions $\mathscr F$ and assigns a real number to each function. (For example, for any function $g$ you can create a distribution that takes $f$ to $\int_{-\infty}^\infty f(x)g(x)dx$.) The Dirac delta distribution is the function
\begin{align}
\delta_0 : \mathscr F &\to \mathbb R \\
f & \mapsto \delta_0(f) = f(0)
\end{align}
which takes an arbitrary well-behaved function $f$ and returns the value of $f$ at the origin, $f(0)$.
In your integral,
$$\int \delta(H(p,q)-E)dpdq,$$
the Dirac delta is being taken over a one dimensional space, and evaluated at the variable $\epsilon = H(p,q)-E$. If your analysis does not account for that, then it is wrong.
This is easier to handle for the specific example of the harmonic oscillator, where you're trying to calculate
$$\mathrm{int}(E) = \iint \delta(\tfrac12(p^2+q^2)-E)dpdq.$$
This now has a Dirac delta evaluated over a one-dimensional energy variable which itself depends on a function of two separate integration variables $-$ which should look like quite a nightmare! To handle this, the thing to do is to make a suitable change of variables (in this case, to radial coordinates) that will separate the two. So, if we define $h=\tfrac12(p^2+q^2)$ and $\phi = \arctan(p/q)$ (so that $q=\sqrt{2h}\cos(\phi)$ and $p=\sqrt{2h}\sin(\phi)$, and $dp\,dq = dh\,d\phi$), which gives us
$$\mathrm{int}(E) = \int_0^{2\pi}\int_0^\infty \delta(h-E)dh\,d\phi= \int_0^{2\pi}1 d\phi \times \int_0^\infty \delta(h-E)dh.$$
This separation shows what's really happening: the Dirac delta is indeed singular, and it is being integrated over one dimension to give a finite result. The other integration variable then gives the circumference of the ring of equal-energy states, which is what we want to calculate with $\Omega(E)$.
