WIP: Replace collision code with faster, vectorisable tests and cache-friendly datastructures.

[qazwsxal]

Hello all, finally getting around to this idea after having it pinging around in my head for a few years. There's some pretty major bottlenecks in the collision code that I think we could massively improve, so I went ahead and clauded up a bunch of optimisations:

BSP node traversal and leaf layout are both cache-unfriendly

First of all, let's look at the structs, and how data gets laid out in them:

// code/model/model.h (master)
struct bsp_collision_node {
    vec3d min;
    vec3d max;
    int   back;
    int   front;
    int   leaf;   // >= 0 if this is a leaf, index into leaf_list
};                // 36 bytes

struct bsp_collision_tree {
    bsp_collision_node *node_list;
    int                 n_nodes;
    bsp_collision_leaf *leaf_list;
    int                 n_leaves;
    model_tmap_vert    *vert_list;
    vec3d              *point_list;
    SCP_vector<vec3d>   poly_centers;
    int                 n_verts;
    bool                used;
};

and the current BSP leaf data format:

// code/model/model.h (master)
struct bsp_collision_leaf {
    vec3d plane_norm;
    int   vert_start;   // index into tree->vert_list
    ubyte num_verts;
    ubyte tmap_num;
    int   next;         // linked list -- next == index+1 in practice, pure pointer-chase
};

BSP traversal is branchy, and doesn't cull when possible

And the BSP bbox traversal code, I've stripped out a bit for brevity:

// code/model/modelcollide.cpp (master), model_collide_bsp() — shape, not verbatim
void model_collide_bsp(int node_index)
{
    bsp_collision_node *node = &tree->node_list[node_index];

    if (!mc_ray_boundingbox(&node->min, &node->max, ...))
        return;

    if (node->leaf >= 0) {
        model_collide_bsp_poly(node->leaf);
        return;
    }

    model_collide_bsp(node->front);   // always both children --
    model_collide_bsp(node->back);    // no tmin vs. best-hit-so-far cutoff
}

So there's a few issues here:

1. The way we load the BSP means that we append nodes breadth first, so node->front and node->back end up indexing deep into the array, vs hopefully one of them being nicely cache local.
2. We're always checking both front and back, why aren't they stored close together in memory? If a node contained its childs bounding boxes, we'd avoid having to index them seperately, and could load all of this into registers .
3. There's no check to see if a closer intersection is actually possible, we shouldn't enter a bounding box if all possible collisions are worse than our current best hit.

Leaf nodes involve pointer chasing and aren't nicely vectoriseable

Testing one leaf means:
2. look up Mc_pm->maps[tmap_num] for the texture, if it's invisible, check submodel[].flags for the invisible-face path before bothering with anything else
3. for each vertex, index into tree->vert_list (its own array)
4. for each of those, index into a separate tree->point_list (a scatter, not a sequential read)
5. run the face/sphereline check against this leaf
6. follow next (pointer chase) to the next leaf in the chain, and repeat

And here's model_collide_bsp_poly() itself, which is what actually does all this per leaf:

// code/model/modelcollide.cpp (master), model_collide_bsp_poly() -- shape, not verbatim
void model_collide_bsp_poly(bsp_collision_tree *tree, int leaf_index)
{
    vec3d *points[TMAP_MAX_VERTS];

    while (leaf_index >= 0) {                          // walk the linked list
        bsp_collision_leaf *leaf = &tree->leaf_list[leaf_index];

        if (leaf->tmap_num < MAX_MODEL_TEXTURES) {
            if (texture_is_invisible(leaf->tmap_num) && !collide_invisible_flag_set())
                return;                                 // bail out of the WHOLE chain, not just this leaf
        }

        for (int i = 0; i < leaf->num_verts; ++i) {
            int vert_num = tree->vert_list[leaf->vert_start + i].vertnum;
            points[i] = &tree->point_list[vert_num];    // scatter: one lookup per vertex
            // ... uvlist[i] filled in from tree->vert_list[...] here too, omitted
        }

        // ... real code branches on tmap_num < MAX_MODEL_TEXTURES here (flat/textured poly),
        // passing either nullptr/-1 or uvlist/tmap_num into the calls below -- omitted
        if (Mc->flags & MC_CHECK_SPHERELINE)
            mc_check_sphereline_face(leaf->num_verts, points, points[0], &leaf->plane_norm, ...);
        else
            mc_check_face(leaf->num_verts, points, points[0], &leaf->plane_norm, ...);

        leaf_index = leaf->next;                        // on to the next leaf in the chain
    }
}

