Binary boolean operation between multi-island polygons A and B specified by contour
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Def: input polygons are A and B for binary ops and A for unary ops
Def: canon: an unary operator on a single input polygon that resolves the
     self intersections and irregularities of that polygon
Def: segment is a line or an arc
Def: curve is an ordered list of segments with no intersection
Def: face: a closed loop of curves (represented by an ordered list of
     curve IDs), without any curve crossing it; may embed other faces but
     the curves of those faces may not touch the curves of the outer face
Def: output tree [of faces]: a tree of faces; each node is a face with
     child nodes that are smaller faces fully within

1. switch to topology
	1.1. import all segments (from both A and B) into a 2D space, remembering
	     which segment is referred by A or B; compute intersections of any two
	     segments
	1.2. split segments at intersections; NOTE: intersection coords should be
	     integer which means spliting a segment in the middle can change the
	     slope of the two new segments and introduce new interscetions.
	1.3. merge fully overlapping segments, remembering how many times A and B
	     segments participated in the overlap (same count field as in 1.1)
	1.4. collect segments into curves between intersections

2. map faces
	- Algo: https://cp-algorithms.com/geometry/planar.html
		TODO: the actual algorithm is slightly different, document from scratch
	- Note: overlapping lines, e.g. in a bone, O--O case may result in
	        having a curve that's not part of any face; prune such curves
	- Note: for curves of faces, each curve should have a list of
	        faces that use the curve (this is at most two faces per curve)
	- Note: calculate the approximate area of each face to throw away the outer
	        one that overlaps all
	- Note: throw away curves that do not have at least one face associated

3. mark input polarity: for each face
	3.1. choose an internal point of the face
	3.2. face->inA := 1 or 0 depending on if the point is inside input polygon A
	3.3. face->inB := 1 or 0 depending on if the point is inside input polygon B
	3.4. compute polarity of faces: for each face compute face->out from
       face->inA and face->inB, depending on the operator, as:
		- A union B:     out = inA or inB
		- A sub B:       out = inA and not inB
		- A intersect B: out = inA and inB
		- A xor B:       out = inA xor inB
		- canon A:       out = inA

	Note: in point 3.2. anmd 3.3. tests the search needs to use the new geometry
	      created in step 1, not the original input vertices. The difference is
	      that the new model may have some extra vertices that are integer-rounded
	      changing the slope of some edges slightly so the "internal point"
	      chosen using this model may not be an internal point in the original
	      input.

	Note: for point 3.2. anmd 3.3. the algorithm is specified in
	      point-in-poly.html; the direction vector is taken to point inside the
	      triangle at the bottom of the right-most face corners.

4. figure which curves to keep: for each curve:
	- if the curve has two adjacent faces, mark the curve "pruned" if the
	  ->out of those two faces are equal
	- if the curve has only one adjacent face, mark the curve "pruned" if
	  face's ->out is 0 (not filled)

5. gather curves into output faces:

	Naive: Re-run step 2 but ignore curves marked "pruned" then re-run step 3.

	Optimized: remove faces that are either side of "pruned" curves and re-run
	           step 3 only for these

	This step can be skipped if there was no pruning in step 4 (and the
	face list from step 2 can be used as output faces).

	The result is an unordered "list of output faces".

6. stack output faces, building the output tree of faces:
	6.0. create an output tree with the root node being an infinitely large
	     negative face
	6.1. take the largest remaining output face "off the list of output faces"
	6.2. check whether the output face is within any of the leaves
		- if yes: insert the output face as a child face of that leaf face,
		  with polarity being the inverse of the leaf face
		- if not: when face->out==1 the output face is a new root in the
		  output tree with positive polarity otherwise discard it
		Note: since there are no intersections, this is real cheap: start a ray
		      from any point of the face and pick the first face it hits in odd
		      number of times
	6.3. go back to 7.1. if "list of output faces" is not empty

	The result is the output tree.

	(This is to figure the polarity of each output face)

7. flatten the tree:
	Iterate over each face in the output tree; for positive faces create
	a polygon island; for negative faces create a polygon hole within the
	island previously created for the parent face.

8. cleanup
	8.1. optional: visit each output contour and look at each interesected 
	     vertex; if the two segments on the two sides can be converted into
	     a single segment (e.g. coaxial lines or arcs on the same circle), merge
	     them to reduce the number of vertices
	8.2. discard and free the output tree (created in 7.)
	8.3. discard and free the list of output face (created in 6., should be empty by now)
	8.4. free all curves and faces used