# Maximize area of intersecting triangle and square with constraints

There is this nice problem on mathoverflow to find a maximum intersection area of a square and a triangle with a constraint that they both have unit areas.

I tried to solve it with a help of Mathematica as follows

iarea[a_, rot_, d_] := Module[{sqr, trg, pA, h, pD, pC},
(*define a unit square*)
sqr = Polygon[{{0, 0}, {0, 1}, {1, 1}, {1, 0}}];
(*assuming the 1st triangle vertex is at origin,
and the 2nd in the direction rot and distance a *)
pA = a {Cos[rot], Sin[rot]};
(*move distance h in the perpendicular direction to find the 3rd vertex*)
h = 2/a;
pD = d {Cos[rot], Sin[rot]};
pC = pD + h {Cos[rot + \[Pi]/2], Sin[rot + \[Pi]/2]};
trg = Polygon[{{0, 0}, pA, pC}];
(* compute intersection and its area*)
Area[RegionIntersection[{sqr, trg}]]
]


I am trying to beat this value

iarea[Sqrt[2], 0.0, 0]
(*0.828427*)


with the help of NMaximise as follows

NMaximize[{iarea[xa, xr, xd],xa < 3 && xd <= xa && xa > 0 && xd > 0 && xr > 0 && xr < \[Pi]/2}, {xa, xr, xd}]


where I put some constraints on the angle rot, the length of one side a, and the d. However, the run takes too long, I was not patient enough to wait. Can I hope for a solution in this way? Is there a better/faster way?

• iarea[Sqrt[2], 0.0, 0] returns RegionIntersection::reg: {Polygon[{{0,0},{0,1},{1,1},{1,0}}],Polygon[{{0,0},{1.41421,0.},{8.65956*10^-17,1.41421}}]} is not a correctly specified region. and the input for me in version 11.2.0.0. Sep 28, 2017 at 15:56
• @user64494 Strange, I am on MA 11.0.1.0 and on a fresh kernel there is no problem Sep 28, 2017 at 16:30
• Area[RegionIntersection[sqr, trg]] worked both version 11.1 and 11.2. Sep 29, 2017 at 6:59
• In 10.4 Area[RegionIntersection[sqr, trg]] gets it right about 50% of the time. For the rest is gets the intersection wrong -- by missing out vertices or getting their order in Polygon. Sep 29, 2017 at 11:52
• @aardvark2012 is it a bug in MA10.4 ? Sep 29, 2017 at 12:05

I don't think you can do better that 2(Sqrt[2] - 1), but here's a way to get NMaximize working on the problem. I'm not sure exactly what it is about your approach that makes it take so long (I ran out of patience, too), but it could be that RegionIntersection is causing trouble. You could try plotting

Plot3D[iarea[a, 0, d], {a, d} \[Element] ImplicitRegion[0 < d < a < 3, {a, d}]]


and see if it looks anything like mine:

It's no surprise that a surface like that takes a long time to navigate! The reason for this (in my case) is that RegionIntersection[square, triangle] does things like this:

RegionIntersection[Rectangle[{0, 0}, {1, 1}], Triangle[{{0, 0}, {1, 0}, {0.5, 2}}]]

(* Polygon[{{0, 0}, {0, 0}, {0, 0}, {0.25, 1}, {0.75, 1.}}] *)


Clearly the Polygon of the intersection is a little messed up, and is missing the point {1, 0}. (Using exact numbers, ie 0.5 -> 1/2 in triangle, removes the problem in this case). So I've had to include my own intersection function. Pretty much all the code in this answer is constructing a working RegionIntersection function, so let me get that out of the way now.

Region Intersections

The first function, lineintersection, tests to see if the lines between two pairs of vertices intersect, and if so, where. The second, triangleintersection calculates all the vertices of the intersection and returns the intersection Polygon. Note that the arguments for this function are set up specifically for this maximization problem, and certainly can't be used in a more general setting without substantial rewriting. Specifically, it calculates intersection of the square Rectangle[{0, 0}, {1, 1}] with a triangle with vertices {{x, 0}, {y, 0}, {z, 2/(y - x)}} (I'll explain why those vertices below):

