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

Question

I am trying to find the point(s) with the highest average ray length at any given point in a closed curve defined as an implicit function, (ex. $x^2+y^2+\sin(4x)+\sin(4y)=4$). I have some idea of mathematically approaching this problem but have little experience with Mathematica. Their documentation section does little to help.

Mathematically, I'd convert the implicit relation to polar coordinates. With a circle I can convert the polar coordinates explicitly into $r=g(\theta)$ but for most implicit shapes this is not possible.

However, assuming we can convert $f(\theta)$ to polar coordinates $f(r\cos(\theta)+u,r\sin(\theta)+v)$ and solve for $r$, the average radius is $\int_{0}^{2\pi} r d\theta$

You must also check this link to understand "my definition" of highest average ray length for a circle https://files.acrobat.com/a/preview/38a501f2-5b63-4b76-b720-6cadb9c3e142. Note I was able to take the area of the polar equation explcitly in this case.

Answer 1

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}} *)

Answer 2

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
     ]]];

Answer 3

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..

Answer 4

EDIT 2: Updated status on Wolfram Technical Support CASE:3501039

EDIT :

rgn1 = ImplicitRegion[
   x^2 + y^2 + Sin[4 x] + Sin[4 y] <= 4. && -2.4 <= x <= 2 && -2.4 <= y <= 
     2, {x, y}];

The radius of the circle with the same area as rgn1 is

rMean1 = Sqrt[Area[rgn1]/Pi]

(*  2.01147  *)

For reasons that are unclear to me, the area of the ImplicitRegion decreases if the bounds on x and y are widened! This appears to be a bug to me (Wolfram Technical Support CASE:3501039 Reply from Wolfram Technical Support: "I was able to reproduce the issue, and consequently I filed a report with our development team raising the issues").

rgn2 = ImplicitRegion[
   x^2 + y^2 + Sin[4 x] + Sin[4 y] <= 4. && -3. <= x <= 3. && -3. <= y <= 
     3., {x, y}];

The radius of the circle with the same area as rgn2 is

rMean2 = Sqrt[Area[rgn2]/Pi]

(*  1.94023  *)

Visually, the larger radius (rMean1) is a better fit.

Show[
 RegionPlot[rgn1],
 ContourPlot[x^2 + y^2,
  {x, -rMean1, rMean1}, {y, -rMean1, rMean1},
  Contours -> {rMean1^2, rMean2^2},
  ContourStyle -> {Directive[Red, Dashed],
    Directive[Red, Dashed, Thick]},
  ContourShading -> None]]

Answer 5

Note that you have defined mean radius as 1/(2 Pi) Integrate[ r[t] ,{t,0,2Pi}] rather than (then more usual) equivalent area definition : Sqrt[1/(2 Pi) Integrate[ r[t]^2 ,{t,0,2Pi}]]

an approximate way to get your center point:

ContourPlot[ 
 x^2 + y^2 + Sin[4 x] + Sin[4 y] == 4, {x, -3, 3}, {y, -3, 3}]
curve = First@Cases[Normal@%, Line[l_List] :> l , Infinity]

p0 = {x, y} /. 
  Last@NMinimize[Total[Norm[{x, y} - # ] & /@ curve], {x, y}]

{-0.187762, -0.227182}

rmean = Mean[ Norm[p0 - # ] & /@ curve ]

2.00439

Graphics[{Line[curve], {Dashed, Circle[p0, rmean] }, 
  Arrow[ {p0, p0 + rmean {Sqrt[3], 1}/2}]}]

now we can get a better approximation to the mean radius (as defined in link):

sa[a_] := Module[{aa = a}, While[aa <= 0, aa += 2 Pi]; aa]
rbar[{x0_?NumericQ, y0_?NumericQ}] := Module[{},
  a0 = ArcTan[#[[1]], #[[2]]] &@({x0, y0} - curve[[1]]);
  ii = Interpolation[({sa[ArcTan[#[[1]], #[[2]]] - a0],
          Norm[#]} & @({x0, y0} - #)) & /@
      curve[[;; -2]]];
  NIntegrate[ (ii)[theta] , {theta, ii["Domain"][[1, 1]],
     ii["Domain"][[1, 2]]}]/(-Subtract @@ ii["Domain"][[1]]) ]
rbar[p0]

1.99073

I planned to apply FindMinimuim to rbar but it didn't readily work.