This is pretty horrible when it comes to vectorisation, what we'd really want to do here is ray/sphereline intersections with a bunch of triangles in a clean, branchless loop to encourage autovectorisation. However we don't have that here unfortunately!
I'd also note here the num_verts field. The existing collision leaf code handles n-gon collisions, but requires specifying the plane that the polygon will be projected onto, this results in quite odd collisions, particularly if the plane_norm ends up with a strange value, in particular I noticed turret02 in SC_Asura.pof has triangles with really odd plane_norm values that doesn't seem consistent between faces either. It doesn't look intentional.

New implementation:

The datastructures are as follows:

// code/model/modelbvh.h (this branch)
constexpr int BVH_N = 4;

struct bvh_node {  // BVH_N=4 ->128 Bytes
    float minx[BVH_N], miny[BVH_N], minz[BVH_N];
    float maxx[BVH_N], maxy[BVH_N], maxz[BVH_N];
    int32_t child[BVH_N];
    int32_t count[BVH_N];   // count[i] > 0 -> leaf; count[i] == 0 -> internal
};

struct bvh_tri_indices {
    uint32_t i0, i1, i2;
};

struct bvh_tree {
    SCP_vector<bvh_node> nodes;   // depth-first / pre-order layout

    SCP_vector<vec3d> verts;              // shared, deduplicated vertex pool -- interleaved x/y/z
    SCP_vector<bvh_tri_indices> tris;     // per-triangle index triple -- interleaved i0/i1/i2
    SCP_vector<int> tmap_num;
    SCP_vector<int> original_index;
    SCP_vector<int> leaf_index;
    SCP_vector<bvh_uv> uv0, uv1, uv2;
    SCP_vector<vec3d> normal;             // precomputed unit face normal per triangle
    SCP_vector<bvh_bsphere> bsphere;      // precomputed per-triangle bounding sphere

    int root = 0;
};

We parse the data the same way as before to pull out all information present in the BSP, and instead of using it for collision, build an N-way bounding volume hierarchy with a more friendly structure. This is a three pass construction over a "triangle soup". I use Surface Area Heuristic to determine how to split (this is a pretty standard offline construction approach).

// code/model/modelbvh.cpp (this branch) -- build shape, not verbatim
bvh_tree bvh_build(SCP_vector<bvh_triangle> triangles)
{
    // 1. binary SAH build: top-down, 
    BinaryNode *binary_root = sah_build_binary(triangles);

    // 2. greedily collapse pairs of binary nodes into BVH_N-wide nodes, lowest
    //    SAH-cost collapse first -- standard production technique, 
    // much simpler to implement correctly vs N-way SAH split
    collapse_to_n_wide(binary_root, BVH_N);

    // 3. flatten depth-first / pre-order -- a node's subtree immediately follows
    //    it in the output array, specifically to avoid the BSP tree's breadth-
    //    first cache problem
    bvh_tree tree;
    flatten_depth_first(binary_root, tree.nodes);
    return tree;
}

There's a few advantages to laying out the BVH tree like this. The major one, is that we can test all 4 children of a single node while their data is fresh in cache. The 4-wide layout and child-bbox packing means one node touches at most 2 cache lines with and is often fetched with a single burst vs. up to 4 separate scattered fetches for the BSP equivalent.

BVH Traversal

// code/model/modelbvh.h, bvh_visit_triangles() (this branch)
template <typename Visitor>
void bvh_visit_triangles(const bvh_tree& tree, const vec3d& origin, const vec3d& dir,
    float t_max, float radius, Visitor&& visit)
{
    struct StackEntry {
        int32_t index;
        int32_t count; // >0: leaf (index is a triangle-array start); ==0: internal node
        float tmin;    // the tmin this entry was queued with -- see below
    };
    StackEntry stack[64];
    int sp = 0;
    stack[sp++] = {tree.root, 0, 0.0f};

    while (sp > 0) {
        StackEntry entry = stack[--sp];
        if (entry.tmin > t_max)
            continue;                          // stale: a nearer hit tightened t_max
                                                 // after this entry was queued
        if (entry.count > 0) {
            visit(entry.index, entry.count, t_max);   // t_max passed by ref --
            continue;                                 // visitor can tighten it
        }

        const bvh_node& node = tree.nodes[entry.index];
        StackEntry candidates[BVH_N];
        int num_candidates = 0;
        for (int i = 0; i < BVH_N; ++i) {
            if (node.child[i] < 0) continue;
            float tmin;
            if (!bvh_detail::ray_aabb_visit_tmin(origin, inv_dir,
                    {node.minx[i], node.miny[i], node.minz[i]},
                    {node.maxx[i], node.maxy[i], node.maxz[i]}, t_max, radius, tmin))
                continue;                       // skip, no descent at all
            candidates[num_candidates++] = {node.child[i], node.count[i], tmin};
        }

        sort_ascending_by_tmin(candidates, num_candidates);      // <=4 elements
        for (int i = num_candidates - 1; i >= 0; --i)            // push farthest-first
            stack[sp++] = candidates[i];                         // so nearest pops next
    }
}

