libigl / libigl

Simple MPL-2.0-licensed C++ geometry processing library.
http://libigl.github.io/libigl/
GNU General Public License v3.0
4.55k stars 1.13k forks source link

[Question] cut_mesh and cut_to_disk #1328

Open xarthurx opened 4 years ago

xarthurx commented 4 years ago

The cut_to_disk function mentioned that the cut mesh can be materialised by the cut_mesh function.

While cut_to_disk function generate cuts as a std::vector<std::vector<int>> type, the cut_mesh function take cuts variable as an Eigen type.

Would it make sense to add an overload function of either cut_to_disk or cut_mesh function so that the user don't need to manually translate between the types?

xarthurx commented 4 years ago

The above functionality can be achieved similar like:

            // cuts -> fCuts: translate the vertex-style cuts into the F oriented cuts representation
            MatrixXi fCuts;
            fCuts.resizeLike(viewer.data().F);
            fCuts.setZero();

            std::vector<std::vector<int>> VF, VFi;
            igl::vertex_triangle_adjacency(viewer.data().V, viewer.data().F, VF, VFi);
            for (auto i : cuts[0]) {
              for (size_t f = 0; f < VF[i].size(); f++) {
                fCuts.row(VF[i][f])[VFi[i][f]] = 1; // VF[i][f] -- the face that the specific v is in; VFi[i][f] -- the idx of the vert in that face
              }
            }

Additionally, since cut_mesh doesn't handle non-edge-manifold cuts, while cut_to_disk can output that, it will also be helpful to add possibility for cut_mesh to cut in that way.

e.g. I have a topology-disk mesh, with an edge trail, I'd like to split the mesh into two parts.

  1. Use cut_to_disk to output the boundary of the the patch;
  2. use cut_mesh to extract the mesh.

In general, I think libigl is lacking of user-friendly function for CAD-like operations...

xarthurx commented 4 years ago

For now, I'm still trying to figure out how the above "looks-simple" task can be achieved from libigl...

xarthurx commented 4 years ago

OK, I think I made a mistake, the cut given to cut_mesh should be #F x 3 of edge idx, not vert idx. But still confused about how are these cuts generated...

The documentation of this part is really not very clear.

BruegelN commented 4 years ago

Maybe have a look at the assert_is_disk()in the unit test of cut_to_disk on how to get from std::vector<std::vector<int>> cuts to Eigen::MatrixXi cut_mask(F.rows(), 3).

xarthurx commented 4 years ago

@BruegelN Yeah, this also works, thanks.

evouga commented 4 years ago

Late to the party here, but here is the code I wrote for performing the cut that is found using cut_to_disk. I didn't commit it into libigl since libigl already had an implementation of mesh-cutting, but maybe it's useful to you (or anyone else with the same question):

