# Is there a numerical method/built-in to calculate the boundary of a set of graphs?

Recently, I encounted a geometry problem in my work. For a curve-family that owns the following parametric equation: $$E(t,\theta)= \begin{pmatrix} x_E(t,\theta) \\ y_E(t,\theta) \end{pmatrix}$$

where, $\theta \in [0,2\pi]$ and $t\in [0,1]$

Traditionally, I could apply the envelope theory to solve the points that located in the boundary.

$$\frac{\partial x_E(t,\theta)}{\partial t}\frac{\partial y_E(t,\theta)}{\partial \theta}-\frac{\partial x_E(t,\theta)}{\partial \theta}\frac{\partial y_E(t,\theta)}{\partial t}=0 \qquad (1)$$

So for the fixed $t_i$, the parameter of envelope-point $\theta_i$ could be solved with equation$(1)$

Here is a normal instance that using above theory(denoted envelope-point with $\color{blue} \square$).

Obvisouly, I can connect $P_{0,1}, P_{1,1}, \dots P_{4,1}$ and $P_{0,2},P_{1,2},\dots P_{4,2}$ in sequence to achieve the boundary/envelope.

However, for the complicated case, there are only some envelope-points(denoted with $\color{blue} \square$) on the boundary/envelope. That is, the envelope theory will unapplicable.

So my question is:

• Is there a numerical method/built-in to calculate the boundary?(and achieve the coordinates of points that located on the boundary)

### Update

Here are some data

coeff = {{{0., -5., 0}, {-5.2203, 0., 1.7945}},
{{-0.4188, -4.9846, 0.1071}, {-5.3218, 0.3923, 2.0267}},
{{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}},
{{-1.3234, -4.8618, 0.2192}, {-5.5122, 1.1672, 2.6475}},
{{-1.8203, -4.7553, 0.2473}, {-5.6022, 1.5451, 3.0486}},
{{-2.3568, -4.6194, 0.2742}, {-5.6897, 1.9134, 3.5173}},
{{-2.9427, -4.455, 0.3147}, {-5.7755, 2.27, 4.0578}},
{{-3.5912, -4.2632, 0.3857}, {-5.8604, 2.6125, 4.6738}},
{{-4.3197, -4.0451, 0.5068}, {-5.9456, 2.9389, 5.368}},
{{-5.1524, -3.802, 0.7017}, {-6.0327, 3.2472, 6.1428}},
{{-6.1237, -3.5355, 1.}, {-6.1237, 3.5355, 7.}}};

coeff2 = {{{0., -5., 0}, {-5.2203, 0., 1.7945}},
{{-0.4188, -4.9846, 0.3754}, {-5.3218, 0.3923, 1.8307}},
{{-0.8583, -4.9384, 0.6792}, {-5.4189, 0.7822, 1.8663}},
{{-1.3234, -4.8618, 0.9146}, {-5.5122, 1.1672, 1.9093}},
{{-1.8203, -4.7553, 1.0855}, {-5.6022, 1.5451, 1.9672}},
{{-2.3568, -4.6194, 1.1959}, {-5.6897, 1.9134, 2.047}},
{{-2.9427, -4.455, 1.2502}, {-5.7755, 2.27, 2.1556}},
{{-3.5912, -4.2632, 1.2528}, {-5.8604, 2.6125, 2.2995}},
{{-4.3197, -4.0451, 1.2087}, {-5.9456, 2.9389, 2.4846}},
{{-5.1524, -3.802, 1.1229}, {-6.0327, 3.2472, 2.7164}},
{{-6.1237, -3.5355, 1.}, {-6.1237, 3.5355, 3.}}};


which are the coefficient of ellipse. Namely, {{a,b,c},{d,e,f}}

$\begin{cases} x=a \sin\theta+b \cos\theta +c \\ y=d \sin\theta +e \cos\theta +f \\ \end{cases}$

