Performance improvements

[GridEyes-2010]

I've made a few code improvements to your CPU version that increase efficiency by about 10 times The test results in the same environment are as follows： 646464--->5ms 128128128--->7ms 256256256--->19ms 512512512--->121ms 102410241024--->728ms 204820482048---->5576ms

[Forceflow]

Cool, would be interesting to see these improvements in a PR?

[GridEyes-2010]

void Voxelizer::cpu_voxelize_mesh_solid( const voxinfo& info, Triangle_mesh& mesh, unsigned int* voxel_table, bool morton_order) {

    const Vec3f& move_min = info.pMin;
    const int uzero = 0;
    int nvert = (int)mesh.n_vertices();
    auto points = mesh.points();
    std::vector<Vec3f> npoints(nvert);
    for (int i = 0; i < nvert; ++i)
    {
        npoints[i] = points[i] - move_min;
    }

    int nface = (int)mesh.n_faces();

    const Vec3f delta_p(info.unit[0], info.unit[1], info.unit[2]);
    const Vec3i vzero(0, 0, 0);
    const Vec3i grid_max(info.gridsize[0] - 1, info.gridsize[1] - 1, info.gridsize[2] - 1);
    size_t yz = static_cast<size_t>(info.gridsize[1]) * static_cast<size_t>(info.gridsize[2]);
    Vec3f invuint(1 / info.unit[0], 1 / info.unit[1], 1 / info.unit[2]);
    const size_t s31 = 31;
    if (morton_order)
    {

pragma omp parallel for

        for (int i = 0; i < nface; ++i) {
            auto fv = mesh.cfv_begin(FaceHandle(i));
            auto& v0 = npoints[fv->idx()]; ++fv;
            auto& v1 = npoints[fv->idx()]; ++fv;
            auto& v2 = npoints[fv->idx()];
            auto e0 = v1 - v0;
            auto e1 = v2 - v1;
            auto e2 = v0 - v2;
            auto n = (e0.cross(e1)).normalize_cond();
            if (std::fabs(n[0]) < float_error) {
                continue;
            }
            n[0] = 1 / n[0];
            Vec2f v0_yz = Vec2f(v0[1], v0[2]);
            Vec2f v1_yz = Vec2f(v1[1], v1[2]);
            Vec2f v2_yz = Vec2f(v2[1], v2[2]);
            if (!checkCCW(v0_yz, v1_yz, v2_yz))
            {
                std::swap(v1_yz, v2_yz);
            }
            Vec2f bbox_max = OpenMesh::max(v0_yz, OpenMesh::max(v1_yz, v2_yz));
            Vec2f bbox_min = OpenMesh::min(v0_yz, OpenMesh::min(v1_yz, v2_yz));
            Vec2i bbox_max_grid = Vec2i(
                (int)std::floor(bbox_max[0] * invuint[1] - 0.5f),
                (int)std::floor(bbox_max[1] * invuint[2] - 0.5f));
            Vec2i bbox_min_grid = Vec2i(
                (int)std::ceil(bbox_min[0] * invuint[1] - 0.5f),
                (int)std::ceil(bbox_min[1] * invuint[2] - 0.5f));

            for (int y = bbox_min_grid[0]; y <= bbox_max_grid[0]; ++y)
            {
                for (int z = bbox_min_grid[1]; z <= bbox_max_grid[1]; ++z)
                {
                    Vec2f point = Vec2f((y + 0.5f) * info.unit[1], (z + 0.5f) * info.unit[2]);
                    int checknum = check_point_triangle(v0_yz, v1_yz, v2_yz, point);
                    if ((checknum == 1 && TopLeftEdge(v0_yz, v1_yz)) ||
                        (checknum == 2 && TopLeftEdge(v1_yz, v2_yz)) ||
                        (checknum == 3 && TopLeftEdge(v2_yz, v0_yz)) || (checknum == 0))
                    {
                        int xmax = int(get_x_coordinate(n, v0, point) * invuint[0] - 0.5f);
                        for (int x = 0; x <= xmax; x++)
                        {
                            size_t location = mortonEncode_LUT((unsigned)x, (unsigned)y, (unsigned)z);
                            size_t int_location = location >> 5;
                            unsigned int bit_pos = size_t(31) - (location & s31);
                            unsigned int mask = 1 << bit_pos;

pragma omp atomic

                            voxel_table[int_location] ^= mask;
                        }
                    }
                }
            }
        }

    }
    else
    {

pragma omp parallel for

        for (int i = 0; i < nface; ++i)
        {
            auto fv = mesh.cfv_begin(FaceHandle(i));
            auto& v0 = npoints[fv->idx()]; ++fv;
            auto& v1 = npoints[fv->idx()]; ++fv;
            auto& v2 = npoints[fv->idx()];
            auto e0 = v1 - v0;
            auto e1 = v2 - v1;
            auto e2 = v0 - v2;
            auto n = (e0.cross(e1)).normalize_cond();
            if (std::fabs(n[0]) < float_error) {
                continue;
            }
            n[0] = 1 / n[0];
            Vec2f v0_yz = Vec2f(v0[1], v0[2]);
            Vec2f v1_yz = Vec2f(v1[1], v1[2]);
            Vec2f v2_yz = Vec2f(v2[1], v2[2]);
            if (!checkCCW(v0_yz, v1_yz, v2_yz))
            {
                std::swap(v1_yz, v2_yz);
            }
            Vec2f bbox_max = OpenMesh::max(v0_yz, OpenMesh::max(v1_yz, v2_yz));
            Vec2f bbox_min = OpenMesh::min(v0_yz, OpenMesh::min(v1_yz, v2_yz));
            Vec2i bbox_max_grid = Vec2i(std::floor(bbox_max[0] * invuint[1] - 0.5f),
                                        std::floor(bbox_max[1] * invuint[2] - 0.5f));
            Vec2i bbox_min_grid = Vec2i(std::ceil(bbox_min[0] * invuint[1] - 0.5f),
                                        std::ceil(bbox_min[1] * invuint[2] - 0.5f));
            size_t ypos = static_cast<size_t>(bbox_min_grid[0] * info.gridsize[1]);
            for (int y = bbox_min_grid[0]; y <= bbox_max_grid[0]; ++y)
            {
                size_t zyz = bbox_min_grid[1] * yz;
                for (int z = bbox_min_grid[1]; z <= bbox_max_grid[1]; ++z)
                {
                    Vec2f point((y + 0.5f) * info.unit[1], (z + 0.5f) * info.unit[2]);
                    int checknum = check_point_triangle(v0_yz, v1_yz, v2_yz, point);
                    if ((checknum == 1 && TopLeftEdge(v0_yz, v1_yz)) ||
                        (checknum == 2 && TopLeftEdge(v1_yz, v2_yz)) ||
                        (checknum == 3 && TopLeftEdge(v2_yz, v0_yz)) || 
                        (checknum == 0))
                    {
                        size_t xmax = size_t(get_x_coordinate(n, v0, point) * invuint[0] - 0.5f);
                        for (int x = 0; x <= xmax; x++)
                        {
                            size_t location = x + ypos + zyz;
                            size_t int_location = location >> 5;
                            unsigned int bit_pos = size_t(31) - (location & s31);
                            unsigned int mask = 1 << bit_pos;

pragma omp atomic

                            voxel_table[int_location] ^= mask;
                        }
                    }
                    zyz += yz;
                }
                ypos += info.gridsize[1];
            }
        }
    }
}

Just code improvement, not algorithm improvement

[Forceflow]

Could you make a PR for this? It's hard comparing code with all the formatting gone.