7
$\begingroup$

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

enter image description here

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

enter image description here

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.

$\endgroup$
1
  • 1
    $\begingroup$ Look at Thinning $\endgroup$
    – yarchik
    Jun 19 at 9:23

4 Answers 4

10
$\begingroup$

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

enter image description here

  • Just for fun
reg = BoundaryDiscretizeRegion[
   BoundaryDiscretizeGraphics[
    Text[Style[ω, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
radius = Abs@
   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", 
          RegionNearest[reg][RegionCentroid[reg]]]],.5*radius]}, {x, y},
      Method -> Automatic]}, 
  Graphics[{reg, Red, Disk[{x, y} /. sol[[2]], Abs[sol[[1]]]], Blue, 
    AbsolutePointSize[3], Point[{x, y} /. sol[[2]]]}]]]

enter image description here

reg = BoundaryDiscretizeRegion[
   BoundaryDiscretizeGraphics[
    Text[Style[ω, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
radius = Abs@
   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]]]]}]]

enter image description here

reg = BoundaryDiscretizeRegion[
   BoundaryDiscretizeGraphics[
    Text[Style[℘, Bold, FontFamily -> "Times"]], _Text]];
dist = SignedRegionDistance[reg];
radius = Abs@
   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}]

enter image description here

$\endgroup$
2
  • $\begingroup$ Thanks. Is there a way to get the curve as a region, or as a set, instead of as an image. $\endgroup$ Jun 19 at 16:01
  • $\begingroup$ +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. :) $\endgroup$
    – Michael E2
    Jun 20 at 0:30
7
$\begingroup$

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

Addendum

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.

$\endgroup$
2
  • $\begingroup$ (+1)Wonderful ! $\endgroup$
    – cvgmt
    Jun 22 at 0:23
  • 1
    $\begingroup$ Simpler, but slower:Show[reg,Epilog->{Select[MeshPrimitives[vor,1],RegionWithin[reg,#]&] }] $\endgroup$
    – yode
    Jun 22 at 5:38
6
$\begingroup$
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}]

enter image description here

An attempt to get points

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

enter image description here

enter image description here

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

$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ Jun 19 at 13:18
4
$\begingroup$
reg = BoundaryDiscretizeRegion[
  BoundaryDiscretizeGraphics[
   Text[Style[\[Omega], Bold, FontFamily -> "Times"]], _Text]]
prange = CoordinateBounds[reg, .3];
s = Show[reg, PlotRange -> prange];

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

enter image description here

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]

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.