How to explicitly compute the integral of the average distance between two points inside the unit disc?

Question

The average distance integral between two points $\boldsymbol x_1=(x_1,y_1)$ and $\boldsymbol x_2=(x_2,y_2)$ inside a unit disc (up to a factor) is $$\int_{-1}^1\int_{-1}^1\int_{-\sqrt{1-x_1^2}}^{\sqrt{1-x_1^2}}\int_{-\sqrt{1-x_2^2}}^{\sqrt{1-x_2^2}}\sqrt{(x_2-x_1)^2+(y_2-y_1)^2}\,\mathrm dy_2\,\mathrm dy_1\,\mathrm dx_2\,\mathrm dx_1$$ The only method I saw they solved it was using probability and solving an ODE resulting from it. Moreover, I did see how they computed explicitly a similar integral but for points on a unit square. My question is basically how can we solve this integral. The first thing coming to mind is maybe setting $x_2-x_1=r\cos\theta$ and $y_2-y_1=r\sin\theta$. Or maybe even $\boldsymbol x_1=r_1(\cos\theta_1,\sin\theta_1)$ and $\boldsymbol x_2=r_2(\cos\theta_2,\sin\theta_2)$. I've tried the first but looks like the change of coordinates isn't one to one so we need to define two other variables...