I want to efficiently compute weighted average of the whole cartesian plane with weight $ \overline{r}(x,y)$

The definition of $ \overline{r}(x,y)$ is the average radius of the closed region $ \mathcal{C}$ from any point.

Suppose we have curve $ \mathcal{C}$ with a ray, rotating by $ \theta$ , which emenates from point $ \mathbf{p}$ intersecting $ \mathcal{C}$ at one or more points. ( We will denote the intersections as $ \left\{\mathbf{q}_i\right\}_{i=1}^{n}$ depending on $ n$ intersections). Hence if one or more intersections exist, $ r = \sum\limits_{i=1}^{n}\|\mathbf {q}_{i} – \mathbf p\|$ , otherwise $ r=0$ .

Combining the previous definitions, the average radius function which is $ \overline{r}(\mathbf{p})=\int_{0}^{2\pi}r \ d\theta$

For star shaped curves, with $ \mathbf{p}$ inside the region, one can use an alternate and simpler definition

$ $ \overline{r}(\mathbf{p})=\frac1{2\pi}\oint_{\mathbf {q}_1\in\mathcal C}\|\mathbf {q}_1-\mathbf p\|\,\mathrm d\theta.$ $

Thanks to @Rahul, I have a code that can implement $ \overline{r}(x,y)$ . For example, if $ \mathcal{C}$ is defined by $ x^2+y^2+\sin(4x)+\sin(4y)=4$

curve = DiscretizeRegion[       ImplicitRegion[        x^2+y^2+sin(4*x)+sin(4*y) == 4, {{x, -3, 3}, {y, -4, 4}}], {{-3, 3},     {-4, 4}}, AccuracyGoal -> 8]     q = MeshCoordinates[curve];     edges = First /@ MeshCells[curve, 1];     signedAngle[a_, b_] := Arg[(Complex @@ a)/(Complex @@ b)]     avgRadius[p_] :=       1/(2 \[Pi]) Abs[Sum[Module[{q1, q2, r, d\[Theta]}, q1 = q[[First@e]];          q2 = q[[Last@e]];          r = EuclideanDistance[p, (q1 + q2)/2];(*midpoint approximation*)          d\[Theta] = signedAngle[q1 - p, q2 - p];          r d\[Theta]], {e, edges}]]    avgRadius[{x, y}] 

We can plot the contour of $ \overline{r}(x,y)$ on $ \mathcal{C}$ .

f = ContourPlot[avgRadius[{x, y}], {x, -3, 3}, {y, -3, 3},    Exclusions -> None, Contours -> 100, PlotLegends -> Automatic] 

Hence the weighted average over the entire cartesian plane is

$ $ \begin{array}{cc} {x_w=\frac{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x\overline{r}(x,y) \ dx \ dy }{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\overline{r}(x,y) \ dx \ dy}} & { y_w=\frac{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}y\overline{r}(x,y) \ dx \ dy}{\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\overline{r}(x,y) \ dx \ dy} }\end{array}$ $

I know the integral converges, since as $ \mathbf{p}$ moves further away from $ \mathcal{C}$ , the average radius approaches zero.

The problem is I am unable to accurately approximate the integral using NIntegrate under 35 minutes.

x_w=(NIntegrate[x*avgRadius[{x,y}],{x,-10,10}, {y,-10,10}])/(NIntegrate[avgRadius[{x,y}],{x,-10,10},{y,-10,10}])  y_w=(NIntegrate[y*avgRadius[{x,y}],{x,-10,10},{y,-10,10}])/(NIntegrate[avgRadius[{x,y}],{x,-10,10},{y,-10,10}]) 

How do I do accurately and efficiently approximate the weighted average of the $ x$ and $ y$ -coordinates over the entire cartesian plane with weights $ \overline{r}(x,y)$ .

“Accurately” means the error should be less than .00001. “Efficently” means computing the integral under $ 15$ minutes? I am unable to do so using NIntegrate