How to compute the boundary of a group of ellipse efficiently?

Question

Update1

I discovered that the bountary of a series of ellipses consists of the following three parts

Part (I): the first ellipse's effective black-segment;

Part (II): the effective envelope-points of the ellipes from the second to last second;

Part (III): the last ellipse's effective red-segment;

Given that there are $n$ ellipses $E_1,E_2,\cdots,E_n$ on the plane.

• For the first ellipse $E_1$, the black segment that outside the ellipse $E_i(i=2,\cdots,n)$ is effective;

• For the envelope-point that on the ellipse $E_i(i=2,\cdots,n-1)$, the envelope-point $\theta_i$ that outsides the ellipse $E_j(j=1,\cdots,i-1,i+1,\cdots,n)$ is effective;

• For the last ellipse $E_n$, the red segment that outside the ellipse $E_i(i=1,\cdots,n-1)$ is effective;

Data

For a ellipse, which owns the following parametric formula

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

where, $\theta \in [0,2\pi]$

matThetaList =
 {{{{-0, -5, 0}, {-5.2203, 0, 1.7945}}, {2.4798, 5.7546}}, 
  {{{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}}, {3.1275, 6.2599}},
  {{{-1.8203, -4.7553, 0.2473}, {-5.6022 , 1.5451, 3.0486}}, {0.7316, 3.3481}},
  {{{-2.9427, -4.4550, 0.3147}, {-5.7755, 2.2700, 4.0578}}, {1.1944, 3.4426}}};

here, the variable matThetaList stores the ellipse $E_i$'s coefficient $\{\{a_i,b_i,c_i\},\{d_i,e_i,f_i\}\}$ and envelope-points $\theta_i^1,\theta_i^2$

Namely, matThetaList= $\{ \\ \{\{\{a_1,b_1,c_1\},\{d_1,e_1,f_1\}\},\{\theta_1^1,\theta_1^2\}\},\\ \{\{\{a_2,b_2,c_2\},\{d_2,e_2,f_2\}\},\{\theta_2^1,\theta_2^2\}\},\\ \cdots \\ \}$

---

I have implemented this in the Answer, However, owing to the function FindBoundary[] will be called many times, the performance of my function is very slow.

So I would like to know:

• Is there other more better/efficient algorithm to solve the boundary of the Ellipses $E_1,\cdots,E_n$?.

Update2

---

For the general case(all the sections are the complete ellipse), RunnyKine's solution works well and it was very fast. However, when the section was a partial ellipse, that solution failed. Here is a partial ellipse case

(*data for ellipse segments*)
(*About ellipsePoints[], please see my answer below*)
ellipseMat = 
{{{0.,-5.,0.},{-5.22027,0.,0.294118}},
 {{-0.418837,-4.98459,0.228686},{-5.32183,0.392295,-0.033668}},  
 {{-0.858274,-4.93844,0.325822},{-5.41893,0.782172,-0.364501}},
 {{-1.32336,-4.86185,0.291034},{-5.51219,1.16723,-0.688098}},
 {{-1.82027,-4.75528,0.123195},{-5.60223,1.54509,-0.994631}},
 {{-2.35676,-4.6194,-0.179982},{-5.68973,1.91342,-1.27478}},
 {{-2.94275,-4.45503,-0.622558},{-5.77547,2.26995,-1.5198}},
 {{-3.59125,-4.2632,-1.2113},{-5.86038,2.61249,-1.72161}},
 {{-4.31974,-4.04509,-1.95715},{-5.94562,2.93893,-1.87293}},
 {{-5.15241,-3.80203,-2.8775},{-6.0327,3.24724,-1.96744}},
 {{-6.12372,-3.53553,-4.00001},{-6.12372,3.53553,-2.00001}}};

ellipseDomain = 
 {{2.38622,7.03856},{2.49067,6.93411},{2.57819,6.84659},{2.65607,6.76871},
  {2.72819,6.69659},{2.79696,6.62782},{2.86409,6.56069},{2.93095,6.49383},   
  {2.99873,6.42605},{3.06856,6.35622},{3.1416,6.28318}};

