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;

§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;

§5. Vectors.

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;
}

§8. 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)));
}

§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);
    }
}

§13.

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;
    }
}