8
$\begingroup$

I'd like to use boundaries between components of a segmented image to create a mesh for finite elements, as described in Mesh for Images for three materials and meshes with multiple regions from 2D images?.

The goal is smooth boundaries that roughly approximate the edges of the segmented components: the exact component boundaries, specified by individual pixels are finer than justified by the blur or noise in the original image.

To create a mesh, the boundaries must not have any gaps. One approach is using the "Contours" property of the components. But these boundaries be very jagged, which leads to unnecessarily fine meshes.

Other approaches are extracting boundaries from ImageMesh of each component, or from a contour plot. Both methods give smoother boundaries, but can have gaps or duplicates.

Is there a way to get boundaries without gaps, as given by the Contours property, but smoother, as given by the ImageMesh or contour plot approaches?

The boundaries should be a form suitable for meshing, e.g., line segments.

Example

Create an image:

img = Blur[Rasterize[Graphics[{Gray, Rectangle[], Blue, Annulus[]}, 
    PlotRangePadding -> None], RasterSize -> 200], 5]

example image for component boundaries

Segment the image (clustering components is an example, the actual segmentation could come from other methods):

components = ClusteringComponents[img, 3];

Method 1: use Cluster property of the segments

boundaries = 
  ComponentMeasurements[components, "Contours", All, 
   "ComponentAssociation"];

Show[Values[
  Graphics[#, PlotRange -> {{0, 200}, {0, 200}}, 
     PlotRangePadding -> 10] & /@ boundaries]]

boundaries from Contours property

The boundaries between components are jagged. Also, each such boundary appears twice, once for each of the two neighboring components. Since the two boundaries use exactly the same points, that duplication doesn't cause a problem with meshing.

Method 2: use ImageMesh boundary of each segment

regions = 
  Table[ImageMesh[
    Image[components /. Thread[DeleteCases[Range[3], i] -> 0]]], {i,3}];

Graphics[{EdgeForm[Black],Riffle[{Red, LightBlue, Green}, regions]}, 
 PlotRangePadding -> None]

ImageMesh of components

The boundaries of these regions are smoother than those from the components "Contours" property. But there are slight gaps between the regions, so their boundaries don't quite match, as seen in this detail:

detail of ImageMesh of components

Method 3: use a contour plot

Create a contour plot on a blurred version of the components, with contour values selected to be between the gray levels specified for the components.

imgC = ColorConvert[
   Blur[Colorize[components, 
     ColorRules -> {1 -> White, 2 -> Black, 3 -> Gray}], 5], 
   "Grayscale"];

ListContourPlot[Reverse@ImageData@imgC, Contours -> {0.4, 0.8}, 
 ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, 
 ContourShading -> None, 
 Frame -> None, DataRange -> {{0, 200}, {0, 200}},
 (* add outer boundary *)
 Epilog -> Line[{{0, 0}, {200, 0}, {200, 200}, {0, 200}, {0, 0}}]
 ]

boundaries from ListContourPlot

The boundaries are smooth, but with two slightly different boundaries between some of the components.

$\endgroup$
2
4
$\begingroup$

One simple approach is to upsize the original image using the "Lanczos" resampling method (the default for ImageResize) and then apply the Douglas-Peucker algorithm to remove the jaggies and simplify the curve (see the linked answer for the definition of dp):

components = ClusteringComponents[ImageResize[img, Scaled[3], Resampling -> "Lanczos"], 3];

boundaries = ComponentMeasurements[components, "Contours", All][[All, 2]];

eps = 1.5;
Graphics[boundaries /. Line[pts_] :> Line[dp[pts, eps]]]

out

(I've tried playing with other resampling methods, but hadn't found anything better than "Lanczos".)