static void cutMesh(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
    // list of cuts, each of which is a list (in order) of vertex indices of one cut.
    // Cuts can be closed loops (in which case the last vertex index should equal the
    // first) or open (in which case the two endpoint vertices should be distinct).
    // Multiple cuts can cross but there may be strange behavior if cuts share endpoint
    // vertices, or are non-edge-disjoint.
    const std::vector<std::vector<int> > &cuts,
    // new vertices and faces
    // **DO NOT ALIAS V OR F!**
    Eigen::MatrixXd &newV,
    Eigen::MatrixXi &newF
)
{
    int ncuts = (int)cuts.size();

    // junction vertices that lie on multiple cuts
    std::set<int> junctions;
    std::set<int> seenverts;
    for (int i = 0; i < ncuts; i++)
    {
        std::set<int> seenincut;
        for (int j = 0; j < cuts[i].size(); j++)
        {
            if (seenverts.count(cuts[i][j]))
                junctions.insert(cuts[i][j]);
            seenincut.insert(cuts[i][j]);
        }
        for (int v : seenincut)
            seenverts.insert(v);
    }

    // "interior" cut vertices: vertices that are part of a cut but not a cut endpoint
    // or junction vertex
    std::vector<std::set<int> > cutints;
    cutints.resize(ncuts);
    for (int i = 0; i < ncuts; i++)
    {
        if (cuts[i].empty())
            continue;
        if (cuts[i].front() == cuts[i].back())
        {
            // closed loop
            for (int v : cuts[i])
            {
                if(!junctions.count(v))
                    cutints[i].insert(v);                
            }
        }
        else
        {
            // open cut
            for (int j = 1; j < cuts[i].size() - 1; j++)
            {
                if(!junctions.count(cuts[i][j]))
                    cutints[i].insert(cuts[i][j]);                
            }
        }
    }

    struct edge
    {
        std::pair<int, int> verts;
        edge(int v1, int v2)
        {
            verts.first = std::min(v1, v2);
            verts.second = std::max(v1, v2);
        }

        bool operator<(const edge &other) const
        {
            return verts < other.verts;
        }
    };

    // maps each edge to incident triangles
    std::map<edge, std::vector<int> > edgeTriangles;
    for (int i = 0; i < F.rows(); i++)
    {
        for (int j = 0; j < 3; j++)
        {
            edge e(F(i, j), F(i, (j + 1) % 3));
            edgeTriangles[e].push_back(i);
        }
    }

    // have we visited this face yet?
    bool *visited = new bool[F.rows()];

    // edges that form part of a cut
    std::set<edge> forbidden;
    for (int i = 0; i < ncuts; i++)
    {        
        for (int j = 0; j < (int)cuts[i].size() - 1; j++)
        {
            // works for both open and closed curves
            edge e(cuts[i][j], cuts[i][j + 1]);
            forbidden.insert(e);
        }
    }

    // connected components of faces adjacent to the cuts
    std::vector<std::vector<std::vector<int> > > components;
    components.resize(ncuts);

    // for each cut
    for (int cut = 0; cut < ncuts; cut++)
    {
        for (int i = 0; i < (int)F.rows(); i++)
            visited[i] = false;

        // find a face we haven't visited yet
        for (int i = 0; i < F.rows(); i++)
        {
            if (visited[i]) continue;
            bool found = false;
            for (int j = 0; j < 3; j++)
            {
                if (cutints[cut].count(F(i, j)))
                {
                    found = true;
                }
            }

            if (found)
            {
                // run a BFS along the cut edges, accumulating one connected component
                // cross only edges that contain a vertex in cutints[cut], but are not forbidden
                std::deque<int> q;
                std::vector<int> component;
                q.push_back(i);
                while (!q.empty())
                {
                    int next = q.front();
                    q.pop_front();
                    if (visited[next])
                        continue;
                    visited[next] = true;
                    component.push_back(next);
                    for (int j = 0; j < 3; j++)
                    {
                        int v1 = F(next, j);
                        int v2 = F(next, (j + 1) % 3);
                        edge e(v1, v2);
                        if (cutints[cut].count(v1) == 0 && cutints[cut].count(v2) == 0)
                        {
                            continue;
                        }
                        if (forbidden.count(e))
                            continue;
                        for (int nb : edgeTriangles[e])
                        {
                            if (!visited[nb])
                            {
                                q.push_back(nb);
                            }
                        }
                    } // end BFS
                }
                components[cut].push_back(component);
            } // end if found
        } // end loop over all faces
    } // end loop over cuts

    std::map<int, std::vector<std::vector<int> > > junctcomponents;

    // for each junction
    for (int junc : junctions)
    {
        for (int i = 0; i < (int)F.rows(); i++)
            visited[i] = false;

        // find a face we haven't visited yet
        for (int i = 0; i < F.rows(); i++)
        {
            if (visited[i]) continue;
            bool found = false;
            for (int j = 0; j < 3; j++)
            {
                if (junc == F(i, j))
                {
                    found = true;
                }
            }

            if (found)
            {
                // run a BFS along the cut edges, accumulating one connected component
                // cross only edges that contain the junction, but are not forbidden
                std::deque<int> q;
                std::vector<int> component;
                q.push_back(i);
                while (!q.empty())
                {
                    int next = q.front();
                    q.pop_front();
                    if (visited[next])
                        continue;
                    visited[next] = true;
                    component.push_back(next);
                    for (int j = 0; j < 3; j++)
                    {
                        int v1 = F(next, j);
                        int v2 = F(next, (j + 1) % 3);
                        edge e(v1, v2);
                        if (v1 != junc && v2 != junc)
                        {
                            continue;
                        }
                        if (forbidden.count(e))
                            continue;
                        for (int nb : edgeTriangles[e])
                        {
                            if (!visited[nb])
                            {
                                q.push_back(nb);
                            }
                        }
                    } // end BFS
                }
                junctcomponents[junc].push_back(component);
            } // end if found
        } // end loop over all faces
    } // end loop over cuts

    int vertstoadd = 0;
    // create a copy of each vertex for each component of each cut
    for (int i = 0; i < ncuts; i++)
    {
        vertstoadd += components[i].size() * cutints[i].size();
    }
    // create a copy of each junction point for each component of each junction
    for (int v : junctions)
    {
        vertstoadd += junctcomponents[v].size();
    }

    Eigen::MatrixXd augV(V.rows() + vertstoadd, 3);
    for (int i = 0; i < V.rows(); i++)
    {
        augV.row(i) = V.row(i);
    }

    // create new faces
    Eigen::MatrixXi augF = F;

    // duplicate vertices and reindex faces
    int idx = V.rows();
    for (int cut = 0; cut < ncuts; cut++)
    {
        for (int i = 0; i < components[cut].size(); i++)
        {
            // duplicate vertices
            std::map<int, int> idxmap;
            for (int v : cutints[cut])
            {
                idxmap[v] = idx;
                augV.row(idx) = V.row(v);
                idx++;
            }
            for (int f : components[cut][i])
            {
                for (int j = 0; j < 3; j++)
                {
                    int v = augF(f, j);
                    if (cutints[cut].count(v))
                        augF(f, j) = idxmap[v];
                }
            }
        }        
    }

    for (int junc : junctions)
    {
        for (int i = 0; i < junctcomponents[junc].size(); i++)
        {

            augV.row(idx) = V.row(junc);
            int newidx = idx;
            idx++;

            for (int f : junctcomponents[junc][i])
            {
                for (int j = 0; j < 3; j++)
                {
                    int v = augF(f, j);
                    if (v == junc)
                        augF(f, j) = newidx;
                }
            }
        }
    }

    assert(idx == augV.rows());

    // some duplicated vertices will not have been used
    Eigen::MatrixXi I;
    igl::remove_unreferenced(augV, augF, newV, newF, I);
    delete[] visited;
}