Inspired by the blog post Find a point in every region defined by a system of linear equations, I have come up with the following method to find a point inside each bounded region. The assumptions are:

• No two lines are parallel.
• No three lines are coincident.

Due to numerical instability, it seems, using floats, the coefficients of the equations of the lines are taken to be integers (they could also have been taken to be fractions, of course). Then the method goes like follows:

1. 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)}:
   

2. 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)};
   

3. 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;
   

4. 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;
   

5. 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);
   

6. 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);
   

7. 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 n: According to Israel the number of regions, bounded and unbounded, is given by 1/2*n^2 + 1/2*n + 1, so the number of bounded regions is given by 1/2*n^2 - 3/2*n + 1, subtracting 2*n. This should equal the number of polygons, as obtained above.