Help finding the point(s) inside "non-star" closed shape with the highest average ray length?

Suppose the curve is star-shaped with respect to your center point $\mathbf p=(u,v)$, so that any ray emanating from $\mathbf p$ meets the curve exactly once, at say point $\mathbf q$. Then $r = \|\mathbf q - \mathbf p\|$, $\theta$ is the angle between $\mathbf q-\mathbf p$ and the $x$-axis, and the average radius is $$\frac1{2\pi}\oint_{\mathbf q\in\mathcal C}\|\mathbf q-\mathbf p\|\,\mathrm d\theta.$$ Conveniently, this integral can also be computed for non-star-shaped curves; for a ray that meets the curve multiple times, it amounts to taking the total length of all segments that lie in the interior of the curve.

Let us discretize the curve so that the integral can be computed as a sum:

curve = DiscretizeRegion[ImplicitRegion[
   x^2 + y^2 + Sin[4 x] + Sin[4 y] == 4, {{x, -3, 3}, {y, -3, 3}}]]

q = MeshCoordinates[curve];
edges = First /@ MeshCells[curve, 1];
signedAngle[a_, b_] := Arg[(Complex @@ a)/(Complex @@ b)]
avgRadius[p_] := 
 1/(2 π) Abs[
   Sum[Module[{q1, q2, r, dθ}, 
     q1 = q[[First@e]]; 
     q2 = q[[Last@e]]; 
     r = EuclideanDistance[p, (q1 + q2)/2]; (* midpoint approximation *)
     dθ = signedAngle[q1 - p, q2 - p]; 
     r dθ], 
    {e, edges}]]

