如何加速三角形与AABB的相交测试

Question

我正在学习三角形网格体素化的算法，其中涉及测试三角形和 AABB 之间的交集。我在网上找到了分离轴定理或者SAT算法，需要将三角形和AABB的顶点投影到13个分离轴上进行判断。算法可以工作，但是效率很低，和Open3D相比可能有一百倍的差别。我应该如何改进它们，或者有什么更好的算法可用。谢谢。

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

这是算法逻辑的问题，还是代码中使用了std::vector，或者复制构造函数的生成导致的