To deal with vectors and cuboids in a three-dimensional integer lattice.
§2. We will store 3-vectors in the obvious way:
typedef struct vector {
int x, y, z;
} vector;
The structure vector is accessed in 3/sm2, 3/hm, 3/em and here.
§3. Some useful constant vectors, including those pointing in each direction. Note that these are not of unit length — rather, they are the ideal grid offsets on the map we will eventually draw.
vector Zero_vector = {0, 0, 0};
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 3/sm2, 3/hm, 3/em and here.
vector Geometry::vec(int x, int y, int z) {
vector R;
R.x = x; R.y = y; R.z = z;
return R;
}
The function Geometry::vec is used in 3/sm2 (§16.1, §8.31, §8.32.1.1, §8.32.1.2, §8.32.1.3), 3/hm (§9, §9.3.1.1, §9.3.2.1, §9.3.3.1).
§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;
}
The function Geometry::vec_eq is used in 3/sm2 (§8.11, §21, §22, §33), 3/hm (§1.3.4, §1.3.4.1, §1.3.4.2, §2.1).
The function Geometry::vec_lateral is used in 3/sm2 (§8.10).
§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;
}
The function Geometry::vec_plus is used in §12, §15, 3/sm2 (§8.21, §8.26, §8.27, §21, §23.1.4.1, §33, §37.1, §8.31, §8.32.1.2, §8.32.1.3), 3/hm (§1.3.4.2, §2, §2.1, §3).
The function Geometry::vec_minus is used in §14, §15, 3/sm2 (§21, §22, §23.1.4.1, §35.1, §8.32.1.2, §8.32.1.3, §8.33), 3/hm (§1.3.4.1).
The function Geometry::vec_negate is used in 3/sm2 (§8.32.1.3.1).
The function Geometry::vec_scale is used in 3/sm2 (§23.1.4.1, §33, §37.1), 3/hm (§1.3.4.1).
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)));
}
The function Geometry::vec_length_squared is used in 3/sm2 (§21).
The function Geometry::vec_length is used in §9.
§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);
}
The function Geometry::vec_angular_separation is used in 3/sm2 (§21, §22).
§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;
}
}
The function Geometry::empty_cuboid is used in 3/sm2 (§8, §12, §8.18, §19, §8.32, §8.34).
The function Geometry::adjust_cuboid is used in §11, 3/sm2 (§19, §8.34).
§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;
}
}
}
The function Geometry::merge_cuboid is used in 3/sm2 (§8.32.1).
§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);
}
}
The function Geometry::cuboid_translate is used in 3/sm2 (§8.21).
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;
}
The function Geometry::within_cuboid is used in §14, 3/sm2 (§16), 3/hm (§2).
§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;
}
The function Geometry::cuboid_index is used in 3/sm2 (§14, §16, §16.1), 3/hm (§1).
The function Geometry::cuboid_volume is used in 3/sm2 (§16.1), 3/hm (§1.1).
§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;
}
}
The function Geometry::thicken_cuboid is used in 3/sm2 (§16.1).