Average distance to a random point in a rectangle from an arbitrary point

In principle the integral is not too difficult using polar coordinates. I will do everything for $x_0=-x_1 ,y_0 =- y_1$ and $p=q=0$ to simplify notation a little bit. Then we don't have to shift and the whole area is given by $4$ times the integral over the first quadrant. We split the remaining integral in two regions.

1.) The triangle with vertices $\{0,0\} ,\{x_1,0\},\{x_1,y_1\}$. Here the upper bound of $r$ goes from $x_1$ to $r_1=\sqrt{x_1^2+y_1^2}$ which means that $r\in [0,\frac{x_1}{\cos(\phi)}]$. The value of $\phi$ is between $0$ and $\phi_1=\arctan(y_1/x_1)$

2.) The triangle with vertices $\{0,0\} ,\{0,y_1\},\{x_1,y_1\}$. Here the upper bound of $r$ goes from $y_1$ to $r_1=\sqrt{x_1^2+y_1^2}$ which means that $r\in [0,\frac{y_1}{\sin(\phi)}]$. The value of $\phi$ is between $\phi_1=\arctan(y_1/x_1)$ and $\pi/2$

We get: $$ \int_{-y_1}^{y_1}\int_{-x_1}^{x_1} dx dy\sqrt{x^2+y^2}=4\int_{0}^{y_1}\int_{0}^{x_1} dx dy\sqrt{x^2+y^2}=\\ 4\int_{0}^{\phi_1}d\phi\int_{0}^{x_1/\cos(\phi_1)}dr r^2+4\int_{\phi_1}^{\pi/2}d\phi\int_{0}^{y_1/\sin(\phi_1)}dr r^2=\\ \frac{4}{3}x_1^3\int_{0}^{\phi_1}d\phi \left(\frac{1}{\cos(\phi)}\right)^3+\frac{4}{3}y_1^3\int_{\phi_1}^{\pi/2}d\phi\left(\frac{1}{\sin(\phi)}\right)^3 $$

The remaining integrals can be done by IBP and the substitutions $y=\arccos{\phi}$ , $y=\arcsin{\phi}$ respectivly.

The result in the end can be simplified a bit but looks still very messy, so i spare it here. But eventually you get a closed form :)

Can you take it from here?