lineintersection[{x : {x1_, x2_}, y : {y1_, y2_}}, {u : {u1_, u2_}, v : {v1_, v2_}}] :=
Module[{denom = (u2 - v2) (x1 - y1) - (u1 - v1) (x2 - y2)},
If[denom == 0, Nothing,
With[{a = (-v2 x1 + u2 (-v1 + x1) + u1 (v2 - x2) + v1 x2)/denom,
b = -(-x2 y1 + u2 (-x1 + y1) + u1 (x2 - y2) + x1 y2)/denom},
If[0 <= a <= 1 && 0 <= b <= 1,
x + a (y - x)(* or u+b(v-u)*), Nothing]
]
]
]

triangleintersection[x1_?NumericQ, y1_?NumericQ, z1_?NumericQ] :=
Module[{x = {x1, 0}, y = {y1, 0}, z = {z1, 2/(y1 - x1)}, allpoints},
allpoints = DeleteDuplicates[Join[
(* Find all points where a square edge intersects a triangle edge *)
lineintersection @@@ Tuples[{Subsets[{x, y, z}, {2}],
{{{0, 0}, {0, 1}}, {{0, 1}, {1, 1}}, {{1, 1}, {1, 0}}, {{1, 0}, {0, 0}}}}],
(* Include vertices of the square that lie inside the triangle, and vice versa *)
Join[
Pick[#, RegionMember[Rectangle[{0, 0}, {1, 1}], #]] &@{x, y, z},
Pick[#, RegionMember[Triangle[{x, y, z}], #]] &@{{0, 0}, {0, 1}, {1, 1}, {1, 0}}]],
Norm[#1 - #2] < 10^-10 &];
Polygon[SortBy[allpoints, ArcTan @@ (# - Mean[allpoints]) &]]
]


Moving on to the actual question...

Problem setup

I've taken the liberty of changing the space of triangles we're searching over.Given a square with vertices {{0, 0}, {0, 1}, {1, 1}, {1, 0}}, I'm going to assume (for the moment) that:

The triangle has one (or possibly two) edge(s) lying along an edge of the square.

Without loss of generality, then, we can assume that the first two vertices are {x, 0} and {y, 0}, with y > x. Then the third vertex is {z, 2/(y - x)}. So it's still a 3D space to search over, as in your iarea, but contains the solution in the thread you linked to plus other maximal solutions, too.

You can explore this space by hand with

Manipulate[
intersection = triangleintersection[x, y, z];
Show[Graphics[{EdgeForm[Black], FaceForm[RGBColor[0.32, 0.5, 0.68]], Rectangle[{0, 0}, {1, 1}],
EdgeForm[Black], FaceForm[RGBColor[0.76, 0.73, 0.87]], Triangle[{{x, 0}, {y, 0}, {z, 2/(y - x)}}],
EdgeForm[{Thick, Black}], FaceForm[RGBColor[0.94, 0.28, 0.18]], intersection,
PointSize[Large], Point[intersection[[1]]]}],
PlotLabel -> "Area = " <> ToString@Area@intersection],
{{x, 1/2 - Sqrt[2]/2}, -1, 0.5}, {{y, 1/2 + Sqrt[2]/2}, 0.5, 2}, {{z, 0.5}, 0, 1}]


As a comparison to the surface I plotted above, a portion of the surface to be searched with a working intersection function looks like this:

Plot3D[Area@triangleintersection[x, y, 0.3], {x, -1, 0.4}, {y, 0.6, 2}]


After all this, finding the maximum area of the intersection is pretty straightforward. Define an area function using triangleintersection:

area[x1_?NumericQ, y1_?NumericQ, z1_?NumericQ] :=
With[{int = triangleintersection[x1, y1, z1]},
If[Length[int[[1]]] < 3, 0, Area@int]
]


and then maximize:

AbsoluteTiming[
sol = NMaximize[{area[x, y, z],
-1. <= x < 0.25 && 0.75 <= y <= 2. && 0 <= z <= 1},
{x, y, z}]
]

(* {2.21181, {0.828427, {x -> -0.233978, y -> 1.18024, z -> 0.564872}}} *)


Which is fairly speedy and successfully finds a triangle with an intersection area of 2 (Sqrt[2] - 1).

triangle = Triangle[{{x, 0}, {y, 0}, {z, 2/(y - x)}} /. sol[[2]]];
intersection = triangleintersection[x, y, z] /. sol[[2]];
Show[Graphics[{EdgeForm[Black], FaceForm[RGBColor[0.32, 0.5, 0.68]], Rectangle[{0, 0}, {1, 1}],
EdgeForm[Black], FaceForm[RGBColor[0.76, 0.73, 0.87]], triangle,
EdgeForm[{Thick, Black}], FaceForm[RGBColor[0.94, 0.28, 0.18]], intersection,
PointSize[Large], Point[intersection[[1]]]}],
PlotLabel -> "Area = " <> ToString@Area@intersection]


Generalising to 5D

We can generalising to the full 5D subspace of triangles with unit with a straightforward adaptation of the area function from above. I'm going to add two more variables: a "displacement", d, which moves the xy-line off the x-axis, and a rotation, \[Theta], about the origin. Specifically, the vertices of the triangle are defined as

$$x = R_\theta.(x_1, d), \quad y = R_\theta.(y_1, d) \quad \text{and} \quad z = R_\theta.\left(z_1, \frac{2}{y_1 - x_1} + d\right)$$

The only differences with the 3D functions I used earlier are in the function arguments and the Module definitions. So to save screen ink, just replace the function and Module lines in the definition of triangleintersection with:

triangleintersection5[x1_?NumericQ, y1_?NumericQ, z1_?NumericQ,
d_?NumericQ, \[Theta]_?NumericQ] := Module[{
x = RotationMatrix[\[Theta]].{x1, d},
y = RotationMatrix[\[Theta]].{y1, d},
z = RotationMatrix[\[Theta]].{z1, 2/(y1 - x1) + d},
allpoints},


and define

area5[x1_?NumericQ, y1_?NumericQ, z1_?NumericQ, d_?NumericQ, \[Theta]_?NumericQ] :=
With[{int = triangleintersection5[x1, y1, z1, d, \[Theta]]},
If[Length[int[[1]]] < 3, 0, Area@int]
]


As before, you can explore this space with:

Manipulate[
triangle = Triangle[RotationMatrix[\[Theta]].# & /@ {{x, d}, {y, d}, {z, 2/(y - x) + d}}];
intersection = triangleintersection5[x, y, z, d, \[Theta]];
Show[Graphics[{EdgeForm[Black], FaceForm[RGBColor[0.32, 0.5, 0.68]], Rectangle[{0, 0}, {1, 1}],
EdgeForm[Black], FaceForm[RGBColor[0.76, 0.73, 0.87]], triangle,
EdgeForm[{Thick, Black}], FaceForm[RGBColor[0.94, 0.28, 0.18]], inersection,
PointSize[Medium], Point[intersection[[1]]]}],
PlotLabel -> ToString@Area@intersection],
{{x, 1/2 - Sqrt[2]/2}, -1, 0.5}, {{y, 1/2 + Sqrt[2]/2}, 0.5, 2},
{{z, 0.5}, 0, 1}, {{d, -0.2}, -1, 1}, {{\[Theta], Pi/4}, -Pi/2, Pi/2}]


Finding a maximal triangle:

AbsoluteTiming[
sol = NMaximize[{area5[x, y, z, d, \[Theta]], -1. <= x < 0.25 &&
0.75 <= y <= 2 && 0 <= z <= 1 && -1 <= d <= 1 &&
-Pi/2 <= \[Theta] <= Pi/2}, {x, y, z, d, \[Theta]}]
]

(* {10.0285, {0.828427, {x -> -0.00236893, y -> 1.41184, z -> 0.00572204,
d -> -5.14455*10^-7, \[Theta] -> 1.23502*10^-6}}} *)

triangle = Triangle[{RotationMatrix[\[Theta]].{x, d},
RotationMatrix[\[Theta]].{y, d},
RotationMatrix[\[Theta]].{z, 2/(y - x) + d}} /. sol[[2]]];
intersection = triangleintersection5[x, y, z, d, \[Theta]] /. sol[[2]];
Show[Graphics[{EdgeForm[Black], FaceForm[RGBColor[0.32, 0.5, 0.68]], Rectangle[{0, 0}, {1, 1}],
EdgeForm[Black], FaceForm[RGBColor[0.76, 0.73, 0.87]], triangle,
EdgeForm[{Thick, Black}], FaceForm[RGBColor[0.94, 0.28, 0.18]], intersection,
PointSize[Medium], Point[intersection[[1]]]}],
PlotLabel -> "Area = " <> ToString@Area@intersection]


Which gives the same area as before, and was in fact in the earlier search space.

• So with your approach you can find a family of solutions with the same maximal area. As TimothyBudd writes on mathoverflow: "If you move the tip as well as the base in such a way that the intersections with the vertical sides are unchanged, then the overlapping area does not change." Interesting! Oct 3, 2017 at 7:42