User Manual: hull algorithm

From FLUX

The core of FLUX is a hull algorithm that identifies the Voronoi cell of the origin, in a 2-D field of points. Each point represents a separate fluxel, projected into the perpendicular plane of the current fluxel. The Voronoi cell is the locus in the plane that is closer to the central point than to any other. That locus is considered to be the cross section of the volume into which the central fluxon's magnetic flux has expanded. The geometry of the locus is used to calculate the magnetic field strength and related quantities such as the magnetic pressure gradient.

Voronoi tesselation of the plane is used for many applications; for a nice description see the textbook by Preparata and Shamos. Finding the Voronoi cell of a single point in a collection is a variant of the convex hull problem. It can be accomplished through a (conceptually) simple algorithm and analytic geometry.

The neighbor algorithm accepts a collection of neighbor candidates, and returns a (generally much smaller) collection of neighbors and associated hull vertices. The hull vertices are a generalization of points: they can represent either the intersection of two lines, or the exact angles of departure for lines that do not intersect in the appropriate quadrant. Hull vertices with no associated point are called "open", and correspond to a range of angles at which the central Voronoi cell extends to infinite distance from the origin.

Voronoi cell construction: an algorithm for sorted points

The original hull algorithm in FLUX (hull_2d) required a collection of points pre-sorted by angle. The algorithm works by constructing the perpendicular bisector of the line segment connecting each candidate neighbor with the central point (at the origin). The perpendicular bisector divides the plane into points which are closer to the origin or closer to the candidate. The smallest hull that can be constructed of such bisectors around the origin is the origin's Voronoi cell.

The hull vertices are the intersections of the perpendicular bisectors of the line segments between the origin and each adjacent neighbor.

The sorted-point algorithm, which is described in Preparata and Shamos's textbook on computational geometry, relies on the fact that subsequent neighbors (in clockwise order) should have bisector intersections that are also in clockwise order: if C is clockwise of B, which is clockwise of A, then the BC vertex should also be clockwise of the AB vertex.

Each point is examined in succession, and the order of its vertex intersections with the following and subsequent current neighbors is compared. If the vertex intersections are in forward order (or one does not exist) then the new point is considered to be a neighbor and is added to the current neighbor list. The previous and following neighbors are then examined for validity and are deleted if necessary, before the next candidate is examined.

The sorted-point algorithm is serviceable and fast but requires very slow preprocessing: all of the candidate points (there are ~300 in a typical calculation within FLUX) must be sorted in order of projected angle in the perpendicular plane of the fluxel being examined, even though all but ~6-8 are ultimately discarded.

Voronoi cell construction: an algorithm for unsorted points

The current Voronoi cell calculation in FLUX accepts unsorted points. It is similar to the sorted-point algorithm, except that at each stage a direct linear search is performed through the existing neighbor-candidate list to find the appropriate position at which the new candidate should be compared. When a new candidate is inserted into the list, the list is linearly copied out one position in memory to make room for the new candidate. These operations are slow, but at any given time there should be no more than ~10 neighbor candidates, so linear scan and direct array insertion are feasible and actually save time when compared with the task of quicksorting 300 neighbor candidates for each hull calculation.

To further save time, the current algorithm uses a crude approximation to the arctangent to calculate projected angles of candidates: the candidate angles are calculated using an octagonal atan2, which yields a monotonically increasing angle compared to the actual arctangent (and hence yields correct sorting order) but is not accurate. The octagonal arctangent is about a factor of 50 cheaper computationally than the true arctan. The more expensive angle calculation is performed at the end of the calcualtion, on the actual neighbors and their hull points only.

Voronoi cell construction: an example

The figure at left illustrates the process of FLUX's unsorted hull algorithm, hull_2d_us(), operating on a central fluxel (at the origin; marked "O" and colored red) and six other neighbor candidate fluxels, projected into the plane (located around the origin; numbered and colored grey).

Step 1

The code examines candidate 1 by constructing the perpendicular bisector of the line segment joining 1 and O. There are no other neighbors in the current neighbor list, so the only hull vertex (H1) is open and encloses an entire half-plane. Candidate 1 is renamed neighbor 1 (indicated by being colored green and labeled N1).

Step 2

The code examines candidate 2 by constructing the perpendicular bisector of the line segment joining 2 and O. N1 is the first existing neighbor widdershins of point 2, so the bisector is compared with N1's bisector. N1 is also the first existing neighbor clockwise of point 2, so the bisector is compared in the other direction as well. On the widdershins (left) side, the the bisectors intersect (at H1) and on the clockwise (right) side they don't intersect at all. Point 2 is kept as a neighbor and labeled N2. The hull vertex clockwise of N2 is labeled H2 and is open.

Step 3

