17
$\begingroup$

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

enter image description here

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.

enter image description here

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]

enter image description here

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ 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?) $\endgroup$
    – Martin Ender
    Commented Mar 17, 2016 at 9:43
  • $\begingroup$ @MartinBüttner No, in my case, the graphs always form a single connected regoin. $\endgroup$
    – xyz
    Commented Mar 17, 2016 at 10:01

5 Answers 5

Reset to default
7
$\begingroup$

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

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

Mathematica graphics

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

enter image description here

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

enter image description here

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

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;
boundary = ImageAdd[
    binary // ColorNegate // DeleteBorderComponents,
    binary
] // EdgeDetect // ColorNegate

How it works

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

enter image description here

Then we binarise it and store that in binary:

enter image description here

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

enter image description here

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:

enter image description here

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:

enter image description here

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

Improve this answer
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – xyz
    Commented Mar 17, 2016 at 10:03
  • $\begingroup$ @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. $\endgroup$
    – Martin Ender
    Commented Mar 17, 2016 at 10:05
10
$\begingroup$

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]

Mathematica graphics

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

MeshCoordinates @ breg
Improve this answer
$\endgroup$
9
  • $\begingroup$ @MartinBüttner this would work: Last[SortBy[ConnectedMeshComponents[reg], RegionMeasure]] $\endgroup$
    – halmir
    Commented Mar 17, 2016 at 13:58
  • $\begingroup$ Thanks. in my laptop, the Mathematica 10.0 gives me this result $\endgroup$
    – xyz
    Commented Mar 17, 2016 at 14:45
  • $\begingroup$ @ShutaoTANG. That's strange. Is that the result you get from RegionBoundary @ alphaShapes2D[pts, 0.5]? $\endgroup$
    – RunnyKine
    Commented Mar 17, 2016 at 18:31
  • $\begingroup$ @RunnyKine Yes, today I ran the code in the other laptop that installed Mathematica 10.3, which still gives me this result $\endgroup$
    – xyz
    Commented Mar 18, 2016 at 1:27
  • 1
    $\begingroup$ @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]] $\endgroup$
    – RunnyKine
    Commented Mar 18, 2016 at 1:33
9
$\begingroup$

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]

enter image description here

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.

Improve this answer
$\endgroup$
0
4
$\begingroup$

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

enter image description here

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

enter image description here

{δ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};

enter image description here

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}

enter image description here

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}

enter image description here

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}

enter image description here

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ +1. It doesn't seem easy to extract the boundary from these curves though. $\endgroup$
    – RunnyKine
    Commented Apr 22, 2016 at 9:18
  • $\begingroup$ @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. $\endgroup$
    – xyz
    Commented Apr 22, 2016 at 9:23
  • $\begingroup$ Ok, that sounds like a nice approach. $\endgroup$
    – RunnyKine
    Commented Apr 22, 2016 at 9:28
  • $\begingroup$ @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. $\endgroup$
    – xyz
    Commented Apr 22, 2016 at 9:31
  • $\begingroup$ Yes, that approach is quite fast. $\endgroup$
    – RunnyKine
    Commented Apr 22, 2016 at 9:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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