Alternative method of producing the boundaries using ListContourPlot (addition of InterpolationOrder -> 2 and Method -> "Spline" in this case doesn't improve the result):

components = ClusteringComponents[img, 3];
lcp = ListContourPlot[Reverse[ArrayPad[components, 2, 1]], Contours -> {1.5, 2.5}, 
  InterpolationOrder -> 1, ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, 
  ContourShading -> None, Frame -> None, DataRange -> All] // EchoTiming
56.9843

out

Applying GaussianFilter allows producing smoother boundaries:

lcp2 = ListContourPlot[
   GaussianFilter[N@ArrayPad[Reverse[components], 4, 1], 2, Padding -> None], 
   Contours -> {1.5, 2.5}, InterpolationOrder -> 1, ContourStyle -> {Blue, Darker@Green}, 
   BaseStyle -> Thick, ContourShading -> None, Frame -> None, DataRange -> All] // 
  EchoTiming
33.9752

out

By adjusting the contour values we can get the boundaries close to each other:

lcp2c = ListContourPlot[
   GaussianFilter[N@ArrayPad[Reverse[components], 4, 1], 2, Padding -> None], 
   Contours -> {1.9, 2.1}, InterpolationOrder -> 1, ContourStyle -> {Blue, Darker@Green}, 
   BaseStyle -> Thick, ContourShading -> None, Frame -> None, DataRange -> All] // 
  EchoTiming
31.8284

out

Increasing the radius for GaussianFilter allows to get even smoother boundaries:

lcp2cg = ListContourPlot[
   GaussianFilter[N@ArrayPad[Reverse[components], 6, 1], 3, Padding -> None], 
   Contours -> {1.9, 2.1}, InterpolationOrder -> 2, Method -> "Spline", 
   ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None, 
   Frame -> None, DataRange -> All] // EchoTiming
118.176

out

Swapping labels 2 and 3 gives other boundaries, what may be useful:

components2 = components /. {2 -> 3, 3 -> 2};

lcp3 = ListContourPlot[
   GaussianFilter[N@ArrayPad[Reverse[components2], 6, 1], 4, Padding -> None], 
   Contours -> {1.9, 2.2}, InterpolationOrder -> 2, Method -> "Spline", 
   ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None, 
   Frame -> None, DataRange -> All] // EchoTiming
306.132

out

The same approach can be applied to the original image data directly, without using ClusteringComponents:

data = ArrayPad[Map[Total, Reverse@ImageData@img, {2}], 2, 3];
lcp4 = ListContourPlot[data, Contours -> {1.25, 2}, InterpolationOrder -> 1, 
   ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None, 
   Frame -> None, DataRange -> All] // EchoTiming
150.602

out

By adjusting the contour values we can get the boundaries close to each other:

data = ArrayPad[Map[# . {4., 4., 0.5} &, Reverse@ImageData@img, {2}], 4, 8.5];
lcp5 = ListContourPlot[data, Contours -> {4, 4.5}, InterpolationOrder -> 1, 
   ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None, 
   Frame -> None, DataRange -> All] // EchoTiming
149.305

out

$\endgroup$
1
  • $\begingroup$ Thanks for the variety of useful approaches. Good to have alternatives in case one method doesn't work well for some images. The pointer to the D-P algorithm is especially helpful: simple and effective. $\endgroup$
    – tad
    Sep 1 at 1:25
3
$\begingroup$

Update

In the Wolfram Function Repository now is available CurveToBSplineFunction. It is a significantly extended and improved version of the function published in the original version of this answer. It can be used as follows:

components = ClusteringComponents[img, 3];
(* remove backgound (the component 1)*)
boundaries = ComponentMeasurements[components, "Contours"][[2 ;;, 2]];

(* maximal polynomial degree is 90 *)
functs = ResourceFunction["CurveToBSplineFunction"][#, 90,
      "CurveClosed" -> True][t] & @@@ Flatten[boundaries, 1];
ParametricPlot[functs, {t, 0, 1}, Axes -> False]

output

Please see additional explanations and examples of use in the online Documentation.


Original answer

Here I present my own implementation of the contour parametrization approach: a function contourToBSplineFunction. It can automatically assign weights to the original points based on how far this point is from adjacent points. This approach allows preserving the sharp corners of the original contour. The function also allows specifying the maximal degree of the underlying polynomial basis.

Here is how it can be used:

components = ClusteringComponents[img, 3];
(* remove backgound (the component 1)*)
boundaries = ComponentMeasurements[components, "Contours"][[2 ;;, 2]];

(* automatic weghts for the control points (preserve sharp corners), maximal polynomial degree is 90 *)
functs = contourToBSplineFunction[#, 90][t] & @@@ Flatten[boundaries, 1];
ParametricPlot[functs, {t, 0, 1}, Axes -> False]

out

(* equal weghts for the knots, maximal polynomial degree is 90 *)
functs = contourToBSplineFunction[#, 90, False][t] & @@@ Flatten[boundaries, 1];
ParametricPlot[functs, {t, 0, 1}, Axes -> False]

out

The fourth parameter allows to specify the scale (by default assumed scale 1 corresponding to the distance 1 between adjacent pixels on the plot):

functs = contourToBSplineFunction[#, 90, True, 10][t] & @@@ Flatten[boundaries, 1];
ParametricPlot[functs, {t, 0, 1}, Axes -> False]

out

Parameterizing an arbitrary polygon:

SeedRandom[1];
pts = RandomPolygon[{"Simple", 30}][[1]];
funct = contourToBSplineFunction[pts, 10, True, .001][t];
ParametricPlot[funct, {t, 0, 1}, PlotStyle -> Red, Prolog -> Polygon[pts]]

out

Definition:

ClearAll[contourToBSplineFunction]
contourToBSplineFunction[points_?MatrixQ, maxDegree_Integer : 50, weights_ : True, 
 scale : _?NumberQ : 1] := 
  Block[{pairs, pair2ListW, pairsW},
   pairs = Partition[points, 2, 1, {1, 1}];
   pair2ListW[{{x1_, y1_}, {x1_, y1_}}] = Nothing;
   pair2ListW[{{x1_, y1_}, {x2_, y2_}}] := 
    Block[{dist = Norm[{x2 - x1, y2 - y1}], n}, n = Floor[dist/scale]; 
     Join[{{{x1, y1}, dist}}, 
      Table[{{x1, y1} + i {x2 - x1, y2 - y1}/n, dist/n}, {i, 1, n - 1}], 
      {{{x2, y2}, dist}}]];
   pairsW = Flatten[pair2ListW /@ pairs, 1];
   BSplineFunction[pairsW[[All, 1]], SplineDegree -> maxDegree, 
    SplineWeights -> If[TrueQ@weights, pairsW[[All, 2]], Automatic], 
    SplineClosed -> True]];
$\endgroup$
1
$\begingroup$

Here is a curve parametrization approach based on the code by Simon Woods and its adaptation by Michael E2 for the cases like this. In our case it requires further adaptation because the contours returned by ComponentMeasurements aren't suitable "as is" for Simon's tocurve function.

components = ClusteringComponents[ImagePad[img, 4, White], 3];

Preliminarily, we must convert our contours to the form suitable for Simon's tocurve function.

First, we delete consecutive duplicate points (and remove the background component 1):

boundaries = ComponentMeasurements[components, "Contours"][[2 ;;, 2]] //.
   Line[{a___, p_, p_, b___}] :> Line[{a, p, b}];
gr = Graphics[Riffle[{Blue, Darker@Green}, boundaries]]

out

(It is really strange that ComponentMeasurements simplifies the contours in the horizontal and vertical directions, but at the same time creates many consecutive duplicate points...)

Secondly, we fill the vertical and horizontal gaps between successive points in order to make the curves suitable for Fourier:

boundaries4Fourier = 
  boundaries //. {Line[{a___, {x1_, y1_}, {x1_, y2_}, b___}] /; Abs[y2 - y1] > 1 :> 
     Line[{a, Sequence @@ Table[{x1, y}, {y, y1, y2, Sign[y2 - y1]}], b}], 
    Line[{a___, {x1_, y1_}, {x2_, y1_}, b___}] /; Abs[x2 - x1] > 1 :> 
     Line[{a, Sequence @@ Table[{x, y1}, {x, x1, x2, Sign[x2 - x1]}], b}]};

Simon's function:

param[x_, m_, t_] := 
  Module[{f, n = Length[x], nf}, 
   f = Chop[Fourier[x]][[;; Ceiling[Length[x]/2]]];
   nf = Length[f];
   Total[(2 Abs[f]/Sqrt[n] Sin[
        Pi/2 - Arg[f] + 2. Pi Range[0, nf - 1] t])[[;; Min[m, nf]]]]];

tocurve[Line[data_], m_, t_] := param[#, m, t] & /@ Transpose[data];

Michael E2's code

pts = Flatten[boundaries4Fourier, 1][[All, 1]];
divider = 9;(*controls smoothing*)
pp = MapThread[(*parametrized paths*)
  tocurve[Line[#1], #2, t] &, {pts, 
   Max[#, 4] & /@ Round[(Length /@ pts)/divider](*no.of modes-controls smoothing*)}];
pp = pp - (pp /. _Sin -> 0) + (Mean /@ pts);
plot = ParametricPlot[pp, {t, 0, 1}, Axes -> False]

out

Not ideal, but fast and smooth.

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