如何在C++中高效将3D三角形光栅化至均匀体素网格？

核心思路 #

把2D扫描线光栅化扩展到3D的核心是按轴分层切片：选择三角形跨度最大的轴作为扫描轴，将3D问题拆解为一系列平行的2D截面光栅化任务。每个截面是三角形与体素网格某一层平面的交线，对这条交线做2D扫描线填充，最终映射到对应的3D体素。这种方法复用成熟的2D扫描线逻辑，同时通过选最长轴最小化切片数量，保证效率。

步骤拆解与代码实现 #

1. 基础定义与预处理 #

先定义顶点、体素网格的基础结构，预处理三角形关键参数：

#include <vector>
#include <algorithm>
#include <cmath>

// 3D顶点结构体
struct Vec3 {
    float x, y, z;
    Vec3(float x=0, float y=0, float z=0) : x(x), y(y), z(z) {}
    Vec3 operator-(const Vec3& other) const { return Vec3(x-other.x, y-other.y, z-other.z); }
    Vec3 cross(const Vec3& other) const {
        return Vec3(
            y*other.z - z*other.y,
            z*other.x - x*other.z,
            x*other.y - y*other.x
        );
    }
    float dot(const Vec3& other) const { return x*other.x + y*other.y + z*other.z; }
    Vec3 normalize() const {
        float len = std::sqrt(x*x + y*y + z*z);
        return Vec3(x/len, y/len, z/len);
    }
};

// 3D体素网格，一维数组模拟（索引：x + y*width + z*width*height）
class VoxelGrid {
public:
    int width, height, depth;
    std::vector<float> data; // 存储密度/颜色值

    VoxelGrid(int w, int h, int d) : width(w), height(h), depth(d), data(w*h*d, 0.0f) {}

    float& get(int x, int y, int z) {
        return data[x + y*width + z*width*height];
    }
};

2. 确定扫描轴与分层范围 #

计算三角形包围盒，选跨度最大的轴作为扫描轴，减少切片数量：

// 获取三角形的包围盒
void getTriangleAABB(const Vec3& v0, const Vec3& v1, const Vec3& v2, 
                     float& minX, float& maxX, float& minY, float& maxY, float& minZ, float& maxZ) {
    minX = std::min({v0.x, v1.x, v2.x});
    maxX = std::max({v0.x, v1.x, v2.x});
    minY = std::min({v0.y, v1.y, v2.y});
    maxY = std::max({v0.y, v1.y, v2.y});
    minZ = std::min({v0.z, v1.z, v2.z});
    maxZ = std::max({v0.z, v1.z, v2.z});
}

// 选择扫描轴（0=X,1=Y,2=Z）
int selectScanAxis(float dx, float dy, float dz) {
    if (dx >= dy && dx >= dz) return 0;
    if (dy >= dx && dy >= dz) return 1;
    return 2;
}

3. 分层计算截面交线 #

对每个体素层，计算三角形与该层平面的交线（若相交）：

// 计算三角形与某平面的交线（planeVal是扫描轴对应的坐标值）
bool getTriangleSlice(const Vec3& v0, const Vec3& v1, const Vec3& v2, int axis, float planeVal,
                      Vec3& p0, Vec3& p1) {
    // 根据轴统一坐标获取逻辑
    auto getCoord = [&](const Vec3& v) {
        switch(axis) {
            case 0: return v.x;
            case 1: return v.y;
            default: return v.z;
        }
    };

    std::vector<std::pair<Vec3, Vec3>> edges = {{v0,v1}, {v1,v2}, {v2,v0}};
    std::vector<Vec3> intersections;

    for (auto& edge : edges) {
        float c0 = getCoord(edge.first);
        float c1 = getCoord(edge.second);
        if ((c0 <= planeVal && c1 >= planeVal) || (c0 >= planeVal && c1 <= planeVal)) {
            // 线性插值计算交点
            float t = (planeVal - c0) / (c1 - c0);
            Vec3 p(
                edge.first.x + t*(edge.second.x - edge.first.x),
                edge.first.y + t*(edge.second.y - edge.first.y),
                edge.first.z + t*(edge.second.z - edge.first.z)
            );
            intersections.push_back(p);
        }
    }

    if (intersections.size() == 2) {
        p0 = intersections[0];
        p1 = intersections[1];
        return true;
    }
    return false;
}

4. 2D扫描线填充截面 #

将3D交点投影到垂直于扫描轴的2D平面，用扫描线算法填充对应体素，更新3D体素值：

