Convexity Coefficient

Author

Eric W. Weisstein
September 21, 2010

This notebook downloaded from http://mathworld.wolfram.com/notebooks/ComputationalGeometry/ConvexityCoefficient.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/ConvexityCoefficient.html.

©2010 Wolfram Research, Inc. except for portions noted otherwise

Monte Carlo

Setup

Graphics`Mesh`MeshInit[]

Internal`SetVisualizationOptions[EnableMeshObject->True]

?Graphics`Mesh`*Q

?Graphics`Mesh`*

France

Graphics[p=Polygon[CountryData["France","Polygon"][[1,1]]]]

Length[p[[1]]]

mesh=Graphics`Mesh`MeshObject[Graphics[p]]

{{xl,xh},{yl,yh}}=mesh@"Range"

makept[]:={RandomReal[{xl,xh}],RandomReal[{yl,yh}]}
makeseg[]:=Line[{makept[],makept[]}]

nseg=50;
Graphics[{Opacity[0.5],p,Opacity[1],Table[(seg=makeseg[];{If[Graphics`Mesh`IntersectQ[Graphics[{p,seg}],AllPoints->False],Red,Yellow],seg}),{i,nseg}]}]

nseg=100;
allseg=Table[makeseg[],{i,nseg}];

pinter=Graphics`Mesh`FindIntersections[Graphics[Flatten[{p,allseg},1]]];

Graphics[{Opacity[0.5],p,Opacity[1],Brown,allseg,Red,Point[pinter]}]

pinter=Flatten[Graphics`Mesh`FindIntersections[Graphics[{p,#}],AllPoints->False]& /@ allseg,1];

Graphics[{Opacity[0.5],p,Opacity[1],Brown,allseg,Red,Point[pinter]}]

Pentagram

poly=Polygon[{{(1/2),(1/2) (Sqrt[1-(2/(Sqrt[5]))])},{(1/2) (-2+(Sqrt[5])),(1/2) (Sqrt[1-(2/(Sqrt[5]))])},{0,Root[1-5 ((#1)^2)+5 ((#1)^4)&,3]},{1-((Sqrt[5])/2),(1/2) (Sqrt[1-(2/(Sqrt[5]))])},{-(1/2),(1/2) (Sqrt[1-(2/(Sqrt[5]))])},{(1/4) (-3+(Sqrt[5])),-(1/(Sqrt[130+58 (Sqrt[5])]))},{(1/4) (1-(Sqrt[5])),Root[1-20 ((#1)^2)+80 ((#1)^4)&,1]},{0,Root[1-25 ((#1)^2)+5 ((#1)^4)&,2]},{(1/4) (-1+(Sqrt[5])),Root[1-20 ((#1)^2)+80 ((#1)^4)&,1]},{(1/4) (3-(Sqrt[5])),-(1/(Sqrt[130+58 (Sqrt[5])]))}}]//N;

mesh=Graphics`Mesh`MeshObject[Graphics[poly]]

InputForm[%]

{{x1,x2},{y1,y2}}=mesh@"Range"

ineq[x_,y_]:=(2+(5+3 (Sqrt[5])) x≥(Sqrt[2 (5+(Sqrt[5]))]) y&&((2+(Sqrt[50+22 (Sqrt[5])]) y≥(5+(Sqrt[5])) x&&(Sqrt[5])+2 (Sqrt[5]) x+2 (Sqrt[5+2 (Sqrt[5])]) y≤1)||((Sqrt[5])+2 (Sqrt[5]) x+2 (Sqrt[5+2 (Sqrt[5])]) y≥1&&1+2 (Sqrt[5]) x≤(Sqrt[5])+2 (Sqrt[5+2 (Sqrt[5])]) y&&2 (Sqrt[5+2 (Sqrt[5])]) y≤1&&5 x+3 (Sqrt[5]) x+(Sqrt[2 (5+(Sqrt[5]))]) y≤2)))||(1+2 (Sqrt[5]) x≥(Sqrt[5])+2 (Sqrt[5+2 (Sqrt[5])]) y&&2+(5+(Sqrt[5])) x+(Sqrt[50+22 (Sqrt[5])]) y≥0&&(5+3 (Sqrt[5])) x+(Sqrt[2 (5+(Sqrt[5]))]) y≤2)||((Sqrt[(2/(65+29 (Sqrt[5])))])+(Sqrt[5]) y≥3 y&&(((5+3 (Sqrt[5])) x+(Sqrt[2 (5+(Sqrt[5]))]) y≥2&&-2 (Sqrt[(2/(65+29 (Sqrt[5])))])+(Sqrt[50-22 (Sqrt[5])]) x+y≤(Sqrt[5]) y)||(2 (Sqrt[(2/(65+29 (Sqrt[5])))])+(Sqrt[50-22 (Sqrt[5])]) x+(Sqrt[5]) y≥y&&2+(5+3 (Sqrt[5])) x≤(Sqrt[2 (5+(Sqrt[5]))]) y)))||(2 (Sqrt[5+2 (Sqrt[5])]) y≥1&&(Sqrt[650-290 (Sqrt[5])])+10 (Sqrt[5-2 (Sqrt[5])]) x≥10 (-2+(Sqrt[5])) y&&10 ((Sqrt[5-2 (Sqrt[5])]) x+(-2+(Sqrt[5])) y)≤(Sqrt[650-290 (Sqrt[5])]))

nineq[x_,y_]:=(2.`+11.70820393249937` x≥3.804226065180614` y&&((2.`+9.95959313953112` y≥7.23606797749979` x&&1.2360679774997898`+4.47213595499958` x+6.155367074350507` y≤0)||(1.2360679774997898`+4.47213595499958` x+6.155367074350507` y≥0&&4.47213595499958` x≤1.2360679774997898`+6.155367074350507` y&&6.155367074350507` y≤1.`&&11.70820393249937` x+3.804226065180614` y≤2.`)))||(4.47213595499958` x≥1.2360679774997898`+6.155367074350507` y&&2.`+7.23606797749979` x+9.95959313953112` y≥0&&11.70820393249937` x+3.804226065180614` y≤2.`)||(0.7639320225002102` y≤0.12410828034667905`&&((11.70820393249937` x+3.804226065180614` y≥2.`&&0.8980559531591694` x≤0.2482165606933581`+1.2360679774997898` y)||(0.2482165606933581`+0.8980559531591694` x+1.2360679774997898` y≥0&&2.`+11.70820393249937` x≤3.804226065180614` y)))||(6.155367074350507` y≥1.`&&1.241082803466792`+7.265425280053607` x≥2.360679774997898` y&&7.265425280053607` x+2.360679774997898` y≤1.241082803466792`)

p=Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^4}];//Timing

