Simple mathematical functions for dealing with vectors and cuboids in a three-dimensional integer lattice.
§1. We will store 3-vectors in the obvious way:
typedef struct vector { int x, y, z; } vector;
- The structure vector is accessed in 4/sm, 4/rhm, 4/rem and here.
§2. Some useful constant vectors:
vector Zero_vector = {0, 0, 0}; vector Geometry::zero(void) { return Zero_vector; }
§3. Now for vectors pointing in each direction. Note that not all of these are of unit length — rather, they are the ideal grid offsets on the map we will eventually draw.
vector N_vector = {0, 1, 0}; vector NE_vector = {1, 1, 0}; vector NW_vector = {-1, 1, 0}; vector S_vector = {0, -1, 0}; vector SE_vector = {1, -1, 0}; vector SW_vector = {-1, -1, 0}; vector E_vector = {1, 0, 0}; vector W_vector = {-1, 0, 0}; vector U_vector = {0, 0, 1}; vector D_vector = {0, 0, -1};
§4. A cuboid is a volume of space with opposing corners at integer grid positions which form a tightest-possible bounding box around a finite number of points (of size population); when this is 0, of course, the corners are meaningless and are by convention at the origin.
typedef struct cuboid { int population; struct vector corner0; struct vector corner1; } cuboid;
- The structure cuboid is accessed in 4/sm, 4/rhm, 4/rem and here.
vector Geometry::vec(int x, int y, int z) { vector R; R.x = x; R.y = y; R.z = z; return R; }
§6. A vector is "lateral" if lies in the \(x\)-\(y\) plane.
int Geometry::vec_eq(vector U, vector V) { if ((U.x == V.x) && (U.y == V.y) && (U.z == V.z)) return TRUE; return FALSE; } int Geometry::vec_lateral(vector V) { if ((V.x == 0) && (V.y == 0)) return FALSE; return TRUE; }
§7. The vector space operations:
vector Geometry::vec_plus(vector U, vector V) { vector R; R.x = U.x + V.x; R.y = U.y + V.y; R.z = U.z + V.z; return R; } vector Geometry::vec_minus(vector U, vector V) { vector R; R.x = U.x - V.x; R.y = U.y - V.y; R.z = U.z - V.z; return R; } vector Geometry::vec_negate(vector V) { vector R; R.x = -V.x; R.y = -V.y; R.z = -V.z; return R; } vector Geometry::vec_scale(int lambda, vector V) { vector R; R.x = lambda*V.x; R.y = lambda*V.y; R.z = lambda*V.z; return R; }
int Geometry::vec_length_squared(vector V) { return V.x*V.x + V.y*V.y + V.z*V.z; } float Geometry::vec_length(vector V) { return (float) (sqrt(Geometry::vec_length_squared(V))); }
§9. Angles. We compute unit vectors in the D and E directions and then the squared length of their difference. This is a fairly sharply increasing function of the absolute value of the angular difference between D and E, and is such that if the angles are equal then the result is zero; and it's cheap to compute. So although it might seem nicer to calculate actual angles, this is better.
float Geometry::vec_angular_separation(vector E, vector D) { float E_distance = Geometry::vec_length(E); float uex = E.x/E_distance, uey = E.y/E_distance, uez = E.z/E_distance; float D_distance = Geometry::vec_length(D); float udx = D.x/D_distance, udy = D.y/D_distance, udz = D.z/D_distance; return (uex-udx)*(uex-udx) + (uey-udy)*(uey-udy) + (uez-udz)*(uez-udz); }
§10. Cuboids. To form a populated cuboid, first request an empty one, and then adjust it for each vector to join the population.
cuboid Geometry::empty_cuboid(void) { cuboid C; C.population = 0; C.corner0 = Zero_vector; C.corner1 = Zero_vector; return C; } void Geometry::adjust_cuboid(cuboid *C, vector V) { if (C->population++ == 0) { C->corner0 = V; C->corner1 = V; } else { if (V.x < C->corner0.x) C->corner0.x = V.x; if (V.x > C->corner1.x) C->corner1.x = V.x; if (V.y < C->corner0.y) C->corner0.y = V.y; if (V.y > C->corner1.y) C->corner1.y = V.y; if (V.z < C->corner0.z) C->corner0.z = V.z; if (V.z > C->corner1.z) C->corner1.z = V.z; } }
§11. The following expands \(C\) minimally so that it contains \(X\).
void Geometry::merge_cuboid(cuboid *C, cuboid X) { if (X.population > 0) { if (C->population == 0) { *C = X; } else { Geometry::adjust_cuboid(C, X.corner0); Geometry::adjust_cuboid(C, X.corner1); C->population += X.population - 2; } } }
§12. Here we shift an entire cuboid over (assuming all of the points inside it have made the same shift).
void Geometry::cuboid_translate(cuboid *C, vector D) { if (C->population > 0) { C->corner0 = Geometry::vec_plus(C->corner0, D); C->corner1 = Geometry::vec_plus(C->corner1, D); } }
int Geometry::within_cuboid(vector P, cuboid C) { if (C.population == 0) return FALSE; if (P.x < C.corner0.x) return FALSE; if (P.x > C.corner1.x) return FALSE; if (P.y < C.corner0.y) return FALSE; if (P.y > C.corner1.y) return FALSE; if (P.z < C.corner0.z) return FALSE; if (P.z > C.corner1.z) return FALSE; return TRUE; }
§14. Suppose we have a one-dimensional array whose entries correspond to the integer grid positions within a cuboid (including its faces and corners). The following returns \(-1\) if a point is outside the cuboid, or returns the index if it is.
int Geometry::cuboid_index(vector P, cuboid C) { if (Geometry::within_cuboid(P, C) == FALSE) return -1; vector O = Geometry::vec_minus(P, C.corner0); int width = C.corner1.x - C.corner0.x + 1; int height = C.corner1.y - C.corner0.y + 1; return O.x + O.y*width + O.z*width*height; } int Geometry::cuboid_volume(cuboid C) { if (C.population == 0) return 0; int width = C.corner1.x - C.corner0.x + 1; int height = C.corner1.y - C.corner0.y + 1; int depth = C.corner1.z - C.corner0.z + 1; return width*height*depth; }
§15. Thickening a cuboid is a little more than adjusting; we give it some extra room. (The result is thus no longer minimally bounding, but we sacrifice that.)
void Geometry::thicken_cuboid(cuboid *C, vector V, vector S) { if (C->population++ == 0) { C->corner0 = Geometry::vec_minus(V, S); C->corner1 = Geometry::vec_plus(V, S); } else { if (V.x < C->corner0.x) C->corner0.x = V.x - S.x; if (V.x > C->corner1.x) C->corner1.x = V.x + S.x; if (V.y < C->corner0.y) C->corner0.y = V.y - S.y; if (V.y > C->corner1.y) C->corner1.y = V.y + S.y; if (V.z < C->corner0.z) C->corner0.z = V.z - S.z; if (V.z > C->corner1.z) C->corner1.z = V.z + S.z; } }