// 2D扫描线填充线段到体素层，更新3D体素
void rasterizeSlice(const Vec3& p0, const Vec3& p1, int axis, int sliceIdx, VoxelGrid& grid, float value, const Vec3& triNormal, const Vec3& triV0) {
    // 投影到2D平面
    float x0, y0, x1, y1;
    switch(axis) {
        case 0: // X轴扫描，投影到Y-Z
            x0 = p0.y; y0 = p0.z;
            x1 = p1.y; y1 = p1.z;
            break;
        case 1: // Y轴扫描，投影到X-Z
            x0 = p0.x; y0 = p0.z;
            x1 = p1.x; y1 = p1.z;
            break;
        default: // Z轴扫描，投影到X-Y
            x0 = p0.x; y0 = p0.y;
            x1 = p1.x; y1 = p1.y;
            break;
    }

    // 保证y0 <= y1
    if (y0 > y1) {
        std::swap(x0, x1);
        std::swap(y0, y1);
    }

    float yStep = 1.0f; // 体素边长，假设为1
    float xSlope = (y1 - y0) < 1e-6 ? 0 : (x1 - x0) / (y1 - y0);
    float currentX = x0;

    // 遍历每个Y对应的体素行
    for (float y = y0; y <= y1; y += yStep) {
        int yVoxel = static_cast<int>(std::floor(y));
        if (yVoxel < 0 || yVoxel >= (axis==0 ? grid.height : (axis==1 ? grid.depth : grid.height))) continue;

        float xStart = std::min(currentX, currentX + xSlope*yStep);
        float xEnd = std::max(currentX, currentX + xSlope*yStep);
        int xStartVoxel = static_cast<int>(std::floor(xStart));
        int xEndVoxel = static_cast<int>(std::floor(xEnd));

        // 填充X范围内的体素
        for (int x = xStartVoxel; x <= xEndVoxel; x++) {
            if (x < 0 || x >= (axis==0 ? grid.width : (axis==1 ? grid.width : grid.width))) continue;

            // 计算体素中心到三角形的距离，加权更新密度
            Vec3 voxelCenter;
            switch(axis) {
                case 0: voxelCenter = Vec3(sliceIdx+0.5f, x+0.5f, yVoxel+0.5f); break;
                case 1: voxelCenter = Vec3(x+0.5f, sliceIdx+0.5f, yVoxel+0.5f); break;
                default: voxelCenter = Vec3(x+0.5f, yVoxel+0.5f, sliceIdx+0.5f); break;
            }
            float distance = std::abs(triNormal.dot(voxelCenter - triV0));
            float weight = 1.0f / (1.0f + distance); // 距离越近权重越高

            // 更新3D体素值
            switch(axis) {
                case 0: grid.get(x, yVoxel, sliceIdx) += value * weight; break;
                case 1: grid.get(x, sliceIdx, yVoxel) += value * weight; break;
                default: grid.get(x, yVoxel, sliceIdx) += value * weight; break;
            }
        }

        currentX += xSlope*yStep;
    }
}

5. 主光栅化函数 #

整合所有步骤，完成3D三角形到体素网格的光栅化：

void rasterizeTriangle3D(const Vec3& v0, const Vec3& v1, const Vec3& v2, VoxelGrid& grid, float densityValue) {
    float minX, maxX, minY, maxY, minZ, maxZ;
    getTriangleAABB(v0, v1, v2, minX, maxX, minY, maxY, minZ, maxZ);

    int axis = selectScanAxis(maxX-minX, maxY-minY, maxZ-minZ);

    // 转换为体素索引范围
    int startIdx = static_cast<int>(std::floor(
        axis==0 ? minX : (axis==1 ? minY : minZ)
    ));
    int endIdx = static_cast<int>(std::floor(
        axis==0 ? maxX : (axis==1 ? maxY : maxZ)
    ));

    // 提前计算三角形法向量，用于距离加权
    Vec3 edge1 = v1 - v0;
    Vec3 edge2 = v2 - v0;
    Vec3 triNormal = edge1.cross(edge2).normalize();

    // 遍历每个体素层
    for (int idx = startIdx; idx <= endIdx; idx++) {
        float planeVal = idx + 0.5f; // 体素中心坐标
        Vec3 p0, p1;
        if (getTriangleSlice(v0, v1, v2, axis, planeVal, p0, p1)) {
            rasterizeSlice(p0, p1, axis, idx, grid, densityValue, triNormal, v0);
        }
    }
}

优化技巧 #

• 整数运算替代浮点：若体素网格是整数坐标对齐，可将所有浮点坐标乘以体素分辨率转为整数，减少浮点运算开销。
• 跳过空切片：提前判断三角形包围盒与体素网格的交集，只处理有重叠的层，避免无效循环。
• 批量连续体素填充：对连续的X范围，直接批量更新数组，减少循环次数，提升缓存命中率。

内容的提问来源于stack exchange，提问作者Spektre