2019-02-05 02:44:07 +02:00
|
|
|
[Geometry::] Spatial Geometry.
|
|
|
|
|
2020-01-24 02:10:36 +02:00
|
|
|
To deal with vectors and cuboids in a three-dimensional integer lattice.
|
2019-02-05 02:44:07 +02:00
|
|
|
|
|
|
|
@ We will store 3-vectors in the obvious way:
|
|
|
|
|
|
|
|
=
|
|
|
|
typedef struct vector {
|
|
|
|
int x, y, z;
|
|
|
|
} vector;
|
|
|
|
|
|
|
|
@ 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};
|
|
|
|
|
2021-06-27 18:04:28 +03:00
|
|
|
=
|
|
|
|
vector Geometry::zero(void) {
|
|
|
|
return Zero_vector;
|
|
|
|
}
|
|
|
|
|
|
|
|
@
|
|
|
|
|
|
|
|
=
|
2019-02-05 02:44:07 +02:00
|
|
|
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};
|
|
|
|
|
|
|
|
@ 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;
|
|
|
|
|
|
|
|
@h Vectors.
|
|
|
|
|
|
|
|
=
|
|
|
|
vector Geometry::vec(int x, int y, int z) {
|
|
|
|
vector R;
|
|
|
|
R.x = x; R.y = y; R.z = z;
|
|
|
|
return R;
|
|
|
|
}
|
|
|
|
|
|
|
|
@ 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 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;
|
|
|
|
}
|
|
|
|
|
|
|
|
@h Lengths.
|
|
|
|
|
|
|
|
=
|
|
|
|
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)));
|
|
|
|
}
|
|
|
|
|
|
|
|
@h 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);
|
|
|
|
}
|
|
|
|
|
|
|
|
@h 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 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;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
@ 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;
|
|
|
|
}
|
|
|
|
|
|
|
|
@ 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;
|
|
|
|
}
|
|
|
|
|
|
|
|
@ 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;
|
|
|
|
}
|
|
|
|
}
|