Smoother boundaries between components of segmented image

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]


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


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


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:

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,
Frame -> None, DataRange -> {{0, 200}, {0, 200}},
Epilog -> Line[{{0, 0}, {200, 0}, {200, 200}, {0, 200}, {0, 0}}]
]


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

• Strongly related: (1), (2), (3), (4). Commented Aug 27, 2021 at 8:46
• Also related: (5), (6), (7), (8). Commented Aug 27, 2021 at 10:08

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


(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


Applying GaussianFilter allows producing smoother boundaries:

lcp2 = ListContourPlot[
Contours -> {1.5, 2.5}, InterpolationOrder -> 1, ContourStyle -> {Blue, Darker@Green},
BaseStyle -> Thick, ContourShading -> None, Frame -> None, DataRange -> All] //
EchoTiming

33.9752


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

lcp2c = ListContourPlot[
Contours -> {1.9, 2.1}, InterpolationOrder -> 1, ContourStyle -> {Blue, Darker@Green},
BaseStyle -> Thick, ContourShading -> None, Frame -> None, DataRange -> All] //
EchoTiming

31.8284


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

lcp2cg = ListContourPlot[
Contours -> {1.9, 2.1}, InterpolationOrder -> 2, Method -> "Spline",
ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None,
Frame -> None, DataRange -> All] // EchoTiming

118.176


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

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

lcp3 = ListContourPlot[
Contours -> {1.9, 2.2}, InterpolationOrder -> 2, Method -> "Spline",
ContourStyle -> {Blue, Darker@Green}, BaseStyle -> Thick, ContourShading -> None,
Frame -> None, DataRange -> All] // EchoTiming

306.132


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


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


• 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.
Commented Sep 1, 2021 at 1:25

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]


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

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]


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


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]


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


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


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


(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*)