# Finding "middle curve" inside region

I have some simply connected region, as shown below.

  reg = DiscretizeRegion[BoundaryDiscretizeGraphics[Text[Style[\[Omega], Bold, FontFamily -> "Times"]], _Text]];
ParametricPlot[{x, y}, Element[{x, y}, reg]]


Now, I wish to find the "middle curve", something like the red curve below.

I have no real need to be mathematically precise in the definition, just some curve complete contained in the region that in a vague sense represents the bulk.

My initial idea was to "slice" the region and then in each portion find the middle point, but it's not easy to do it because a straight line will intersect the region twice or thrice.

• Look at Thinning Jun 19 at 9:23

Maybe use Pruning and SkeletonTransform as an alternative。

reg = BoundaryDiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[ω, Bold, FontFamily -> "Times"]], _Text]];
pic = Image[
Pruning[SkeletonTransform[ColorNegate@Graphics@reg,
Method -> Automatic], {35}], ImageSize -> 400];
g = Graphics[{Opacity[.3], reg}, ImageSize -> 400];
Overlay[{ColorNegate@pic, g}]


• Just for fun
reg = BoundaryDiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[ω, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
First@NMinimize[{dist@{x, y}, {x, y} ∈ reg}, {x, y}];
Dynamic[With[{sol =
NMinimize[{dist@{x, y}, {x, y} ∈ reg, {x, y} ∈
Disk[RegionNearest[reg][
MousePosition["Graphics",
Method -> Automatic]},
Graphics[{reg, Red, Disk[{x, y} /. sol[[2]], Abs[sol[[1]]]], Blue,
AbsolutePointSize[3], Point[{x, y} /. sol[[2]]]}]]]


reg = BoundaryDiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[ω, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
First@NMinimize[{dist@{x, y}, {x, y} ∈ reg}, {x, y}];
findPoint[{x0_, y0_}] :=
NMinimize[{dist@{x, y}, {x, y} ∈ reg, {x, y} ∈
Disk[RegionNearest[reg][{x0, y0}], .5*radius]}, {x, y},
Method -> Automatic][[2, ;; , 2]];
DynamicModule[{list = {}},
EventHandler[
Dynamic[Graphics[{{Opacity[.2], reg},
Point@list}]], {"MouseClicked" :>
AppendTo[list,
findPoint[
MousePosition["Graphics",
RegionNearest[reg][RegionCentroid@reg]]]]}]]


reg = BoundaryDiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[℘, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
First@NMinimize[{dist@{x, y}, {x, y} ∈ reg}, {x, y}];
findPoint[{x0_, y0_}] :=
NMinimize[{dist@{x, y}, {x, y} ∈ reg, {x, y} ∈
Disk[RegionNearest[reg][{x0, y0}], .2*radius]}, {x, y},
Method -> Automatic];
list = {};
DynamicModule[{},
EventHandler[
Dynamic[Deploy@
Graphics[{{Opacity[.2], reg},
Point[list[[;; , 2, ;; , 2]]]}]], {"MouseMoved" :>
AppendTo[list,
findPoint[
MousePosition["Graphics",
RegionNearest[reg][MousePosition["Graphics"]]]]]}]]


After collect all the data in list, we draw the animation.

list;
ani = Animate[
Graphics[{reg, Red,
Table[Disk[{x, y} /. list[[i]][[2]], list[[i]][[1]] // Abs], {i,
1, k}]}, PlotRange -> 4], {k, 1, Length@list}]


• Thanks. Is there a way to get the curve as a region, or as a set, instead of as an image. Jun 19 at 16:01
• +1 for saving me the time of trying to figure out how to do the red disk thing. I don't think I'd ever get to it. :) Jun 20 at 0:30

We can approximate the medial axis with the interior edges of the mesh coordinate's Voronoi diagram. Sort of similar to the answer here, but in 2D.

reg = BoundaryDiscretizeGraphics[Text[Style[ω, Bold, FontFamily -> "Times"]], _Text];

vor = VoronoiMesh[MeshCoordinates[reg]];

keepQ = Thread[And @@ RegionMember[reg] /@ Transpose[MeshPrimitives[vor, 1][[All, 1]]]];

keepinds = Pick[Range[MeshCellCount[vor, 1]], keepQ];

medialaxis = MeshRegion[MeshCoordinates[vor], MeshCells[vor, {1, keepinds}]];

Show[
reg,
MeshRegion[medialaxis, MeshCellStyle -> {0 -> None, 1 -> Black}]
]


Note that this is a geometrically unstable approach and results can vary wildly depending on the discretization of the region. Based on papers of Amenta and Attali, a better approximation can be found, but we need to be able to find intersections of spheres efficiently.

• (+1)Wonderful ! Jun 22 at 0:23
• Simpler, but slower：Show[reg,Epilog->{Select[MeshPrimitives[vor,1],RegionWithin[reg,#]&] }]
– yode
Jun 22 at 5:38
reg = DiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[\[Omega], Bold, FontFamily -> "Times"]], _Text]];
pp = ParametricPlot[{x, y}, Element[{x, y}, reg],
ImageSize -> {300, 200}];
th = Colorize[Thinning[MorphologicalComponents[pp]],
ColorRules -> {0 -> White, _ -> Black}, ImageSize -> {300, 200}];
Rasterize@Overlay[{th, pp}]


An attempt to get points

Getting points is more complicated. Here is a ListPlot of my unsuccessful, but interesting attempt:

My code is all over the place and I will not include it, but the approach started with the following:

mc = MorphologicalComponents[Graphics[reg]];
sp = SequencePosition[#, {Repeated[0]}, Overlaps -> False] & /@ mc;


mc is a list with 1 in the background and 0 in the character. Looking for sequence of zeros with SequencePosition returns all the {start, end} positions of the sequence of 0 in a row. This corresponds to the width of the character stroke. What I did is then compute the mean between start and end for each sequence (the x coordinate) and added the y coordinate (row number). The net result is a list of points that can be plotted with ListPlot. Clearly, this works for the vertical strokes in the character but breaks down when we get to horizontal strokes at the base of the character (we would need to measure our zeros in columns at that point).

• Thanks. It looks nice, but, is there a way to smooth out the curve ? I mean, the part in the left has an awkward straight segment. Jun 19 at 13:18
reg = BoundaryDiscretizeRegion[
BoundaryDiscretizeGraphics[
Text[Style[\[Omega], Bold, FontFamily -> "Times"]], _Text]]
prange = CoordinateBounds[reg, .3];
s = Show[reg, PlotRange -> prange];

img = Pruning[
Thinning[
ColorNegate[
Binarize[
LaplacianGaussianFilter[DistanceTransform[ColorNegate@s],
4]]]]], {35}];


take pixel values and create a graph.

pos = PixelValuePositions[img, 1];
g = NearestNeighborGraph[pos, 2];

Counts[VertexDegree[g]]


<|2 -> 1334, 3 -> 6|>

Trim the graph.

deg3 = Pick[VertexList[g], VertexDegree[g], 3 ]


{{589, 447}, {590, 446}, {591, 445}, {238, 438}, {381, 303}, {338,
137}}

HighlightGraph[g, deg3, VertexLabels -> Automatic]


NeighborhoodGraph[g, {338, 137}, VertexLabels -> Automatic]


w = First[
ConnectedGraphComponents[
EdgeDelete[
VertexDelete[g,
Most[deg3] ], {338, 137} \[UndirectedEdge] {339, 138}]]]


Convert it to a mesh region.

cc = GraphEmbedding[w];
idim = ImageDimensions[img];
x = Rescale[cc[[All, 1]], {0, idim[[1]]}, prange[[1]]];
y = Rescale[cc[[All, 2]], {0, idim[[2]]}, prange[[2]]];

res = MeshRegion[Transpose[{x, y}], Line[Range[VertexCount[w]]]];

Show[{reg, res}]