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

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.

• Re. the bounty: If you read the second paragraph of my answer carefully, you'll find that that's exactly what my answer already does. That's all I'll say about it since I'm not interested in commenting on 2-year-old posts. – user484 Dec 29 '17 at 0:26
• @Rahul Why are you not interested in 2-year old posts...... – Arbuja Dec 29 '17 at 0:42
• @Rahul For star shapes, if there is symmetry, the point with the highest average ray length is on the axis. For non-star shapes if there is no symmetry, the point with the highest average ray length is off the axis. – Arbuja Dec 29 '17 at 0:46
• @Rahul I'll stop have nice day...... – Arbuja Dec 29 '17 at 0:46
• Arbuja, since it utilizes a signed length element, I am pretty sure that Rahul's code does precisely what you aim for, even for nonstarshaped domains. If you do not think so, please give a concrete and detailed counterexample. – Henrik Schumacher Dec 31 '17 at 11:45

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

• I inputed your data but my output was First::normal : Nonatomic expression expected at position 1 in First[curve] . First::normal : Nonatomic expression expected at position 1 in First[1]  – Arbuja Dec 31 '15 at 22:07
• Sorry, I forgot to include the definition of the curve. – user484 Jan 1 '16 at 0:45
• I still got the same output as before when I inputted your data. – Arbuja Jan 1 '16 at 14:29
• How did you combine and format all parts of your calculations into one thing? I think this is reason I am not getting any output. – Arbuja Jan 1 '16 at 15:11
• I tried copy and pasting your calculation but its not gives me $RecursionLimit::reclim2 : Recursion depth of 1024 exceeded during evaluation of MeshCoordinates[curve] .$RecursionLimit::reclim2 : Recursion depth of 1024 exceeded during evaluation of MeshCoordinates[curve] . I need to find a way of doing this to apply this to many other implcit relations – Arbuja Jan 2 '16 at 22:41

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[DeveloperToPackedArray[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[
DeveloperToPackedArray[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 = InternalBag[{0}];
Do[If[(r2squared[[i]] < 1. 10^-14) || ((Abs[z1[[i]]] < 1. 10^-14) && (Abs[z2[[i]]] < 1. 10^-14)),
InternalStuffBag[bag, i]], {i, 1, Length[r2squared]}
];
idx = InternalBagPart[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
]]];

• I might just be confused today. What happens when you use the curve x^2+y^2==4? Suppose you compute the averaged radial distance for the point {2,0} (or use {1.990} if there are issues using poits on the boundary? I would expect that average to be larger than 2. Am I misunderstanding something? – Daniel Lichtblau Jan 13 '18 at 21:32
• @DanielLichtblau Hm. For a point on the boundary, at most half of the rays emanating from it actually intersect with the disk. The maximum ray length is $4$. And $4 \cdot \frac{1}{2} = 2$, so the average lenght must be less than $2$. Or not? – Henrik Schumacher Jan 13 '18 at 22:19
• If it i just inside then all rays hit the disk, and it seems the result should be on the large side of 2. (Maybe I shouldn't think too hard on weekends.) – Daniel Lichtblau Jan 13 '18 at 23:44
• @DanielLichtblau All of them hit but but most of them still have very short ray length... – Henrik Schumacher Jan 13 '18 at 23:52
• I now believe your notion of "most" is in fact what the original poster had in mind. Yours basically gets its average as radius per radian (angular) measure while mine is more radius per unit boundary length. And regardless, mine is not going to be correct for interior points where the region is not star shaped. – Daniel Lichtblau Jan 14 '18 at 15:26

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

• Daniel, I really like how you constructed the parameterization! That would never have come to my mind. – Henrik Schumacher Jan 14 '18 at 20:37
• I tried your code but it didn't work. I'm using Mathematica 11.1. – Arbuja Jan 15 '18 at 15:59
• It works for me. Did you start from a fresh kernel session? If not, my guess is some variables need to be cleared, e.g. x,y,t. – Daniel Lichtblau Jan 15 '18 at 16:52
• Nevermind, your code works. Thank You! – Arbuja Jan 17 '18 at 14:56
• The parametrization is, in effect, by theta, from 0 to 2 pi. We then consider theta as a function of t, hence that arctan derivative to get dtheta = dtheta/dt * dt. The averaging then involves dividing out by that parameter interval length. – Daniel Lichtblau Jan 21 '18 at 21:25

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]},


• So how do I find the point $(u,v)$ with the highest average radius. – Arbuja Dec 31 '15 at 20:00
• I think I confused you instead point with the highest average radius how about "average highest ray length". – Arbuja Dec 31 '15 at 23:08

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.

• When I checked using rgn = ImplicitRegion[ (x+.187762)^2 + (y+.227182)^2 + Sin[4 (x+.187762)] + Sin[4 (y+.227182)] <= 4 && -3. <= x <= 3. && -3. <= y <= 3., {x, y}]; rMean = Sqrt[Area[rgn]/Pi] I got the radius to be 1.93819. Is there some error in accuracy? – Arbuja Dec 31 '15 at 20:43
• The definition of mean radius (in your paper) is different from the area-equivalent radius the others are calculating. – george2079 Dec 31 '15 at 21:02
• Ohh I'll put this in my post – Arbuja Dec 31 '15 at 21:06