How to accelerate the intersection test between triangles and AABB

Question

I am learning the algorithm for voxelization of triangular meshes, which involves testing the intersection between triangles and AABB. I found the Separating Axis Theorem or SAT algorithm online, which requires the vertices of triangles and AABB to be projected onto 13 separation axes for judgment. Algorithms can work, but their efficiency is very low, and there may be a hundred times difference compared to Open3D. How should I improve them, or what better algorithms are available?

bool checkTriangleIntersectAABB(const std::vector<Eigen::Vector3d>& vTriangle, const Eigen::Vector3d& vCenter, const Eigen::Vector3d& vExtent)
{
    std::vector<Eigen::Vector3d> Axis;
    std::vector<Eigen::Vector3d> XYZNormal = { Eigen::Vector3d(1, 0, 0), Eigen::Vector3d(0, 1, 0), Eigen::Vector3d(0, 0, 1) };
    std::vector<Eigen::Vector3d> Point = { vTriangle[0] - vCenter, vTriangle[1] - vCenter, vTriangle[2] - vCenter };
    std::vector<Eigen::Vector3d> Edge = { Point[0] - Point[1], Point[1] - Point[2], Point[2] - Point[0] };

    Axis.push_back(XYZNormal[0]);
    Axis.push_back(XYZNormal[1]);
    Axis.push_back(XYZNormal[2]);
    for (int i = 0; i < XYZNormal.size(); i++)
    {
        for (int j = 0; j < Edge.size(); j++)
        {
            Axis.push_back(XYZNormal[i].cross(Edge[j]).normalized());
        }
    }
    Axis.push_back(Edge[0].cross(Edge[1]).normalized());

    for (int i = 0; i < Axis.size(); i++)
    {
        std::vector<double> PointProjection = { Point[0].dot(Axis[i]), Point[1].dot(Axis[i]), Point[2].dot(Axis[i]) };
        double CubeProjection = vExtent.x() * abs(XYZNormal[0].dot(Axis[i]))
            + vExtent.y() * abs(XYZNormal[1].dot(Axis[i]))
            + vExtent.z() * abs(XYZNormal[2].dot(Axis[i]));
        double MinPoint = std::min({ PointProjection[0], PointProjection[1], PointProjection[2] });
        double MaxPoint = std::max({ PointProjection[0], PointProjection[1], PointProjection[2] });
        if (std::max(-MaxPoint, MinPoint) > CubeProjection) return false;
    }
    return true;
}

std::shared_ptr<open3d::geometry::VoxelGrid> voxelizeInSurface(const open3d::geometry::TriangleMesh& vTriangleMesh, double vVoxelSize)
{
    _ASSERTE(vVoxelSize > 0);

    Eigen::AlignedBox3d BoundBox;
    for (auto& Vertex : vTriangleMesh.vertices_) BoundBox.extend(Vertex);
    Eigen::Vector3d MinBound = BoundBox.min(), BoundSize = BoundBox.max() - BoundBox.min();

    auto VoxelGrid = std::make_shared<open3d::geometry::VoxelGrid>();
    VoxelGrid->origin_ = MinBound;
    VoxelGrid->voxel_size_ = vVoxelSize;
    double HalfVoxelSize = vVoxelSize / 2;
    for (auto& Indexs : vTriangleMesh.triangles_)
    {
        std::vector<Eigen::Vector3d> Triangle = { vTriangleMesh.vertices_[Indexs[0]],
                                                 vTriangleMesh.vertices_[Indexs[1]],
                                                 vTriangleMesh.vertices_[Indexs[2]] };
        Eigen::AlignedBox3d TriangleBound;
        TriangleBound.extend(Triangle[0]);
        TriangleBound.extend(Triangle[1]);
        TriangleBound.extend(Triangle[2]);
        Eigen::Vector3i MinIndex = Eigen::floor((TriangleBound.min() - MinBound).array() / vVoxelSize).cast<int>();
        Eigen::Vector3i MaxIndex = Eigen::floor((TriangleBound.max() - MinBound).array() / vVoxelSize).cast<int>();
        for (int x = MinIndex[0]; x <= MaxIndex[0]; x++)
        {
            for (int y = MinIndex[1]; y <= MaxIndex[1]; y++)
            {
                for (int z = MinIndex[2]; z <= MaxIndex[2]; z++)
                {
                    Eigen::Vector3i Index(x, y, z);
                    Eigen::Vector3d Center = (Index.cast<double>() * vVoxelSize + MinBound).array() + HalfVoxelSize;
                    Eigen::Vector3d Extent(HalfVoxelSize, HalfVoxelSize, HalfVoxelSize);
                    if (checkTriangleIntersectAABB(Triangle, Center, Extent))
                    {
                        Eigen::Vector3d Color = (Center - MinBound).array() / BoundSize.array();
                        VoxelGrid->voxels_[Index] = open3d::geometry::Voxel(Index, Color);
                    }
                }
            }
        }
    }
    return VoxelGrid;
}

Is this a problem with algorithm logic, or is it caused by the use of std:: vector in the code, or the generation of copy constructors

Answer 1

The Akenine-Moller solution for triangle:AABB intersection using SAT is quite inefficient and also fails on degenerate triangles. But it's the solution that always comes up in searches so it's easy to focus on it.

An improvement can be found in the work by Schwarz and Siedel here: https://michael-schwarz.com/research/publ/files/vox-siga10.pdf

Their approach is quite elegant:

1. Check for an intersection between the plane of the triangle and the AABB.
2. Project the triangle and AABB onto the three coordinate planes (XY, YZ and ZX) and test for intersection.

You might find the rest of the paper interesting too, since it covers the topic that you're inquiring about.

However this is still has the limitation that positive intersections are only detected after all the potential disjoint cases have been checked. There is no early out in the cases where the triangle and the AABB intersect. So I use another approach which is faster than both the Akenine-Moller and Schwarz-Siedel:

1. Check for intersection between the AABB and the bounds of the triangle. Exit early if disjoint.
2. Test each of the three triangle edges against the AABB. Exit early with an intersection.
3. Check for intersection between the plane of triangle and the AABB. Exit if disjoint.
4. Test each of the four internal diagonal axes of the AABB against the triangle, for the cases where the box intersects the face of the triangle without touching any of its edges.

A Unity project with C# code for this solution as well as the other two can be found here: https://github.com/creijon/intersections

Answer 2

In order to reduce complexity you can use a rasterization algorithm.

• First, raster the longest edge: If it is the first, for example: from vTriangleMesh.vertices_[Indexs[0]] to vTriangleMesh.vertices_[Indexs[1]]
• Find the vertices' cell coordinates. For example: start at (i0,j0,k0) and end at (i1,j1,k1).
• You can raster from start cell to end cell by Bresenham's line drawing algorithm or similar (in their 3d version).
• After coloring the longest edge, raster line from each of the colored cells to the opposite vertex (vTriangleMesh.vertices_[Indexs[2]] if the first is the longest).
• To overcome the problem of rasterization algorithms, that they do not guarantee to color all cells that the line intersects, check the closets environment (8 cells) of each colored cell, for intersection with the triangle.