This comes from an old post over on the SE about Kobon Triangles, but I've added some new functionality and thought I might share it here.
I enjoyed working on this, I learned how to construct a MeshRegion
from a set of points, and to find polygons by using FindCycle
. First I will give the code and then explain,
kobonTriangle[k_] :=
Module[{r0, r1, r2, pts, ilns, lines, edges, vertices, triangles, mesh},
r0 := RandomReal[{-1, 1}];
r1 := RandomReal[{-1, 0}];
r2 := RandomReal[{0, 1}];
pts = Transpose[
{Array[{r0, r1} &, k - 1],
Array[{r0, r2} &, k - 1]
}];
ilns = InfiniteLine /@ pts~Join~{{{0, 0}, {1, 0}}};
lines = Flatten[
Partition[Sort@#, 2, 1] & /@ Table[
Flatten[List @@@ (RegionIntersection[
ilns[[n]], #] & /@ Delete[ilns, n]), 1],
{n, Length@ilns}], 1];
vertices = Flatten[lines, 1] // DeleteDuplicates;
edges = lines /. MapIndexed[#1 -> First@#2 &, vertices];
triangles = FindCycle[Graph[#1 <-> #2 & @@@ edges], {3}, All];
mesh = MeshRegion[
vertices, {Line /@ edges,
triangles /. {a_ <-> b_, b_ <-> c_, c_ <-> a_} :>
Polygon[{a, b, c, a}]}];
(* Uncomment this section to have the number of triangles reported on the image *)
(* Labeled[
mesh,
Row[{"Number of lines = ", k, ", Number of Triangles = ",
Length@triangles}]] *)
mesh
];
Here are a couple of examples,
{kobonTriangle[5], kobonTriangle[8]}
<img src="http://i.stack.imgur.com/MsaOA.png" alt="enter image description here"/>
In any iteration, chances are you won't find the optimal solution. For example, for 5 and 8 lines, there are solutions with 5 and 15 triangles, respectively, rather than the 3 and 9. But if you run the code enough times, you can often find a near-optimal solution. There are most definitely better algorithms to search for these, but I like this one. I let it run for an hour and got these results:
<img src="http://i.stack.imgur.com/TJmmc.png" alt="enter image description here"/>
How it works
I was inspired by Trevor Simonton's javascript code here. The idea is to generate k
random lines that intersect so as to get a decent number of triangles. To that end, we start with one line that is oriented horizontally, and then generate k-1
lines that cross this line.
Here is the code to do this,
r0 := RandomReal[{-1, 1}];
r1 := RandomReal[{-1, 0}];
r2 := RandomReal[{0, 1}];
pts = Transpose[
{Array[{r0, r1} &, k - 1],
Array[{r0, r2} &, k - 1]
}];
ilns = InfiniteLine /@ pts~Join~{{{0, 0}, {1, 0}}};
You can see the lines via,
Graphics[ilns]
<img src="http://i.stack.imgur.com/IbQHV.png" alt="enter image description here"/>
We need to zoom out to see all the intersections
Graphics[ilns, PlotRange -> {{-2, 2.0}, {-2, 2.0}}]
<img src="http://i.stack.imgur.com/tietE.png" alt="enter image description here"/>
Now I would like to cut off the lines after the intersection points to create a closed shape. First I will use RegionIntersection
to find all the intersection points. Then I create line segments between each intersection point, but first I sort the intersection points to make sure that we don't have any overlapping line segments.
lines = Flatten[
Partition[Sort@#, 2, 1] & /@ Table[
Flatten[List @@@ (RegionIntersection[
ilns[[n]], #] & /@ Delete[ilns, n]), 1],
{n, Length@ilns}], 1];
vertices = Flatten[lines, 1] // DeleteDuplicates;
Graphics[{Line /@ lines, {Red, PointSize[Medium], Point /@ vertices}}]
<img src="http://i.stack.imgur.com/iwGXI.png" alt="enter image description here"/>
So we have our basic shape, but how to find the triangles, and only the non-overlapping triangles? By making a Graph
that is isomorphic to the shape above, we can take advantage of the Graph
functions in Mathematica
edges = lines /. MapIndexed[#1 -> First@#2 &, ipts]
Graph[edges, VertexLabels -> "Name"]
(* {{1, 2}, {2, 3}, {3, 4}, {5, 3}, {3, 6}, {6, 7}, {8, 9}, {9,
5}, {5, 4}, {9, 10}, {10, 1}, {1, 7}, {8, 10}, {10, 2}, {2, 6}} *)
<img src="http://i.stack.imgur.com/eP38F.png" alt="enter image description here"/>
Now we can find the triangles easily enough, and only non-overlapping triangles will be found because we've cut the lines into non-overlapping segments already.
triangles = FindCycle[Graph[#1 <-> #2 & @@@ edges], {3}, All]
Length@triangles
(* {{8 <-> 9, 9 <-> 10, 10 <-> 8}, {2 <-> 3, 3 <-> 6,
6 <-> 2}, {1 <-> 10, 10 <-> 2, 2 <-> 1}, {3 <-> 4, 4 <-> 5, 5 <-> 3}} *)
(* 4 *)
Now we just wrap it all up into a MeshRegion
for display purposes,
Labeled[
MeshRegion[
vertices, {Line /@ edges,
triangles /. {a_ <-> b_, b_ <-> c_, c_ <-> a_} :>
Polygon[{a, b, c, a}]}],
Row[{"Number of lines = ", k, ", Number of Triangles = ",
Length@triangles}]]
<img src="http://i.stack.imgur.com/Jd4sb.png" alt="enter image description here"/>
So this code is perhaps not efficient - I imagine that FindCycles
and the routine to find the intersection both scale at or worse than
$\mathcal{O}(n^2)$ but
$n$ is small so that is no worry.
One slight problem
For every time we get a decent shape like
kobonTriangle[14]
<img src="http://i.stack.imgur.com/eLTS8.png" alt="enter image description here"/>
we will get 5 that look like this,
<img src="http://i.stack.imgur.com/thIfx.png" alt="enter image description here"/>
where some of the lines are so long as to make the shape hard to view. The best plan I could come up with for discriminating against the latter shapes is by comparing the area of the triangles to the total area of a square encasing all the points. The shape with the best relative triangle area wins. Here is a routine that goes through 2000 9-line shapes and displays the best result alongside the current result,
evaluate[mesh_] := {(Area /@ MeshPrimitives[mesh, 2] //
Total)/(Times @@
Abs@*Subtract @@@ (MinMax /@
Transpose[MeshPrimitives[mesh, 0][[All, 1]]])),
Length@MeshPrimitives[mesh, 2]}
bestN = 0;
bestA = 0;
bestN = 0;
bestA = 0;
Dynamic[{bestN, bestTriangle, currentTriangle}]
With[{nt = 9},
Do[
currentTriangle = kobonTriangle[nt];
{area, ntriangles} = evaluate@currentTriangle;
If[ntriangles >= bestN,
bestN = ntriangles;
If[area > bestA,
bestA = area;
bestTriangle = currentTriangle];
];
, {2000}]
]