(Actually I think the integral $\int r\,\mathrm d\theta$ over a line segment can be evaluated exactly in closed form, but I haven't got around to implementing that.)

Now the average radius can be evaluated pretty quickly:

avgRadius[{0, 0}]
(* 1.99725 *)

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

Then you can find the maximum using, well, FindMaximum:

FindMaximum[avgRadius[{x, y}], {{x, 0}, {y, 0}}]
(* {1.99801, {x -> -0.0525994, y -> -0.0525505}} *)

Edit: Added trapezoidal rule for higher accurary.

Rahul's code does what the OP asked for. However, it is still not too fast. The bottlenecks are signedAngle and Sum. Both can be vectorized away. The following is a function that maps a MeshRegion to a CompiledFunction that maps points to their average radius. In contrast to Rahul, we can also use the trapezoidal for integration. It comes at virtually free cost (we employ RotateRight) and should improve the accuracy significantly.

avgRadiusFunction[curve_MeshRegion] := 
  Module[{q, e1, e2}, q = MeshCoordinates[curve];
   {e1, e2} = Transpose[Developer`ToPackedArray[MeshCells[curve, 1][[All, 1]]]];
   With[{q11 = q[[e1, 1]], q12 = q[[e1, 2]], q21 = q[[e2, 1]], q22 = q[[e2, 2]]},
    Compile[{{p, _Real, 1}},
     Module[{u11, u12, u21, u22, r, dθ, r2squared, r2, p1, p2},
      p1 = Part[p, 1];
      p2 = Part[p, 2];
      u11 = q11 - p1;
      u12 = q12 - p2;
      u21 = q21 - p1;
      u22 = q22 - p2;
      r2squared = (u21^2 + u22^2);
      dθ = ArcTan[(u11 u21 + u12 u22)/r2squared, (u12 u21 - u11 u22)/r2squared];
      r2 = Sqrt[r2squared];
      Abs[(RotateRight[r2] + r2).dθ/(4. Pi)]
      ],
     RuntimeAttributes -> {Listable},
     Parallelization -> True]
    ]
   ];

Using

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]

from this post by the OP, we can create an averaged radius function like this:

f = avgRadiusFunction[curve]; // AbsoluteTiming

{0.050776, Null}

Comparison to avgRadius:

avgRadius[{0.1, 0.}] // AbsoluteTiming
f[{0.1, 0.}] // AbsoluteTiming

{0.999227, 1.99519}

{0.002347, 1.99519}

So, this is already 425 times faster.

Moreover, the resulting function f threads over list:

pp = RandomReal[{-10, 10}, {1000, 2}];
f[pp]; // AbsoluteTiming

{0.590299, Null}

Per point, this is 1800 times faster than avgRadius.

Update

In the meantime I was able to produce a version of the function above that produces CompiledFunctions that are more tolerant towards evaluation on the boundary - at the cost of some speed, though.

avgRadiusFunction[curve_MeshRegion] := Module[{q, e1, e2},
   q = MeshCoordinates[curve];
   {e1, e2} = 
    Transpose[
     Developer`ToPackedArray[MeshCells[curve, 1][[All, 1]]]];
   With[{
     q11 = q[[e1, 1]], q12 = q[[e1, 2]], q21 = q[[e2, 1]], 
     q22 = q[[e2, 2]]
     },
    Compile[{{p, _Real, 1}},
     Module[{u11, u12, u21, u22, r, dθ, r2squared, r2, p1, p2, idx, bag, z1, z2, j},
      p1 = Part[p, 1];
      p2 = Part[p, 2];
      u11 = q11 - p1;
      u12 = q12 - p2;
      u21 = q21 - p1;
      u22 = q22 - p2;
      r2squared = (u21^2 + u22^2);
      z1 = (u11 u21 + u12 u22);
      z2 = (u12 u21 - u11 u22);
      bag = Internal`Bag[{0}];
      Do[If[(r2squared[[i]] < 1. 10^-14) || ((Abs[z1[[i]]] < 1. 10^-14) && (Abs[z2[[i]]] < 1. 10^-14)), 
        Internal`StuffBag[bag, i]], {i, 1, Length[r2squared]}
       ];
      idx = Internal`BagPart[bag, All];
      If[Length[idx] > 1,
       Do[j = idx[[i]]; r2squared[[j]] = 1.; z1[[j]] = 1.; z2[[j]] = 1.;, {i, 2, Length[idx]}];
       dθ = ArcTan[z1/r2squared, z2/r2squared];
       Do[dθ[[idx[[i]]]] = 0.;, {i, 2, Length[idx]}];
       ,
       dθ = ArcTan[z1/r2squared, z2/r2squared];
       ];
      r2 = Sqrt[r2squared];
      Abs[(RotateRight[r2] + r2).dθ/(4. Pi)]
      ],
     CompilationTarget -> "WVM",
     RuntimeAttributes -> {Listable},
     Parallelization -> True
     ]]];

One approach is to parametrize the boundary, and use that parametization for the innermost integral that defines the averaged radius. This might or might not correspond to the definition you have in mind though, depending on what specifically is wanted when the region is not convex and rays from interior to boundary can go outside.

I will use the suggested example region for all computations.

expr = x^2 + y^2 + Sin[4*x] + Sin[4*y] - 4;

The idea is to set up and solve a system of differential equations using the arc-length (unit speed) parametrization. Use implicit differentiation on the defining curve to get one equation, and the unit speed condition to get another.

exprt = expr /. {x -> x[t], y -> y[t]};
odes = {D[exprt, t], x'[t]^2 + y'[t]^2 - 1}

(* Out[837]= {4 Cos[4 x[t]] Derivative[1][x][t] + 
  2 x[t] Derivative[1][x][t] + 4 Cos[4 y[t]] Derivative[1][y][t] + 
  2 y[t] Derivative[1][y][t], -1 + Derivative[1][x][t]^2 + 
  Derivative[1][y][t]^2} *)

We require a starting point so we'll pick a value for x and get one solution for y.

starty = y /. FindRoot[(expr /. x -> 1) == 0, {y, 0}]

(* Out[869]= -2.13092 *)

From that start there must be two solutions (one going in each direction). We'll take just the first and use Quiet to avoid being told there are others.

Quiet[soln = NDSolveValue[
   Join[Thread[odes == 0], {x[0] == 1, y[0] == starty}], {x[t], 
    y[t]}, {t, 0, 25}];]

{solnX, solnY} = soln;

Let's have a look.

ParametricPlot[soln, {t, 0, 20}]

To get the correct upper bound for integrating around this curve we need to find the t value where the loop closes. Probably this could have been done with WhenEvent inside NDSolveValue but I lost patience with that approach. So we'll find the first time after t=0 where we next hit the samexandyvalues as the starting points. A plot shows it is aroundt=14.5` or so.

Plot[(solnX - 1)^2 + (solnY - starty)^2, {t, 0, 20}]

The actual computation can be done with NMinimize.

{min, val} = 
 NMinimize[(solnX - 1)^2 + (solnY - starty)^2, {t, 10, 25}]

(* Out[879]= {3.19829*10^-14, {t -> 14.6198}} *)

stopt = t /. val;

--- edited part below ---

Based on rereading the original question, some responses, and comments from @HenrikSchumacher, I gather the idea is to average with respect to radian measure. So average radius is defined as below. We require explicit numeric input so that the later region integrations do not complain.

r[x_?NumberQ, y_?NumberQ] := 
  NIntegrate[Sqrt[(solnX - x)^2 + (solnY - y)^2]*
  D[ArcTan[solnX - x, solnY - y], t], {t, 0, stopt}]/(2*Pi)

Quick check: from the picture, we might expect the average radius from {1,0} to the curve to be around 2.

In[1110]:= r[1, 0]

(* Out[1110]= 1.84107 *)

This is likely to be slower than some other methods shown, but it might be slightly more accurate. Not saying that justifies the speed hit, I'm just trying to get an implementation that works correctly from basic principles, that is to say, avoiding explicit discretization..