Defects are a local region of edges and faces that require additional edges and faces to repair.
The code that adds these includes the apparently difficult task of adding an edge whose projection onto the surface of the sphere does not add any intersections with existing edges. The previous general purpose code does this by doing an extensive and expensive test against other existing edges, involving dot and cross products.
The new code speeds this code up in three ways, but they all exploit a geometrical property. When you project a straight line within the sphere onto the surface of the sphere, it creates a portion of a great circle, which I refer to here as a "great arc". A set of such edges project onto a set of great arcs. If those arcs are in part of one hemisphere, and they are projected onto a plane that is parallel to that hemisphere's equatorial plane, they project onto a set of line segments, and this set will have the same intersections as the great arcs have!
Based on this observation, the code proceeds as follows. It finds the center of the vertexs of the edges, and uses the line from the center of the sphere through this center as the hemisphere's axis. It creates two more axes to make a three orthogonal axis set. It creates a plane through these two more axes and projects the edges onto it - thus creating the set of line segments. The code often tries to add the edge VO-V1 and the edge V1-V0. Such duplicates are discarded.
This projection looks remarkably like a spider web with a hole in it.
It partitions the plane into a grid of cells, centered on the center of these ends of the line segments, and wide enough so that extra line segments can be added around the existing ones without going outside the grid. One last cell is made for all the vertexs that fall outside the grid, called the 'universal cell'.
Each line segment is noted in all the cells it might intersect. Since there are about 100 cells, each cell only gets a few of the thousand or so line segments in a defect.
The code is determining whether any of these great arcs intersect with a great arc of a trial edge that is not yet in the set. The trial edge is similarly projected into the cells, and all the line segments in all cells that it might intersect are tested against it. Again, since line segment might appear in more than one cell, duplicates are detected and not tested again.
This reduces the number of intersection tests dramatically. The old code could be used to do the test, and is for line segments where one goes beyond the grid into the universal cell. However for the lines that are projecting into the grid, there are two much faster ways of doing the test.
It the lines are not almost parallel, a piece of math described in the code turns each segment into (x,y) = Ap + B and (x,y) = Cq + D where p and q are between 0 and 1. These two linear equations are then solved for p and q, and if and only if both are between 0 and 1, the line segments intersect. This can be done remarkably efficiently, and with only one comparison and no divides!
If the lines are almost parallel, then the simultaneous equations are numerically unstable or even involve a multiply by zero. This does not correctly detect whether the line segments intersect. Instead a coordinate system is built using one of the lines as the x-axis, (0,0) being on one end of the line. Both lines are projected into this coordinate space and some simple compares and a few multiplies and adds are needed to decide whether the two intersect.
The code contains support for comparing each of the above approaches and the old approach to see whether they get the same answer. They do except for 4 of the 70M pairs that the fast mris_fix_topology test does. In these 4 cases two line segments pass so close to each other than numerical differences in the accuracy of the answer seems to be important.
It is possible to add further edges to the set without having to reproject the existing lines, and this is utilized to avoid unnecessary projections as the defect is fixed.To allow these lines to be outside the initial area but still within the grid and hence using the above approaches, the initial grid is sized larger than needed for the initial segments.