Double integration with Indicator function The integral of interest is:
$$Q = \int_{\mathbb{R}^2}\int_{\mathbb{R}^2} I\left(\frac{1}{2}\frac{x_2^2 - x_1^2 + y_2^2 - y_1^2}{x_2-x_1} \in [0,1]\right) \nonumber \\
\times I\left(2\arcsin\left(\frac{\sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}}{2\sqrt{y_1^2 + \left(x_1 - \frac{1}{2}\frac{x_2^2 - x_1^2 + y_2^2 - y_1^2}{x_2-x_1}\right)^2}}\right) > \theta\right) \\ \times \exp\left(-C\left(y_1^2 + \left(x_1 - \frac{1}{2}\frac{x_2^2 - x_1^2 + y_2^2 - y_1^2}{x_2-x_1}\right)^2\right)\right)\mathrm{d}Z_1\mathrm{d}Z_2,$$
where $C > 0$, $Z_1 = (x_1, y_1)$, $Z_2 = (x_2,y_2)$, $I(\cdot)$ is an indicator function, and $\theta \in (0, \pi]$.

Start with the stretch-rotation transformations $$\begin{pmatrix} z_1 \\ z_2\end{pmatrix} = \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2\end{pmatrix} \qquad \text{and} \qquad \begin{pmatrix} z_3 \\ z_4\end{pmatrix} = \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix} \begin{pmatrix} y_1 \\ y_2\end{pmatrix}$$
so that ${\rm d}x_1 {\rm d}x_2 {\rm d}y_1 {\rm d}y_2=\frac{1}{4} {\rm d}z_1 {\rm d}z_2 {\rm d}z_3 {\rm d}z_4$ and the integral becomes $$4Q = \int_{\mathbb{R}^4} I\left(\frac{1}{2}\left(z_1+\frac{z_3z_4}{z_2}\right) \in [0,1]\right) \times I\left(2\arcsin\left(\sqrt{\frac{{z_2^2 + z_4^2}}{{z_2^2+z_4^2 + z_3^2\left(1 + \frac{z_4^2}{z_2^2}\right)}}}\right) > \theta\right) \times \exp\left[-\frac{C}{4}\left( z_2^2 + z_4^2 + z_3^2 \left(1+\frac{z_4^2}{z_2^2} \right) \right)\right] \, \mathrm{d}z_1\mathrm{d}z_2\mathrm{d}z_3\mathrm{d}z_4$$ which eliminates the first indicator function by carrying out the integral over $z_1$ and gives a factor of 2. Furthermore the second indicator function proposes the transformation to polar coordinates i.e. $(z_2,z_4)=(r\cos\phi,r\sin\phi)$ and we have $$2Q=\int_{\mathbb{R}}{\rm d}z_3\int_0^\infty r\,{\rm d}r \int_0^{2\pi}{\rm d}\phi \, I\left(2\arcsin\left(\frac{r}{\sqrt{r^2+\frac{z_3^2}{\cos^2\phi}}}\right)>\theta\right) \exp\left[ -\frac{C}{4} \left( r^2 + \frac{z_3^2}{\cos^2\phi} \right) \right]  .$$ Next we carry out the $z_3$ integration which is 2 times the integral from $0$ to $\infty$. For $z_3=0$ the argument of the arcsine is $1$ and so $\pi>\theta$ is definitely true. The other bound follows from solving $$\frac{r}{\sqrt{r^2+\frac{z_3^2}{\cos^2\phi}}}=\sin\frac{\theta}{2} \qquad \Longrightarrow \qquad z_3=r|\cos \phi| \cot \frac{\theta}{2} \, .$$ Hence
$$Q= \sqrt{\frac{\pi}{C}} \int_0^\infty r \, {\rm d}r \int_0^{2\pi} {\rm d}\phi \, |\cos\phi| \, \exp\left[-\frac{Cr^2}{4}\right] \, {\rm erf}\left( \frac{r\sqrt{C} \, \cot \frac{\theta}{2}}{2} \right) \\ = 4 \sqrt{\frac{\pi}{C}} \int_0^\infty {\rm d}r \, r \, \exp\left[-\frac{Cr^2}{4}\right] \, {\rm erf}\left( \frac{r\sqrt{C} \, \cot \frac{\theta}{2}}{2} \right) \\ 
= \frac{8\cot \frac{\theta}{2}}{C} \int_0^\infty {\rm d}r \, \exp\left[-\frac{Cr^2}{4\sin^2 \frac{\theta}{2}} \right] \\
=\frac{8 \sqrt{\pi} \cos \frac{\theta}{2}}{C^{3/2}} \, .$$
