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

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:

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]

We need to zoom out to see all the intersections

Graphics[ilns, PlotRange -> {{-2, 2.0}, {-2, 2.0}}]

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

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}} *)

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

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]

we will get 5 that look like this,

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