Graphics[Line[Append[#, First@#]] & /@ 
  MapThread[ellipsePoints, {ellipseMat, ellipseDomain}]]

When I sampling more sections($300$), I discovered that the boundary should be as below:

Answer 1

Solution for Update1

The ellipses described by matThetaList can be plotted by

Show[ParametricPlot[#[[1]].{Sin[θ], Cos[θ], 1}, {θ, 0, 2 Pi}] & /@ matThetaList, 
    PlotRange -> All]

To describe each of these four curves as an ImplicitRegion, first eliminate θ from the parametric equations given in the question,

h = Total[#^2 & /@ ({Sin[θ], Cos[θ]} /. Solve[{x - c == a Sin[θ] + b Cos[θ], 
    y - f == d Sin[θ] + e Cos[θ]}, {Sin[θ], Cos[θ]}] // Flatten)]
(* (c d - a f - d x + a y)^2/(b d - a e)^2 + 
   (-c e + b f + e x - b y)^2/(b d - a e)^2 *)

then use h < 1 to define each ImplicitRegion, combine the four into a single region with RegionUnion, improve resolution with BoundaryDiscretizeRegion, and plot with RegionPlot.

RegionPlot[BoundaryDiscretizeRegion[RegionUnion[ImplicitRegion[
    (h /. Thread[{a, b, c, d, e, f} -> Flatten@#[[1]]]) < 1, {x, y}] & /@ 
    matThetaList], MaxCellMeasure -> .25], PlotStyle -> White]

which is the desired envelope of the four ellipses. The computation takes less than two seconds on my PC.

In the same way, the curve for matThetaList2, added in the OP's answer, can be computed quite accurately in three seconds.

Solution for Update2

The more difficult problem of finding the envelop of several truncated ellipses is treated similarly, the difference being that each ImplicitRegion is to be truncated. To do so, first compute the lines that truncate each ellipse. This is accomplished, first, by determining the points on the ellipses where truncation occurs, and next determining the line that connects these two points.

lims = Simplify[#[[1, 2]] + (#[[2, 2]] - #[[1, 2]])/(#[[2, 1]] - #[[1, 1]]) 
    (x - #[[1, 1]])] & /@  
    (# & /@ MapThread[{#1.{Sin[θ], Cos[θ], 1} /. θ -> #2[[1]], #1.{Sin[θ], Cos[θ], 1} /. 
    θ -> #2[[2]]} &, {ellipseMat, ellipseDomain}])
(* {-3.28469 + 1.06453*10^-6 x, -3.26025 - 0.0787005 x, -3.27957 - 0.158383 x, ... *)

Finally, this line is used in a second condition of ImplicitRegion.

RegionPlot[BoundaryDiscretizeRegion[RegionUnion[MapThread[ImplicitRegion[
    (y > #2) && ((h /. Thread[{a, b, c, d, e, f} -> Flatten@#1]) < 1), {x, y}] &, 
    {ellipseMat, lims}]], MaxCellMeasure -> .25], PlotStyle -> White]

as desired. This computation requires just under four seconds on my PC, measured with AbsoluteTiming. If, as seems likely, timing scales linearly with the number of ellipses, the 200-ellipse computation mentioned in the question would take about 70 seconds.

Addendum

To verify this timing estimate, and also to verify that no catastrophe occurs as the number of ellipses is increased, I generated 101 truncated ellipses by interpolation from the eleven ellipses given in the question.

interpDomain = Interpolation[#] & /@ Transpose[ellipseDomain];
tableDomain = Table[Through[interpDomain[i]], {i, 1, 11, 1/10}];
interpMat = Interpolation[#] & /@ Transpose[Flatten[#] & /@ ellipseMat];
tableMat = Table[Partition[Through[interpMat[i]], 3], {i, 1, 11, 1/10}];

Here is a plot of the truncated ellipses in Red and their truncating lines in Black.

Then, applying the same code used above provides the envelope.

in just under fifty seconds, again measuring runtime with AbsoluteTiming. My Windows PC has four processors, each with two threads. Interestingly, the computation used 50% or more of the total CPU capacity, indicating very effective parallelization for these geometric functions. A 301 ellipse calculation produced essentially the save envelope curve, but running just over three times as long. Strangely, a 201 ellipse calculation produced the same envelope curve, but with a small glitch in one location. I have not pursued why.

Answer 2

Here is an approach that seems to work for the cases presented, I think it may be general if you're only dealing with two curves. For more than two curves, you'll just need to extend the approach a little.

getSurf[pt1_, pt2_] := 
 Module[{gr1 = Graphics[Line[pt1 ~ Join ~ {pt1[[1]]}]], 
   gr2 = Graphics[Line[pt2 ~ Join ~ {pt2[[1]]}]], reg},
  reg = BoundaryDiscretizeGraphics /@ {gr1, gr2};
  DeleteCases[MeshCoordinates@RegionUnion@reg, 
        Alternatives @@ (MeshCoordinates@RegionIntersection@reg)]]

The function just takes the points of the two curves and returns the points of the trimmed boundary. For the first case:

GraphicsRow[ListPlot[getSurf[pts1, pts2], AspectRatio -> 1, Joined -> #] & /@ 
           {False, True}]

For the second case:

getSurf[pts3, pts4] // ListPlot[#, AspectRatio -> 1] &

Answer 3

Here, I applied the discrete strategy(sampling $400$ points in a period $2\pi$) to calculate the boundary of ellipses $E_1,\cdots,E_n$. The algorithm mainly consists of four steps as follows :

• Using the ellipses $E_2,\cdots,E_n$ to trim the black segment of the first ellipse $E_1$;

• Using the ellipses $E_1,\cdots,E_{n-1}$ to trim the red segment of the last ellipse $E_n$;

• Delete the envelope-points of ellipses $E_2,\cdots,E_{n-1}$ that insides the all other ellipses(apart from itself).

• With the help of built-in ListCurvePathPlot[] to find the boundary

---

Implementation

FindBoundary[matThetaList_, firstBlackDom_, lastRedDom_] :=
 Module[{n, trimPolygon, envelope, validEnvelope, firstBlack, 
   lastRed, validFirstBlack, validLastRed, temp},
  n = Length[matThetaList];
  trimPolygon = 
   ellipsePoints[#, {0., 2 π}] & /@ matThetaList[[All, 1]];
  (*compute the coordinates of the ellipse*)
  envelope =
   Transpose[#1.{Sin[#2], Cos[#2], ConstantArray[1, Length@#2]}] & @@@
     matThetaList;
  (*compute the black segment of the first ellipse and red-segment of last ellipse*)
  firstBlack = ellipsePoints[matThetaList[[1, 1]], firstBlackDom];
  lastRed = ellipsePoints[matThetaList[[-1, 1]], lastRedDom];
  (*extract the effective envelope-point*)
  validEnvelope =
   Flatten[
    Pick[
     envelope,
     MapIndexed[validEnvelopeQ[trimPolygon, #1, First@#2] &, 
     envelope], True], 1];
  (*trim the in-effctive cut-points*)
  temp = firstBlack;
  Do[
   temp = Select[temp, ! inPolyQ[trimPolygon[[i]], #] &],
   {i, 2, n}
  ];
  validFirstBlack = temp;
  (*trim the in-effctive cutting-points*)
  temp = lastRed;
  Do[
   temp = Select[temp, ! inPolyQ[trimPolygon[[i]], #] &],
   {i, 1, n - 1}
  ];
  validLastRed = temp;
  (*visulize the boundary*)
  Show[
   Graphics[{Red, PointSize[Medium], Point[validEnvelope]}],
   ListCurvePathPlot[
    Join[validFirstBlack, validLastRed, validEnvelope], 
    AspectRatio -> Automatic]
  ]
 ]

here, firstBlackDom and lastRedDom are the domain/interval list of the black segment of the first ellipse $E_1$ and red segment of last ellipse $E_n$, respectively.

Auxiliary function

ellipsePoints[mat_, dom : {a_?NumericQ, b_?NumericQ}] := 
  mat.{Sin[#], Cos[#], 1} & /@ Range[a, b, Pi/200.0]
ellipsePoints[mat_, dom : {{_?NumericQ, _?NumericQ} ..}] := 
  mat.{Sin[#], Cos[#], 1} & /@
   Join @@ (Range[#[[1]], #[[2]], Pi/200.0] & /@ dom)

validEnvelopeQ[trimPolygon_, pt : {x_?NumericQ, y_?NumericQ}, idx_] :=
 Module[{trimPoly},
  trimPoly = Delete[trimPolygon, idx];
  And @@ (! inPolyQ[#, pt] & /@ trimPoly)
 ]
validEnvelopeQ[trimPolygon_, pts : {{_?NumericQ, _?NumericQ} ..}, idx_] := 
  validEnvelopeQ[trimPolygon, #, idx] & /@ pts

inPolyQ1 = Compile[{{poly, _Real, 2}, {x, _Real}, {y, _Real}},
   Block[{Xi, Yi, Xip1, Yip1, u, v, w},
    {Xi, Yi} = Transpose@poly;
    Xip1 = RotateLeft@Xi;
    Yip1 = RotateLeft@Yi;
    u = UnitStep[y - Yi];
    v = RotateLeft@u;
    w = UnitStep[-((Xip1 - Xi) (y - Yi) - (x - Xi) (Yip1 - Yi))];
    Total[(u (1 - v) (1 - w) - (1 - u) v w)] != 0]
   ];

 inPolyQ[poly_, pt : {x_, y_}] := inPolyQ1[poly, x, y]

Here, the inPolyQ[] came from Simon Woods's answer

TEST

matThetaList =
  {{{{-0, -5, 0}, {-5.2203, 0, 1.7945}}, {2.4798, 5.7546}}, 
   {{{-0.8583, -4.9384, 0.1765}, {-5.4189, 0.7822, 2.3088}}, {3.1275, 6.2599}},
   {{{-1.8203, -4.7553, 0.2473}, {-5.6022 , 1.5451, 3.0486}}, {0.7316,3.3481}},
   {{{-2.9427, -4.4550, 0.3147}, {-5.7755, 2.2700, 4.0578}}, {1.1944, 3.4426}}};

firstBlackDom = {{5.7546, 2 π}, {0, 2.4798}};
lastRedDom = {{3.4426, 2 π}, {0, 1.1944}};

FindBoundary[matThetaList, firstBlackDom, lastRedDom]

Another test

matThetaList2 =
{{{{0., -5., 0.}, {-5.22027, 0., 1.79454}}, {1.41134, 5.00261}}, 
 {{{-0.418837, -4.98459, 0.375413}, {-5.32183, 0.392295, 1.83067}}, {1.35018, 5.13206}}, 
 {{{-0.858274, -4.93844, 0.679152}, {-5.41893, 0.782172, 1.86633}}, {1.22638, 5.26839}}, 
 {{{-1.32336, -4.86185, 0.914622}, {-5.51219, 1.16723, 1.90933}}, {1.00134, 5.41161}}, 
 {{{-1.82027, -4.75528, 1.08554}, {-5.60223, 1.54509, 1.96716}}, {0.678781, 5.56818}}, 
 {{{-2.35676, -4.6194, 1.19594}, {-5.68973, 1.91342, 2.047}}, {0.293121, 2.34024, 2.95781, 5.76617}}, 
 {{{-2.94275, -4.45503, 1.25017}, {-5.77547, 2.26995, 2.15563}}, {2.14864, 3.23955}},
 {{{-3.59125, -4.2632, 1.25283}, {-5.86038, 2.61249, 2.29949}}, {2.0766, 3.38628}}, 
 {{{-4.31974, -4.04509, 1.20875}, {-5.94562, 2.93893, 2.48455}}, {2.04293, 3.47937}}, 
 {{{-5.15241, -3.80203, 1.12285}, {-6.0327, 3.24724, 2.71637}}, {2.02746, 3.54081}}, 
 {{{-6.12372, -3.53553, 0.999988}, {-6.12372, 3.53553, 2.99999}}, {2.02227, 3.58014}}};

FindBoundary[
  matThetaList2, {{0, 1.411339}, {5.00261, 2 π}}, 
                 {{3.580142, 2 π}, {0, 2.022272}}]

---

Here is two screenshots that applied my method.