MedialGraph5[L_] := 
Module[{g, h, i, j, k, m, n, x, eps, intervalValues, bNode, vlabel, intlabel, 
fixedValues, geodesicVertices, usedVertices, firstValues, firstPairs, 
geodesicPaths, lastLine}, 

(

Off[Transpose::nmtx];

Off[NullSpace::matrix];

Off[General::stop];

modulo[x_,m_] := If[Mod[x,m] == 0, x, Mod[x,m]];

Attributes[modulo]={Listable};

RofA[i_,j_] := j - i - Floor[j/2] + Floor[(i+1)/2] - 1 - 
( If[Range[Floor[j/2] + 1, Floor[(i+1+2*Length[L])/2] - 1] == {},
0, 
(Length[Transpose[L[[modulo[Range[Floor[i/2] + 1, Floor[(j+1)/2] - 1], 
Length[L]], modulo[Range[Floor[j/2] +1, Floor[(i + 1 + 2*Length[L])/2] - 1],
Length[L]]]]]] 
- Length[NullSpace[L[[modulo[Range[Floor[i/2] + 1, 
Floor[(j+1)/2] - 1], Length[L]], modulo[Range[Floor[j/2] + 1, 
Floor[(i + 1 + 2*Length[L])/2] - 1], Length[L]]]]]])]);
 
intervalValues = Table[RofA[i,j], {i,1,2*Length[L]}, {j, i+1, 2*Length[L]}];

bNode = Table[{PointSize[0.025], Point[{Cos[(2*k*Pi)/Length[L]],
- Sin[(2*k*Pi)/Length[L]]}]}, {k, 0, Length[L] - 1}];

vlabel = Table[
Text[Subscript[v, k+1] , {1.1*Cos[(2*k*Pi)/Length[L]], -1.1*Sin[(2*k*Pi)/Length[L]]}],
{k, 0, Length[L] - 1}];

intlabel = Table[
Text[Subscript[I,k+2], {1.3*Cos[2*k*Pi/(2*Length[L])], -1.3*Sin[2*k*Pi/(2*Length[L])]}], 
{k, -1, 2*Length[L] - 2}];

geoV[n_] := Table[{Cos[(2*Pi*Floor[k/2] +eps*(-1)^(k+1))/n],
-Sin[(2*Pi*Floor[k/2] + eps*(-1)^(k+1))/n]} 
/. eps -> Random[Real, {.5, 1.5}], {k, 0, 2*n - 1}];

fixedValues = geoV[Length[L]];

geodesicVertices = Thread[Point[fixedValues]];

connection[g_] := Table[Part[intervalValues, g, h]
- Part[intervalValues, g+1, h-1] != 0,
{h, 2, Length[intervalValues] - g}];

usedVertices = {};

skipTest[k_] := If[(MemberQ[connection[k], True] && 
FreeQ[usedVertices, Part[fixedValues, k]]), k, skipTest[k+1]]
/. usedVertices -> {AppendTo[usedVertices, Part[fixedValues, k]]};

firstValues = Table[skipTest[k], {k, 1, Length[L] - 1}];

firstPairs = Table[{Part[fixedValues, Part[firstValues, k]],
 Part[fixedValues, Length[Take[connection[Part[firstValues, k]],
  Position[connection[Part[firstValues, k]], True, 1, 1] [[1, 1]] ] ]
   + Part[firstValues, k]]}, {k, 1, Length[L] -1}];

FP=firstPairs;

LP = {Complement[fixedValues, Intersection[fixedValues, 
Flatten[firstPairs, 1]]]};

AP = Union[FP,LP];
 
geodesicPaths = Thread[Line[AP]];

LLt[x_,y_] := Product[y-AP[[k,1,2]]
-(AP[[k,2,2]]- AP[[k,1,2]])*(x-AP[[k,1,1]])/(AP[[k,2,1]]
-AP[[k,1,1]]),{k,1,Length[L]}];

s = Sign[LLt[1,0]];

LL[x_,y_] := s*LLt[x,y];

LLF[x_,y_] := If[x*x+y*y>1, 1, LL[x,y]];

CC=ContourPlot[LLF[x,y], {x,-1,1}, {y,-1,1}, Contours -> {0}, 
PlotPoints -> 100, Frame -> False];

MM=Graphics[{bNode, intlabel, vlabel, geodesicPaths,
{PointSize[0.015], geodesicVertices}, {Thickness[.010], Circle[{0, 0},1]}}, 
AspectRatio -> 1, PlotLabel -> "Medial Graph"];

Show[CC, MM];
)

]

