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