The code examines candidate 3 by constructing the perpendicular bisector of the line segment joining 3 and O. N2 is the first existing neighbor widdershins of point 3, and N1 is the first existing neighbor clockwise of point 3, so the bisector is compared with the bisectors of N2 and N1. On the widdershins (left) side, N2 and 3's bisectors intersect (at H2), and 3's bisector does not intersect anything else. Hence N2 remains a valid neighbor, and 3 remains a candidate neighbor. On the clockwise side 3's and H1's bisectors do not intersect, so H1 remains a valid neighbor and 3 is a genuine neighbor. 3 is kept and labeled N3, and the intersection and non-intersection of the two respective bisectors are labeled H2 and H3 (which is open).

Step 4

The code examines candidate 4 by constructing the perpendicular bisector of the line segment joining 4 and O. N2 is the first existing neighbor widdershins of point 4, so the O-4 bisector is compared to the O-N2 bisector. The O-N2 bisector and the O-4 bisector intersect widdershins of H1 even though 4 is clockwise of N2. Hence N2 is no longer a valid neighbor and is eliminated; and 4 remains a candidate neighbor. N3 is the first existing neighbor clockwise of point 4, so the O-4 bisector is compared to the O-N3 bisector. They intersect and the O-4 bisector doesn't intersect anything else on the clockwise side. The clockwise-side intersection is clockwise of the widdershins-side intersection, so 4 is a valid neighbor. It supersedes N2, and H1 and H2 are updated accordingly.

Step 5

The code examines candidate 5 by constructing the perpendicular bisector of the line segment joining 5 and O. N2 is the first existing neighbor widdershins of point 5, so the O-5 bisector is compared to the O-N2 bisector. They intersect, and the O-5 bisector intersects the O-N1 bisector widdershins of the intersection with the O-N2 bisector, so N2 remains a valid neighbor and 5 remains a neighbor candidate. On the clockwise side, N3 is the first existing neighbor. The O-5 bisector intersects the O-N3 bisector widdershins of the O-N2/O-5 intersection, so 5 is not' a valid neighbor (the O-N3/O-5 intersection should be clockwise of the O-N2/O-5 intersection because N3 is clockwise of N2). 5 is rejected.

Step 6

The code examines candidate 6 by constructing the perpendicular bisector of the line segment joining 6 and O. The next neighbor widdershins of 6 is N3. The bisector of O-6 doesn't intersect the bisector of O-N3 on the widdershins side of 6, so 6 is a valid neighbor candidate. The next neighbor clockwise of 6 is N1. The bisector of O-N1 intersects the O-6 bisector clockwise of the O-6/O-N2 intersection, so N1 is no longer a valid neighbor. The O-6 bisector intersects the O-N3 bisector clockwise of the O-N2 bisector, so N2 remains a valid neighbor. Point 6 replaces N1 as a neighbor, and H3 and H1 are updated accordingly.

Use in PDL

It is possible in PDL, to print out the position of the hull points of any point. These data are not stored in the structure, but are computed every time you call the subroutine.

Regular hull

example:

p $world->vertex(2434)->hull

[
 [ -0.72915395  -0.54688343  -0.19717383  -0.49663858            0   -1.9487283   -1.9487283]
 [  0.11791907  -0.92819837   0.31985308  -0.43095506            0  -0.93230859  -0.93230859]
 [  0.85319265  -0.37114575   0.53154265  0.055678526            0   0.10436832   0.10436832]
 [  0.74031682   0.54766535   0.20707544   0.49428311            0    1.1740689    1.1740689]
 [ -0.12237895   0.91696011  -0.32490744   0.42328384            0    2.2254588    2.2254588]
 [ -0.85738216   0.38857475   -0.5289488 -0.026929217            0   -3.0907258   -3.0907258]
 [ -0.85774149   0.38829047   -0.5366718 -0.043989481            0   -3.0598083   -3.0598083]
]

The columns are as follows: neighbor_x, neighbor_y, hull_x, hull_y, open?, a_l, a_r Note that all of the positions are in the frame rotated so that the fluxel is in the z-direction and centered at (0,0) a_r and a_l are the same unless the hull vertex is an open one. The x,y coordinates of an open hull vertex do not mean anything. This routine should work for any vertex that is not an end-vertex.

Photospheric hull

example:

p $world->vertex(-5782)->photohull;

[
 [ -20.14462  14.439259          0          1]
 [    386.25 -77.483333          0          0]
 [-9.8476744  10.538372          0          0]
 [-13.338338  9.2293732          0          0]
 [-13.514913   6.639604          0          0]
]

Columns are as follows: hull_x, hull_y, hull_z, open? The positions are all absolute, and are not rotated or translated. At this time, we do not return the neighbor vertices. As with hull, the x,y coordinates of an open hull vertex do not mean anything.

Use: only use ->photohull on a vertex that is a begin or end vertex that is located on the photosphere that is in photosphere1. It is your job to make sure that you input a valid vertex, the routine will not check for you.