ellipsePoints[{mat1_, mat2_}] :=
{mat1.{Sin[#], Cos[#], 1},
mat2.{Sin[#], Cos[#], 1}} & /@ Range[0, 2 Pi, 0.02 Pi]

points = Flatten[ellipsePoints /@ coeff, 1];
points2 = Flatten[ellipsePoints /@ coeff2, 1];


Thanks for RunnyKine's alphaShapes2D[] with diferent threshold :1,3

reg = RegionBoundary@alphaShapes2D[points, 1];
Show[{reg, ListPlot[points, AspectRatio -> Automatic]}, Axes -> True]

reg2 = RegionBoundary@alphaShapes2D[points, 3];
Show[{reg2, ListPlot[point2s, AspectRatio -> Automatic]}, Axes -> True]


• Would you ever have a case where your graphs don't form a single connected region? (If so, do you want to find two separate boundaries in that case?) Commented Mar 17, 2016 at 9:43
• @MartinBüttner No, in my case, the graphs always form a single connected regoin.
– xyz
Commented Mar 17, 2016 at 10:01

Here is another answer inspired by Rahul's answer that also uses only built-in functions:

RegionBoundary @ DiscretizeGraphics @ Graphics[Polygon /@ ellipsePoints /@ coeff]


RegionBoundary@
DiscretizeGraphics@Graphics[Polygon /@ ellipsePoints /@ coeff2]


RegionBoundary@
DiscretizeGraphics@
Graphics[Polygon /@
Table[
Table[
RotationMatrix[m].{2 + 5 Cos[x], 3 + 6 Sin[x]},
{x, 0, 2 Pi, 0.02 Pi}], {m, 0, Pi, Pi/20}]]


• Hi, RunnyKine. I would like to know that is it possible to know the theory/algorithm behind the built-in RegionBoundary[] and DiscretizeGraphics[]. THX:)
– xyz
Commented Apr 16, 2016 at 8:45
• @ShutaoTANG. I would like to know that myself. Commented Apr 16, 2016 at 9:17
• Maybe the theory/algorithm belong to the part of computational geometry. However, I have no knowledge about that field.
– xyz
Commented Apr 16, 2016 at 9:23

Here is one way to do it, although it gives you a rasterised result.

graph = Table[
ParametricPlot[
RotationMatrix[m].{2 + 5 Cos[x], 3 + 6 Sin[x]}, {x, 0, 2 Pi},
PlotRange -> {{-10, 10}, {-10, 10}}, Axes -> None], {m, 0, Pi,
Pi/30}] // Flatten // Show

binary = graph // Binarize;
binary // ColorNegate // DeleteBorderComponents,
binary
] // EdgeDetect // ColorNegate


### How it works

First, we render the graph without the axes (like you already did):

Then we binarise it and store that in binary:

We want to find the bounary of the outer black region. To do so, we first negate the binary image and remove the outer region with DeleteBorderComponents (this is essentially a flood fill with black from the edges):

If we add these two images together, both the lines of the graph as well as all the regions it encloses will become white, because the lines are white in the original binary image and the inner regions are white in the negated image. However, the outer region is black in both so it remains black:

Now it's trivial to detect the edge. Because EdgeDetect shows edges in white on black, we're negating the image again to get black on white:

I'm looking forward to someone coming up with a solution that gives vectorised output!

• Thanks for your solution:). So I would like to know is it possible to achieve the coordinates of the boundary. Because I need to use these coordinates to do other thing.
– xyz
Commented Mar 17, 2016 at 10:03
• @ShutaoTANG Not easily with the rasterised image. You'd need to figure out the scale of the pixels and then you could just check which pixels are black. That's why I'm hoping someone else will come up with a solution in vectorised form. Please edit your question though to clarify exactly what information you need about the boundary. Commented Mar 17, 2016 at 10:05

Since you require the coordinates, a good starting point is to use my alphaShapes2D code from here. Now it's just a matter of the following one-liner:

breg = First @ ConnectedMeshComponents @ RegionBoundary @ alphaShapes2D[pts, 0.5]


Here pts is as defined in the OP. Since you want the boundary coordinates. you can obtain it through Meshcoordinates:

MeshCoordinates @ breg

• @MartinBüttner this would work: Last[SortBy[ConnectedMeshComponents[reg], RegionMeasure]] Commented Mar 17, 2016 at 13:58
• Thanks. in my laptop, the Mathematica 10.0 gives me this result
– xyz
Commented Mar 17, 2016 at 14:45
• @ShutaoTANG. That's strange. Is that the result you get from RegionBoundary @ alphaShapes2D[pts, 0.5]? Commented Mar 17, 2016 at 18:31
• @RunnyKine Yes, today I ran the code in the other laptop that installed Mathematica 10.3, which still gives me this result
– xyz
Commented Mar 18, 2016 at 1:27
• @ShutaoTANG. OK I just tested in Mathematica 10.2 and I get same result as you. It must be a bug that was fixed in 10.4. Anyways, since you only need the outer boundary, you should still be able to get it with ConnectedMeshComponents[reg][[1]] Commented Mar 18, 2016 at 1:33

Just take the union of the polygons defined by the curves. This uses only built-in functions and doesn't require parameter tuning. I've added RegionBoundary to make it look like the other answers.

RegionBoundary@
BoundaryDiscretizeRegion@
RegionUnion[Polygon /@ ellipsePoints /@ coeff]


If I'm not mistaken, the $\alpha$-shape slightly smooths out the concave corner at the lower left of the figure, while the union doesn't.

Today, I discovered an example from Wolfram Documentation Center about plotting a parametric region:

ParametricPlot[
r^2 { Sqrt[t] Cos[t], Sin[t]}, {t, 0, 3 Pi/2}, {r, 1, 2}]


Namely, ParametricPlot[] could generate the region of family: $f(\theta,t)$

For my question, it is a family of ellipses with respect to time variable $t$.

{δxp, δyp, δzp} = {0, 0, 77.5};
{δxw, δyw, δzw} = {-13.123, 3.631, 42.5};
zOmega = 6;
R = 5;

ParametricPlot[
{R Cos[θC] Cos[θ] + R Sin[θC]/Cos[θA] Sin[θ] +
Sin[θC]/Cos[θA] (yM - δyp) +
Tan[θA] Sin[θC] (zOmega + δzw - δzp) +
Cos[θC] (xM - δxp) + δxp - δxw,

-R Sin[θC] Cos[θ] + R Cos[θC]/Cos[θA] Sin[θ] +
Cos[θC]/Cos[θA] (yM - δyp) +
Tan[θA] Cos[θC] (zOmega + δzw - δzp) -
Sin[θC] (xM - δxp) + δyp - δyw},
{θ, 0, 2 π}, {t, 0, 1}, AspectRatio -> Automatic]


Case 1

{xM, yM, zM, θA, θC} =
{13.123 + 2.9665 t, -13.5318 - 9.5375 t, 65.0342 + 6.201 t,
-0.2915 - 0.6638 t, -3.1416 + 0.7854 t};


Case 2

{xM, yM, zM, θA, θC} =
{13.123 +0.1381 t,-13.5318-5.455 t,65.0342 +5.6236 t,
-0.2915-0.6638 t,-3.1416+0.7854 t}


Case 3

{xM, yM, zM, θA, θC} =
{13.123 -5.5188 t,-13.5318-8.721 t,65.0342 +10.2424 t,
-0.2915-0.6638 t,-3.1416+0.7854 t}


Case 4

{xM, yM, zM, θA, θC} =
{13.123 +4.3807 t,-13.5318-7.9045 t,65.0342 +9.0877 t,
-0.2915-0.6638 t,-3.1416+0.7854 t}


• +1. It doesn't seem easy to extract the boundary from these curves though. Commented Apr 22, 2016 at 9:18
• @RunnyKine Yes, extracting the boundary is hard. Now my thought is converting the first ellipse($t=0$) to a polygon, then using the RegionUnion[] operation.
– xyz
Commented Apr 22, 2016 at 9:23
• Ok, that sounds like a nice approach. Commented Apr 22, 2016 at 9:28
• @RunnyKinen By far, I think your second approach RegionBoundary @ DiscretizeGraphics @ Graphics[Polygon /@ ellipsePoints /@ coeff] is the fastest method, just needing 0.2~0.6 s  to calculate the boundary of $11$ discretized ellispses. In addition, ParametricPlot[] just samples a set of ellipse to achive the region.
– xyz
Commented Apr 22, 2016 at 9:31
• Yes, that approach is quite fast. Commented Apr 22, 2016 at 9:53