Judging from the setup of the problem, I'll assume $x_1, t > 0$ and only study that case.
Let $\lambda = \frac12x_1$ and $\mu = \frac12 ct$. Introduce a bunch of variables:
$$
\begin{cases}
(x,y) = (z_1 - \lambda, z_2)\\
r_1 = \sqrt{(x+\lambda)^2 + y^2} = \sqrt{z_1^2+z_2^2}\\
r_2 = \sqrt{(x-\lambda)^2 + y^2} = \sqrt{(z_1-x_1)^2 + z_2^2}\\
(u,v) = \left(\frac{r_1 - r_2}{2},\frac{r_1 + r_2}{2}\right)
\end{cases}$$
The actual integral you want to integral can be rewritten as
$$\mathcal{I} = \int_{[-\lambda,\infty)\times(-\infty,\infty)}\rho(x+\lambda)\delta\left(\frac2c(u+\mu)\right)\frac{dxdy}{r_1r_2}$$
Since the integrand is symmetric with respect to the $y$-axis, we can replace
the integral by one over the upper $(x,y)$ plane:
$$\int_{[-\lambda,\infty)\times(-\infty,\infty)}\cdots\quad\leftrightarrow\quad 2\int_{[-\lambda,\infty)\times [0,\infty)}\cdots$$
Change variables to $(u,v)$, the upper $(x,y)$ plane corresponds to the
strip $[-\lambda,\lambda] \times [\lambda,\infty)$ in the $(u,v)$ plane.
Notice
$$\begin{align}
4uv &= (r_1-r_2)(r_1+r_2) = r_1^2 - r_2^2 = 4\lambda x\\
u^2 + v^2 &= \frac14\left((r_1-r_2)^2+(r_1+r_2)^2\right)
= \frac12\left(r_1^2+r_2^2\right)
= x^2+y^2+\lambda^2
\end{align}
$$
A point $(x,y) \in [-\lambda,\infty) \times [0,\infty)$ is equivalent to corresponding $(u,v)$ belongs to following region:
$$\Omega \stackrel{def}{=} \{ (u,v) \in [ -\lambda,\lambda ] \times [ \lambda, \infty ) : uv \ge -\lambda^2 \}$$
In terms of $(u,v)$, the area element has the form
$$dx \wedge dy = dx \wedge \frac{dy^2}{2y} =
\frac{1}{2\lambda y} d(\lambda x)\wedge d(x^2+y^2+\lambda^2)
= \frac{1}{2\lambda y} d(uv)\wedge d(u^2+v^2)\\
= \frac{1}{\lambda y} (udv + vdu)\wedge(udu + vdv)
= \frac{v^2-u^2}{\lambda y} du\wedge dv
= \frac{r_1r_2}{\lambda y} du\wedge dv
$$
Since
$$y^2 = u^2 + v^2 - x^2 - \lambda^2 = u^2 + v^2 - \frac{u^2v^2}{\lambda^2} - \lambda^2
= \frac{1}{\lambda^2}(v^2 - \lambda^2)(\lambda^2-u^2)$$
We find
$$\frac{dxdy}{r_1r_2} = \frac{dudv}{\lambda y} =
\frac{dudv}{\sqrt{(v^2-\lambda^2)(\lambda^2-u^2)}}
$$
and hence
$$\mathcal{I} = 2\int_\Omega \rho(x+\lambda)\frac{\delta(\frac{2}{c}(u+\mu))dudv}{\sqrt{(v^2-\lambda)^2(\lambda^2-u^2)}}$$
Apply the formula mentioned in other answers, we have
$$\delta\left(\frac2c(u+\mu)\right) = \frac{c}{2}\delta(u + \mu)$$
Integrate the delta function over $u$, the $u$ get fixed to $-\mu$.
When $\mu > \lambda$, the line $u = -\mu$ lies outside the strip
$[-\lambda,\lambda] \times [\lambda,\infty)$ and $\mathcal{I}$ vanishes.
For $0 < \mu < \lambda$, the integral reduces to
$$\mathcal{I} = \frac{c}{\sqrt{\lambda^2-\mu^2}} \int_\lambda^{\lambda^2/\mu}\frac{\rho\left(\lambda - \frac{\mu v}{\lambda}\right) d v}{\sqrt{v^2-\lambda^2}}
= \frac{c}{\sqrt{\lambda^2-\mu^2}}\int_0^{\lambda-\mu} \frac{\rho(z_1)dz_1}{\sqrt{(\lambda-z_1)^2-\mu^2}}
$$
The rest of the integral depends on the actual form of $\rho(\cdot)$ and I'll leave that for you.
