- Generate the equations (here
8, as originally posed in the above-mentioned blog post):
n := 8:
generator := rand(1000) - 500:
a := [seq(generator(),i=1..n)]:
b := [seq(generator(),i=1..n)]:
c := [seq(generator(),i=1..n)]:
eqs := {seq(a[i]*x + b[i]*y + c[i] = 0,i = 1..n)}:
- Find all intersection points:
dropNames := (s) -> [rhs(select(has,s,x)[]),rhs(select(has,s,y)[])]:
points := {seq(seq({i,j} =
dropNames(solve({eqs[i],eqs[j]},{x,y}))
,j = i+1..n),i = 1..n)};
- Find all nearest neighbour pairs:
pairs := {}:
for i1 in {$1..nops(points)} do
for i2 in {$1..nops(points)} minus {i1} do
# Check that the line segment joining the two points obey:
# 1. It is part of some line (cond1)
# 2. It is not crossed by any line (cond2)
x1,y1 := rhs(points[i1])[1],rhs(points[i1])[2];
x2,y2 := rhs(points[i2])[1],rhs(points[i2])[2];
l1 := [seq](eval(lhs(eqs[j]),{x = x1,y = y1}),j = 1..nops(eqs));
l2 := [seq](eval(lhs(eqs[j]),{x = x2,y = y2}),j = 1..nops(eqs));
cond1 := evalb(lhs(points[i1]) intersect lhs(points[i2]) <> {});
cond2 := not member(true,
{seq}(evalb(l1[j] * l2[j] < 0),j = 1..nops(eqs))
);
if cond1 and cond2 then pairs := pairs union {[i1,i2]} end if
end do
end do:
pairs;
- Find all 'convex' triples:
triples := {}:
for i1 in {$1..nops(pairs)} do
for i2 in {$1..nops(pairs)} minus {i1} do
v1 := rhs(points[pairs[i1][2]]) - rhs(points[pairs[i1][1]]);
v2 := rhs(points[pairs[i2][2]]) - rhs(points[pairs[i2][1]]);
if
# Check for common middle point, and 'convexity'
evalb(pairs[i1][2] = pairs[i2][1]) and
evalb(v1[1] * v2[2] - v1[2] * v2[1] > 0)
then triples := triples union {[pairs[i1][],pairs[i2][2]]}
end if
end do
end do:
triples;
- By iteration, find all convex paths:
evolvePath := proc(path)
local triple;
# Closed path
if member([path[-1],path[1]],pairs) then path
# Not closed path
else map(
x -> `if`(path[-2..-1] = x[1..2],[path[],x[3]],NULL),
triples
)[];
end if
end proc:
evolvedPaths := proc(paths)
local pathsEvolved;
pathsEvolved := map(evolvePath,paths);
`if`(pathsEvolved <> paths,
evolvedPaths(pathsEvolved),pathsEvolved
)
end proc:
convexPaths := evolvedPaths(triples);
- Find all convex polygons by first removing all non-closed convex paths, and, subsequently, converting all lists to sets to remove cyclically redundancy:
polygons := map(convert,map(
x -> `if`(member([x[-2],x[-1],x[1]],triples),x,NULL),convexPaths
),set);
- Find the arithmetic mean centers of the polygons:
centers := map(p -> `+`(map(i -> rhs(points[i]),p)[]) / nops(p),polygons);
A consistency check for any given 