N.B. MC_COLLIDE_ALL is handled properly in the code, these snippets are modified not to include it to focus on runtime collision detection.

Visitor here is just a template argument for what kind of leaf intersection we're testing, ray or sphere. It means we don't have to check what kind every time we hit a leaf node, branch prediction is usually OK with this, but better safe than sorry. The main model_collide() function does dispatch based on the existing flags. Function pointers work too, but I'd rather the potential for inlining here for a little performance boost.
I've intentionally switched to using an explicit LIFO stack rather than relying on a recursive definition, this works fine and avoid the overhead from recursive calls. The max stack depth is 64 here, which is pretty big compared to what the real world use cases show (19 on BPC ships, worth checking on BTA's and MediaVPs though), but it's data dependent so I've left a large overhead in case someone builds pathological ships.

I investigated swapping out the looped box test for an auto-vectorisable form, but the differences were negilible.bvh_detail::ray_aabb_visit_tmin also contains an early out and doesn't count a bbox as a "hit" if t_max is already shorter than the bbox intersection.

I also sort the 4 child bboxes so that they're added to the stack nearest last, only 4 nodes here, so sorting is fast, and forces us to explore the best potential node first. This ensures we first visit the node that if hit, will result in the greatest culling of other nodes.

Leaf Traversal

Ray-triangle intersections in ray_triangle_leaf_simd use Möller-Trumbore over a SIMD_WIDTH-wide batch. This is unconditional/branchless so it autovectorizes; I checked the disassembly and got real packed xmm/ymm with widths of 4 and 8 respectively.

The sphereline (radius-swept) leaf test in mc_check_triangle_sphereline_face is scalar, not SIMD-batched — every attempt at vectorizing this one measured slower and was reverted (see below); it still gets the same early-out discipline the ray path gets from front-to-back traversal:

The bvh_tree struct above already shows the leaf-side representation. There's a few points to note here:

• Leaf size isn't fixed, SAH construction can decide that splitting a node isn't worth it, there's a couple of very high (80-120) triangle count leaf nodes in edge-case geometry where large flat triangles next to a bunch of small ones don't align with axes.
• Might be worth checking if there's any improvements to the algorithm that can be made there, this is a pretty standard arrangement when modelling.
• Any node with <=16 triangles in it is automatically converted into a leaf node. At a certain point, linear scans are just worth doing, and 16 seemed the sweet spot for sphereline intersections.
• N-gons are triangulated, we test ray and sphere collisions against the raw geometry, not n-gons projected onto .POF surface norms, triangulation matches generation for rendering.
• This does lead to some minor hit-location differences.
• Tree nodes store child and count, if count > 0, then the triangles are taken from tris in the [child, child+count) range.
• tris corresponds to a densely packed vector of vertex indices, three per tri.
• These then point into the verts array, three float32 vals per entry. We do a deduplication pass here, sharp edge rendering might need split edges and duplicated verts, but for collision purposes, duplicating data is pointless, and as these

Submodel walk: full linear search vs. top-level LBVH

• Whole-model queries walk every submodel via a plain linked-list recursion, there's no spatial index or bounds check before recursing — aibig.cpp self-documents this as ~10% of AI frametime
• A spatial index would be nice here! However, submodels can rotate/translate at runtime (turrets, animations, Lua), so a static once-at-load spatial tree isn't safe.
• So instead, I use a top-level LBVH (Morton-code centroid BVH tree) rebuilt fresh once per frame per ship instance, This takes between 0.5-8.8μs to rebuild depending on submodel count (tested up to 200), and is re-used across all collision checks.
• Synthetic tests showed a good performance, boost here, but real-mission testing needs to be done.
• I also added skips to avoid doing a collision test entirely if a submodel's AABB is already farther than an existing hit.

Shields: just another subobject

• Currently shields have a bespoke byte-blob walker mc_check_sldc() that uses hard-coded offset arithmetic, this has an order-dependent bug (front/back || short-circuit skips the back subtree)
• Shields are pure triangles with no submodel hierarchy, so I build a bvh_tree as if they were just another submodel and go through the exact same code (model_collide_bvh_triangle())
• No more collision bug, and now we dont have two subtly different implementations of tree-walking collision code