-
Notifications
You must be signed in to change notification settings - Fork 0
User_Manual:_Code:_geometry.c
(Return to Code or the User Manual)
geometry.c is a geometry manipulation package for FLUX. The core of FLUX is coordinate free and deals with vectors as such; the code in geometry.c uses Cartesian coordinates and implements many basic, familiar vector functions from norm and dot product, through complex algorithm such as the hull finder. Here's a list, sorted in rough order of complexity, of all the routines.
These functions return the answer directly in the return value.
NUM norm2_2d(NUM *x)
squared norm of a 2-vector
NUM norm_2d(NUM *x)
norm of a 2-vector (involves a square root)
NUM norm2_3d
squared norm of a 3-vector
NUM norm_3d
norm of a 3-vector (involves a square root)
NUM inner_2d(NUM *p0, NUM *p1)
Inner (dot) product of two 2-vectors
NUM inner_3d(NUM *p0, NUM *p1)
Inner (dot) product of two 3-vectors
NUM cross_2d(NUM *p0, NUM *p1)
Cross product of two 2-vectors (this is an antisymmetric two-tensor, ie a pseudoscalar).
int above_plane(POINT3D A, POINT3D B, POINT3D C, POINT3D X)
Given three noncolinear points A, B, and C, and a test point X, determine whether X is "above" the oriented plane defined by the triangle ABC. ("above" here means displaced in the same direction as the vector AB × AC).
int in_simplex(POINT3D A, POINT3D B, POINT3D C, POINT3D D, POINT3D X)
Determine whether a test point X is inside the simplex whose vertices are given by the noncoplanar points A, B, C, and D. If ABCD is not a simplex (the points are coplanar or worse), the result is undefined.
These functions generally stuff values into a pre-allocated first argument.
void cross_3d(NUM *out, NUM *p0, NUM *p1)
Cross product of two 3-vectors (this is an antisymmetric two-tensor, ie a pseudovector).
void scale_3d(NUM *out, NUM *a, NUM b)
Scalar multiplication of a 3-vector and a scalar. The output and input vectors may be the same 3-vector if desired.
void scale_2d(NUM *out, NUM *a, NUM b)
Scalar multiplication of a 2-vector and a scalar. The output and input vectors may be the same 2-vector if desired.
void sum_3d(NUM *out, NUM *a, NUM *b)
Sum of two 3-vectors. Either a or b may be duplicated into out if desired.
void diff_3d(NUM *out, NUM *a, NUM *b)
Sum of a and -b. Either a or b may be duplicated into out if desired.
void cp_3d(NUM *out, NUM *in)
Copies in into out.
Both 2x2 and 3x3 matrices are supported, as 4-arrays and 9-arrays respectively. The rows run fast, and the columns run slow. As with the vector functions above, output arrays must be pre-allocated.
void rotmat_2d(NUM *out, NUM alpha)
Returns a rotation matrix (4 elements) from an angle (alpha in radians).
void rotmat_2d_fr_slope(NUM *out, NUM dy, NUM dx)
Returns a rotation matrix (4 elements) from a slope rather than an angle. All four quadrants are supported, but dx and dy must not both be zero.
NUM det_2d(NUM *mat)
Determinant of a 2x2 matrix (4 elements), returned as a scalar.
NUM det_3d(NUM *mat)
Determinant of a 3x3 matrix (9 elements), returned as a scalar.
mat_mult_2d(NUM *out, NUM *a, NUM *b)
Returns the product of (4-element) matrices a and b, in a 4-element matrix. You cannot duplicate the inputs into the output, because cross terms would be overwritten.
void transpose_2x2(NUM *mat)
Transposes a 2x2 (4-element) matrix in place.
void transpose_3x3(NUM *mat)
Transposes a 3x3 (9-element) matrix in place.
mat_mult_3d(NUM *out, NUM *a, NUM *b)
Returns the product of (9-element) matrices a and b, in a 9-element matrix. You cannot duplicate the inputs into the output, because cross terms would be overwritten.
mat_vmult_3d(NUM *out, NUM *mat, NUM *v)
Multiplies a 3x3 matrix (9 elements) by a column 3-vector (3 elements) on the right, returning a column 3-vector (3 elements).
vec_mmult_3d(NUM *out, NUM *mat, NUM *v)
Multiplies a 3x3 matrix (9 elements) by a row 3-vector (3 elements) on the left, returning a row 3-vector (3 elements).
void projmatrix(NUM *out, NUM *x0, NUM *x1)
Given two points in 3-space, returns a projection matrix that rotates the directed line segment x0->x1 into the Z axis; this is useful for projecting 3-vectors down into the perpendicular plane of a line segment (such as a fluxel). The output array should have room for 9 elements.
void reflect(NUM *out, NUM *point, PLANE *plane)
Given a point and a plane (which contains an origin and a normal vector), reflect the point through the plane. This is used for generating image fluxels, e.g. for photospheric boundaries.
void points2plane(PLANE *plane, NUM *p0, NUM *p1, NUM *p2)
Given three noncolinear points, returns the plane that contains them, as a PLANE structure. If the points are colinear, then some of the cross products yield 0. The PLANE is still formed but may have NaN components or a zero-length normal vector.
int p_l_intersection(NUM *out, PLANE *plane, NUM *p0, NUM *p1)
Given a plane and a line (specified by two points that lie on it), return the location of the intersection between them. Returns a flag indicating whether the intersection is between the two points. If the line is in the plane, then a valid point is returned. If the line is parallel to, but outside, the plane, then one or more components of the output vector are NaN.
int xy_l_intersection(NUM *out, NUM *p0, NUM *p1)
Special case of p_l_intersection, above, for when the plane is the XY plane.
int p_inside_tri(NUM *tri0, NUM *tri1, NUM *tri2, NUM *p)
Given a (2-D) triangle and a (2-D) point, determine whether the point is strictly inside the triangle, or outside it. Return value is 1 if the point is in the triangle, 0 if it is on the boundary or outside.
int trivloop(FLUXON *f)
Given a fluxon, determines whether it is a "trivial" loop -- that is, whether it keeps itself to itself. Returns 1 if the fluxon can be contained in a 3-ball whose diameter is smaller than 1/3 the distance to the nearest foreign neighbor of any fluxel in the fluxon.
NUM cart_2d(NUM *x1, NUM *x2)
Returns the Cartesian distance between two points in a plane (expressed as 2-vectors).
NUM cart_3d(NUM *x1, NUM *x2)
Cartesian distance between two points in 3-space (expressed as 3-vectors).
NUM p_l_dist(NUM *p0, NUM *x0, NUM *x1)
Cartesian distance of closest approach between a point in 3-space (p0) and a line that passes through two other points (x0 and x1).
NUM p_ls_dist(NUM *p0, NUM *x0, NUM *x1)
Cartesian distance of closest approach between a point in 3-space (p0) and a line segment connecting two other points (x0 and x1). The distance of closest approach may be between p0 and one of the endpoints of the line segment, and therefore not perpendicular to the line (as in p_l_dist, above).
NUM l_l_dist(NUM *a0, NUM *b0, NUM *c0, NUM *d0)
Cartesian distance of closest approach between two arbitrary lines in 3-space. The lines are defined by two points each: a0-b0 and c0-d0. Note that the distance is the closest approach of the two complete lines, not merely the line segments between the defining points.
void p_ls_closest_approach(NUM *out, NUM *a0, NUM *b0, NUM *c0)
Location of closest approach between a point (c0) and a line segment (a0-b0) in 3-space.
void ls_closest_approach(NUM *out0, NUM *out1, NUM *a0, NUM *b0, NUM *c0, NUM *d0)
Returns the line segment of closest approach between two line segments (a0-b0 and c0-d0) in 3-space. On return, out0 contains the point on a0-b0, and out1 contains the point on c0-d0.
NUM ls_ls_dist(NUM *a0, NUM *b0, NUM *c0, NUM *d0)
Cartesian distance of closest approach between two line segments (see also ls_closest_approach, above).
Retained for historical purposes. Unused and deprecated.
NUM fl_segment_deluxe_dist(NUM *p0, NUM *p1, VERTEX *v0, VERTEX *v1)
Returns the wonky non-Cartesian distance of closest approach between two fluxels, and the points of closest Cartesian approach. The points p0 and p1 are stuffed with the points of closest Cartesian approach on the fluxels in question. This works by calling ls_closest_approach and then scaling the distance by the squared cosecant of the out-of-plane angle from v0's perpendicular plane.
This is the main distance calculating engine used by the hull algorithm and thus the magnetic pressure calculator.
These routines deal with lines in the plane and their intersections; they are building blocks for the hull algorithm that is the basis of the fluxon magnetic-pressure calculation.
Lines in 2-space are represented as triplets (3-vectors, sort of). The 0 and 1 element of the triplet contain the origin of the line (this is useful because it retains a datum from which the line is derived, rather than simpy a Y-intercept). The 2 element of the triplet contains the slope (dy/dx). This works OK because IEEE numbers represent the extended reals (including Infinity), so non-finite slope implies a vertical line.
int perp_bisector_2d(NUM *L_out, NUM *P, NUM *Q)
Takes two 2-points, and returns the origin and slope of their perpendicular bisector in the 3-array out On exit, L_out (3 elements) contains the 2-line specification for the bisector.
int intersection_2d(NUM *out, NUM *L1, NUM *L2)
Given two 2-D lines, return a 2-vector containing their intersection. If the lines are parallel, then the elements of out are left with the value NaN.
These routines deal with location and orientation of planes and points. They are building blocks for the find_simplex_by_location routine that is used for field value interpolation.
int above_plane(POINT3D A, POINT3D B, POINT3D C, POINT3D X)
Takes three noncolinear 3-points (A, B, and C) and treats them as describing an oriented plane segment (ABC widdershins). Returns 1 if X is above the plane, 0 if it is below it. (If X is coplanar, it is treated as above).
int in_simplex(POINT3D A, POINT3D B, POINT3D C, POINT3D D, POINT3D X)
Takes four noncoplanar 3-points (A, B, C, and D) and treats them as describing a simplex in 3-space (a tetrahedron). Returns 1 if X is inside the simplex, and 0 if it is outside. If it is on one of the faces it has about an even chance of being called inside or outside.
void project_n_fill(VERTEX *v, DUMBLIST *horde)
Given a VERTEX and a DUMBLIST of VERTICES, project each vertex in the DUMBLIST down to 2-D, using fl_segment_deluxe_dist to find the projected radius. The angle, radius, and scr fields of the members of the DUMBLIST are filled in, with the 2-D projected values. This is used as a precursor to hull finding: all candidate neighbors are projected into the perpendicular plane of the central fluxel being considered.
The scr field gets the 2-point that comes from projecting the vertex (in the requisite non-Euclidean manner, using fl_segment_deluxe_dist, above) into the central vertex's perpendicular plane. The a field gets a sleazy approximation to the angle in that plane (see atan2_oct, below), and the r field gets the radius of the projected point.
void sort_by_angle_2d(DUMBLIST *horde)
Given a projected-and-filled list of vertices, sort them by the a (in-plane angle) field. This was a key step in early hull calculations, but is no longer used. Sorting turned out to be a rate-limiting step, and the new hull algorithm does just fine without it.
void hull_2d(HULL_VERTEX *out, DUMBLIST *horde, DUMBLIST *rejects)
This is the original central hull-finding algorithm, now unused (but left in the code for historical and reference purposes). It requires the list of neighbor candidates (in horde) to be sorted by increasing anglein the plane. It accepts a large DUMBLIST containing vertices filled with (a,r) projection fields relative to a central vertex (which is considered to be at the origin), and sorts the elements into neighbors (left in the DUMBLIST) and rejects (deleted from horde and inserted into rejects). If you pass in NULL instead of a valid rejects DUMBLIST, then the rejects are simply ignored.
You must also prepare an array of HULL_VERTEX structures that is large enough to contain the output -- in practice, this means it must be as large as the input horde. On return, the horde contains only neighbors, in increasing angular order, and the HULL_VERTEX array contains corresponding hull vertices.
The main algorithm is based on a description by Preparata and Shamos, and relies on accepting or rejecting neighbor candidates depending on the relative angular order of the constructed candidate hull points. Potential neighbors whose hull vertices are in the wrong angular order are rejected. The algorithm is fast, given a sorted list of neighbor candidates -- but the sorting is a real dog.
void hull_2d_us(HULL_VERTEX *out, DUMBLIST *horde, DUMBLIST *rejects)
This is the current main hull-finding algorithm, hence the hot spot for the simulation as a whole. It has identical calling conventions to hull_2d above, but does not require that the horde be presorted. Instead, the neighbors are sorted on the fly using insertion sort. Even though insertion sort is notoriously slow on arrays (because it requires frequent recopying of all the data in the array), this is a big win because out of hundreds to thousands of candidate neighbors only (on average) about 7 are kept, so the copying is not a big hit compared to sorting the list of hundreds to thousands of candidates.
The hull algorithm is straightforward in concept but tricky to parse out of the code because of all the special cases.
int check_hullpoint(VERTEX *v, VERTEX *pv, HULL_VERTEX *phv, VERTEX *nv, HUL_VERTEX *nhv, HULL_VERTEX *scr)
This is a helper routine for hull_2d_us. It determines whether a particular candidate vertex (v) should be kept or rejected, based on the geometry of it and the next and previous vertices in the current hull-under-construction. pv and phv contain the previous (highest angle less then v->a) vertex and its corresponding hull point, and nv and nhv contain the next vertex (lowest angle greater than or equal to v->a) and its corresponding hull point. scr is a scratch hull point that gets populated with some fields if the point is to be kept.
The reason for breaking check_hullpoint out of the main hull_2d_us routine is that it is convenient to perform the (somewhat complex) test in several different branches of the main selection algorithm.
For rendering or interfacing with other modules, sometimes it's necessary to interpolate off the fluxon grid. This is accomplished by the method of simplexes: to find the value of a computed number at an arbitrary point in space, the value is interpolated from nearby VERTEXes. These routines aid that process.
VERTEX *find_vertex_by_location(POINT3D x, WORLD *w, VERTEX *v_hint, int global)
Walks through the WORLD to find the vertex closest to a given point in 3-space. You specify the point (as x) and the WORLD. Normally, a minimization-walk is done - the routine simply steps through the grid from a starting point in the direction of the desired point, until no local links can get it any closer. In that case, you can specify a v_hint location to start the search. (If you pass in 0, it instead starts at the root of the all-VERTEX tree). If you want an explicit global minimization search, feed a nonzero value into global, and the code will examine every VERTEX in the simulation and find the closest to the point. This is generally a lot slower, and in practice seems to be no better at finding the correct VERTEX, though in some pathological cases it might improve matters.
DUMBLIST *find_simplex_by_location(POINT3D x, WORLD *w, VERTEX *v, int global)
Given a location in 3-space, attempt to find a minimal simplex (tetrahedron) of vertices around it. Returns a DUMBLIST that contains the vertices of the simplex. If no such simplex is possible (e.g. the point is outside the convex hull of the simulation's VERTEX cloud) then the DUMBLIST gets fewer than four elements. The returned DUMBLIST is statically allocated and remains valid until the next call of find_simplex_by_location.
find_simplex_by_location works in a pretty boneheaded manner - it is not a true hull finder. It works as follows:
-
Identify the point to be sampled (call it X);
-
Find the closest vertex to the desired sample point (by calling
find_vertex_by_location, above), called P0; -
Among the neighbors of the closest vertex, find the one (called P1) such that the angle P0-X-P1 is as large as possible;
-
Among the nearest neighbors of the closest vertex, find the one that is farthest from the plane defined by the P0-X-P1 triangle (called P2)
-
Among neighbors of P0, P1, and P2, find the one that (A) forms a simplex (with P0, P1, and P2) that encloses X, and (B) minimizes the volume of the simplex, across neighbors.
That generally finds acceptable simplexes, but is not as good as a true Voronoi minimizer.
NUM interpolate_value_simplex(POINT3D x, DUMBLIST *dl, int val_offset)
Given a location and a simplex of VERTEXes that (hopefully) surround it, linearly interpolate a NUM value between the vertices onto the new location. Any NUM stored in a VERTEX can be interpolated; you supply the offset (in bytes) from the start of the VERTEX structure.
NUM interpolate_value(POINT3D x, WORLD *w, VERTEX *v, int global, int val_offset)
Convenience routine calls find_simplex_by_location and interpolate_value_simplex all in one go.
A few routines exist only to trade off unnecessary computational accuracy for speed in the hot spot. At least, for small enough values of "a few". Seriously, project_n_fill and hull_2d_us are the major hot spots of the whole simulation and need to be optimized for speed. They currently aren't, but simple drop-in replacements for expensive tasks (such as square root or trig functions) can potentially pay off big.
NUM atan2_oct( NUM y, NUM x )
Because the a field of each VERTEX is only used to sort by angle, and not for any quantitative purpose, any monotonically increasing function of angle will do. The arctan is computationally expensive, so atan2_oct gives a sleazy approximation by treating the unit circle as an octagon -- rather than integrated distance along the circumference of a unit circle, it returns integrated distance along the circumference of an octagon with vertices on the X and Y axes and the 45-degree lines, and total perimeter 2π. This is much faster because no convergent series is required.
atan2_oct isn't called directly; instead it is invoked with the macro ATAN2, which can be switched between atan2 and atan2_oct by a #define in geometry.h.