This is called disk line picking. For demonstration, let's use the unit circle. Essentially, to ensure a uniform distribution, we choose $x_i = \sqrt{r_i}\cos(\theta_i)$ and $y_i = \sqrt{r_i}\sin(\theta_i)$, and then we can calculate the integral $$\frac{\int_0^{2\pi}\int_0^1\int_0^{2\pi}\int_0^1 \sqrt{r_1+r_2-2\sqrt{r_1r_2}\cos(\theta_1-\theta_2)}\,\mathrm{d}r_1\,\mathrm{d}\theta_1\,\mathrm{d}r_2\,\mathrm{d}\theta_2}{\int_0^{2\pi}\int_0^1\int_0^{2\pi}\int_0^1 \mathrm{d}r_1\,\mathrm{d}\theta_1\,\mathrm{d}r_2\,\mathrm{d}\theta_2}$$
Edit: The transformation we choose is sort of explained on this page about disk point picking. Basically, if you uniformly randomly pick an $r\in (0, 1]$ and a $\theta\in [0, 2\pi)$, you'll get equal concentrations of points in each annulus, which you don't want (because an annulus further from the center will have the same number of points on average as the one closer to the center, even though the further one has a greater area). To control for this, we use the transformation described above, since annular area scales with $r^2$. We can calculate this explicitly by ensuring that $\mathrm{d}A = \mathrm{d}x\wedge \mathrm{d}y$ is constant in $r$. For $x = r\cos(\theta)$, $y = r\sin(\theta)$, we get $\mathrm{d}A = r\,\mathrm{d}r\,\mathrm{d}\theta$, but for $x = \sqrt{r}\cos(\theta)$, $y = \sqrt{r}\sin(\theta)$, we get $\mathrm{d}A = \frac{1}{2}\,\mathrm{d}r\,\mathrm{d}\theta$.