Find all partitions of a plane created by a set of lines
I have a set of lines (linear functions of the form y = mx + b) (120 of them!) and if I graph them all, then they partition the R^2 plane. The lines do not necessarily go through the origin.
What is the most efficient way to find all partitions created by a set of such lines? Personally, I'm having a hard time coming up with any way at all, let alone an efficient one. To be more clear, I include the following image of just 4 lines:
An example of a partition would be the set {(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7} , which is the partition created by the red, yellow and green lines in the first quadrant. Another example would be {(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2}, which is the triangle in the first quadrant that is bounded by the blue, red and green lines.
An example of a non partition would be {(x,y)|5x+3 <= y <= -30x+28}, which is the set bounded by the green line above and the blue line below. This is not a partition because there are several partitions contained within it (like the second set above, for instance), or overlapping it. The set {(x,y)|5x+3 <= y <= -30x+28 && 90x+7 <= y}, however, would be a partition.
The desired output would be a complete list of such sets:
{(x,y)|y <= -30x+28 && 5x+3 <= y <= 60x+2},{(x,y)| -30x+28<= y && 60x+2 <= y <= 90x+7}... etc. They don't have to be given in this notation, of course.
I'm not sure how to approach this problem, and so, unfortunately, can't provide much in the way of what I have tried. Ideally, I would like to do this in R, Python, Mathematica or MATLAB, but I'm open to any option at this point.
EDIT: Since there seems to be an issue with notation, I'll clarify a bit. It is sufficient to simply get a list of conditions on points, such that all points meeting that condition would exactly define a partition. For example, a long list of intersections, would be fine: y <= 5x+3 && y >= 90x+7 && y<= -30x+28 is perfectly good output defining a partition. The desired output, of course, is a complete list of such partitions (as defined above).
Here is a solution in Mathematica. The method involves finding the line intersection points, the line fragments and then the partitions, while keeping track of which lines the points are connected to.
y[1] = 5 x + 3;
y[2] = 90 x + 7;
y[3] = -30 x + 28;
y[4] = 60 x + 2;
findpoints[n_] := Module[{},
xp = DeleteCases[{{First[#1], First[#2]},
Solve[Last[#1] == Last[#2], x]} & @@@
DeleteCases[
Tuples[Array[{#, y[#]} &, n], 2],
{{x_, _}, {x_, _}}], {_, {}}];
yp = y[#[[1, 1]]] /. #[[2, 1]] & /@ xp;
MapThread[{#1, {#2, #3}} &,
{xp[[All, 1]], xp[[All, 2, 1, 1, 2]], yp}]]
xyp = findpoints[4];
{xmin, xmax} = Through[{Min, Max}@
xyp[[All, 2, 1]]] + {-0.7, 0.7};
outers = Flatten[Array[Function[n,
MapThread[List[{{n, 0}, {##}}] &,
{{xmin, xmax}, y[n] /.
List /@ Thread[x -> {xmin, xmax}]}]], 4], 2];
xyp = Join[outers, xyp];
findlines[p_] := Module[{},
pt = DeleteCases[
Cases[xyp, {{p[[1, 1]], _}, _}], p];
n = First@Cases[pt,
{_, First@Nearest[Last /@ pt, Last[p]]}];
{{First[p], First[n]}, {Last[p], Last[n]}}]
lines = Map[findlines, xyp];
(* boundary lines *)
{ymin, ymax} = Through[{Min, Max}@outers[[All, 2, 2]]];
{lbtm, rbtm, ltop, rtop} = {{xmin, ymin},
{xmax, ymin}, {xmin, ymax}, {xmax, ymax}};
xminlines = Partition[Union@Join[{ymin, ymax},
Cases[xyp, {_, {xmin, _}}][[All, 2, 2]]], 2, 1] /.
x_Real :> {xmin, x};
xmaxlines = Partition[Union@Join[{ymin, ymax},
Cases[xyp, {_, {xmax, _}}][[All, 2, 2]]], 2, 1] /.
x_Real :> {xmax, x};
lines2 = Join[Last /@ lines, xminlines, xmaxlines,
{{lbtm, rbtm}}, {{ltop, rtop}}];
ListLinePlot[lines2]
(* add vertex points *)
xyp2 = Join[xyp, {
{{SortBy[Cases[outers, {_, {xmin, _}}],
Last][[-1, 1, 1]], -1}, ltop},
{{SortBy[Cases[outers, {_, {xmax, _}}],
Last][[-1, 1, 1]], -1}, rtop},
{{SortBy[Cases[outers, {_, {xmin, _}}],
Last][[1, 1, 1]], -1}, lbtm},
{{SortBy[Cases[outers, {_, {xmax, _}}],
Last][[1, 1, 1]], -1}, rbtm}}];
anglecalc[u_, v_] := Mod[(ArcTan @@ u) - (ArcTan @@ v), 2 π]
getlineangles[] := Module[{},
(* find the angles from current line
to all the linked lines *)
angle = Map[
anglecalc[{c, d} - {g, h}, # - {g, h}] &,
union = DeleteCases[Union@Join[
Last /@ Cases[lines2, {{g, h}, _}],
First /@ Cases[lines2, {_, {g, h}}]],
{c, d}]];
Sort[Transpose[{N@angle, union}]]]
getpolygon[pt_, dir_] := Module[{},
Clear[p];
p[n = 1] = {{a, b}, {c, d}} = pt;
(* find the angles from vector (0, -1) or (0, 1)
to all the linked lines *)
angle = Map[anglecalc[If[dir == 1, {0, -1}, {0, 1}], # - {c, d}] &,
union = Union@Join[
Last /@ Cases[lines2, {{c, d}, _}],
First /@ Cases[lines2, {_, {c, d}}]]];
lineangles = Sort[Transpose[{N@angle, union}]];
(* next point *)
p[++n] = {{e, f}, {g, h}} = First@
Cases[xyp2, {_, lineangles[[1, 2]]}];
While[Last[p[n]] != Last[p[1]],
lineangles = getlineangles[];
(* reset first point *)
{{a, b}, {c, d}} = {{e, f}, {g, h}};
(* next point *)
p[++n] = {{e, f}, {g, h}} = First@
Cases[xyp2, {_, lineangles[[1, 2]]}]];
Array[p, n]]
len = Length[xyp];
polygons = Join[Array[(poly[#] = getpolygon[xyp[[#]], 1]) &, len],
Array[(poly[# + len] = getpolygon[xyp[[#]], 2]) &, len]];
graphics = DeleteDuplicates /@ Array[Last /@ poly[#] &, 2 len];
sortedgraphics = Sort /@ graphics;
positions = Map[Position[sortedgraphics, #] &,
DeleteDuplicates[sortedgraphics]][[All, 1, 1]];
unique = poly /@ positions;
poly2 = unique[[All, All, 2]];
poly2 = Delete[poly2,
Array[If[Length[Intersection[poly2[[#]],
Last /@ Take[xyp2, -4]]] == 4, {#}, Nothing] &,
Length[poly2]]];
len2 = Length[poly2];
poly3 = Polygon /@ Rest /@ poly2;
Array[(centroid[#] = RegionCentroid[poly3[[#]]]) &, len2];
Show[Graphics[Array[{ColorData[24][#],
poly3[[#]]} &, len2], AspectRatio -> 1/GoldenRatio],
Graphics[Array[Text[#, centroid[#]] &, len2]]]
unique2 = Extract[unique,
Flatten[Position[unique[[All, All, 2]], #] & /@ poly2, 1]];
makerelations[oneconnection_, areanumber_] := Module[{i},
If[Intersection @@ oneconnection == {} ||
(i = First[Intersection @@ oneconnection]) < 1,
Nothing,
centroidx = First[centroid[areanumber]];
linepos = y[i] /. x -> centroidx;
relation = If[linepos < Last[centroid[areanumber]],
" >= ", " < "];
string = StringJoin["y", relation, ToString[y[i]]]]]
findrelations[n_] := Module[{},
areanumber = n;
onearea = unique2[[areanumber]];
connections = Partition[First /@ onearea, 2, 1];
strings = DeleteDuplicates@
Map[makerelations[#, areanumber] &, connections];
StringJoin["Area ", ToString[areanumber],
If[areanumber > 9, ": ", ": "],
StringRiffle[strings, " &&\n "]]]
Show[Plot[Evaluate@Array[y, 4], {x, -1, 1.5},
PlotLegends -> "Expressions", Axes -> None],
Graphics[Array[Text[#, centroid[#]] &, len2]]]
Column@Array[findrelations, len2]
Area 1: y >= 28 - 30 x && y < 3 + 5 x Area 2: y >= 2 + 60 x && y >= 28 - 30 x && y < 7 + 90 x Area 3: y < 28 - 30 x && y < 7 + 90 x && y < 2 + 60 x && y < 3 + 5 x Area 4: y >= 3 + 5 x && y >= 28 - 30 x && y < 2 + 60 x Area 5: y < 3 + 5 x && y >= 2 + 60 x && y < 7 + 90 x Area 6: y < 28 - 30 x && y >= 2 + 60 x && y >= 3 + 5 x && y < 7 + 90 x Area 7: y < 28 - 30 x && y >= 3 + 5 x && y < 2 + 60 x Area 8: y < 28 - 30 x && y >= 7 + 90 x && y >= 3 + 5 x Area 9: y < 2 + 60 x && y >= 7 + 90 x Area 10: y >= 28 - 30 x && y >= 7 + 90 x Area 11: y < 3 + 5 x && y >= 7 + 90 x && y >= 2 + 60 x