How to compute the boundary of a group of ellipse efficiently?
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.
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] &
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.
