Given a set of segments, I would like to compute the set of closed polygons inside the convex hull of the set of the end of those segments. The vertices of the polygons are the intersections of the segments. For example, if you draw the 6 lines restricted which equations are: $x=-1$, $x=0$, $x=1$, $y=-1$, $y=0$, $y=1$, I would like the algorithm to output the four unit squares around the origin.
You're trying to compute something called the arrangement of the line segments. The arrangement is a planar straight-line graph whose vertices are either endpoints or intersection points of line segments. The vertices partition the lines segments into smaller segments, which are the edges of the arrangement. Finally, the vertices and edges partition the plane into several components; these are the faces of the arrangement.
Suppose you are given $n$ line segments, and $k$ pairs of those line segments intersect. (You don't know $k$ in advance.) The intersection points can be computed in $O((n+k)\log n)$ time using the Bentley-Ottmann sweep-line algorithm. (The running time can be reduced to $O(n\log n + I)$ using a more complicated algorithm, where $I\le k$ is the number of intersection points, but the additional effort is probably not worth it.) The algorithm can be modified to report, for each intersection point, which segments intersect at that point. Then it is simply a matter of careful bookkeeping to construct a standard data structure that represents the arrangement.
A few caveats to keep in mind:
The faces of the arrangement are not necessarily simple polygons. If the union of the line segments is connected, each bounded face is a weakly simple polygon, but a single edge can be incident to the face no both sides. If the union of the segments is disconnected, then faces can have holes, which are either trees or weakly simple polygons.
Most textbook (and Wikipedia) descriptions of the Bentley-Ottmann algorithm assume that the input segments are in general position. In this context, "general position" means that at most two segments intersect at any point; no segment endpoint lies on another segment (in particular, the intersection of two segments is either empty or a single point); and no two endpoints have the same $x$-coordinate. In practice, however, this assumption does not hold. As usual in computational geometry, handling degenerate cases correctly and robustly is the most difficult part of implementing the algorithm. (What does your implementation do when the input contains 42 copies of the segment [(0,0),(1,1)], 17 copies of the segment [(0,0),(2,2)], and 9 copies of the zero-length segment [(1,1),(1,1)]?)
For small numbers of segments (at most a few dozen), it's probably faster to use a brute-force $O(n^2)$-time algorithm to compute the intersection points.