((Graphics`Mesh`PolygonArea[poly])/((x2-x1)(y2-y1)))

(p2=Select[p,ineq@@#&])//Length//Timing

(p2=Select[p,nineq@@#&])//Length//Timing

(p2=Select[p,Graphics`Mesh`InPolygonQ[poly,#]&])//Length//Timing

(p2=Select[p,Graphics`Mesh`InPolygonQ[poly,#]&])//Length//Timing

Graphics[{poly,Green,Point[p2]}]

lines=Line/@Partition[p2,2];

Take[lines,10]

Graphics[{poly,{If[Graphics`Mesh`IntersectQ[Graphics[{poly,#}],AllPoints->False],Red,Green],#}&/@lines
}]

good=Count[lines,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,#}],AllPoints->False]&)]

good/Length[lines]//N

(10^4)

Length[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^4}],nineq@@#&],2]]//Timing

Count[l,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly//N,Line[#]}],AllPoints->False]&),{1}]/Length[l]//N//Timing

(10^5)

Length[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^5}],nineq@@#&],2]]//Timing

Count[l,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly//N,Line[#]}],AllPoints->False]&),{1}]/Length[l]//N[#,10]&//Timing

(Length[l]-Count[l,_?(Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}])/Length[l]//N[#,10]&//Timing

c1=Count[l,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}]

c2=Count[l,_?(Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}]

{c1+c2,Length[l]}

(10^6) ineq

Length[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^6}],nineq@@#&],2]]//Timing

Count[l,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}]/Length[l]//N[#,10]&//Timing

MaxMemoryUsed[]/10^6.

(10^7)

Length[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^7}],nineq@@#&],2]]//Timing

Count[l,_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}]/Length[l]//N[#,10]&//Timing

(1./(Sqrt[Length[l]]))

(10^8)

Length[l=]//Timing

c=0;
n=0;
t=AbsoluteTime[];
Do[
c+=Count[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^6}],nineq@@#&],2],_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}];
n+=Length[l];
Print[{i,AbsoluteTime[]-t,c,n,c/n//N}],
{i,100}]

{o,dt,c,n,f}={192,11729.210294`10.04369251163323,27393641,31318103,0.8746903029216042`};
t=AbsoluteTime[]-dt;
Do[
c+=Count[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^6}],nineq@@#&],2],_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}];
n+=Length[l];
Print[{i+o,AbsoluteTime[]-t,c,n,c/n//N}],
{i,50}]

Quit

{o,dt,c,n,f}={242,14811.170232`10.020074375694094,34529582,39474061,0.8747410609716594`};
t=AbsoluteTime[]-dt;
Do[
c+=Count[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^6}],nineq@@#&],2],_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}];
n+=Length[l];
Print[{i+o,AbsoluteTime[]-t,c,n,c/n//N}],
{i,50}]

{o,dt,c,n,f}={292,17544.967505`9.996727557356174,41664265,47632161,0.8747086868471073`};
t=AbsoluteTime[]-dt;
Do[
c+=Count[l=Partition[Select[Table[{RandomReal[{x1,x2}],RandomReal[{y1,y2}]},{10^6}],nineq@@#&],2],_?(!Graphics`Mesh`IntersectQ[Graphics[{poly,Line[#]}],AllPoints->False]&),{1}];
n+=Length[l];
Print[{i+o,AbsoluteTime[]-t,c,n,c/n//N}],
{i,50}]

Quit

Exact

Annulus

As[y_,θ_]:=(((a^2)-(b^2))((π/2)-ArcSin[(b/y)])+(a^2)ArcCos[(b/a)]-b Sqrt[(a^2)-(b^2)])

Assuming[0<b<a,FullSimplify[(1/(((π((a^2)-(b^2))))^2))Integrate[As[y,θ]y ,{y,b,a},{θ,0,2π}]]]//Timing

expr1=((-2 b (Sqrt[(a^2)-(b^2)])+2 (a^2) ArcCos[(b/a)])/((a^2) π-(b^2) π));

expr2=((π (a^2)-2a b(Sqrt[1-((((b/a)))^2)])-2b(Sqrt[(a^2)-(b^2)])+2(a^2)(ArcCos[(b/a)]-ArcSin[(b/a)]))/(2π((a^2)-(b^2))));

FullSimplify[expr1-expr2,0<b<a]

FullSimplify[expr2,0<b<a]

Limit[expr1,b->0]

Pentagram

A=(1/4)(Sqrt[650-290(Sqrt[5])]);
A0=With[{n=5,a=(Sqrt[5])-2},(1/4)n (a^2)Cot[(π/n)]];
A1=With[{b=(Sqrt[5])-2,ρ=(Sqrt[(1/10)(5-(Sqrt[5]))]),r=(1/2)(Sqrt[(1/5)(5-2(Sqrt[5]))])},(1/2)b(ρ-r)];

N[{A,A0,A1}]

FullSimplify[A-(A0+5A1)]