trimesh.cpp

#include "trimesh.h"

// needed for  implementation
#include <cassert>
#include <set>
#include <iostream>
#include <algorithm>

namespace trimesh {
    typedef std::map<std::pair<trimesh_t::index_t, trimesh_t::index_t>, trimesh_t::index_t> directed_edge2index_map_t;

    trimesh_t::index_t directed_edge2face_index(const directed_edge2index_map_t& de2fi, trimesh_t::index_t vertex_i, trimesh_t::index_t vertex_j) {
        assert(!de2fi.empty());
        
        directed_edge2index_map_t::const_iterator it = de2fi.find(std::make_pair(vertex_i, vertex_j));
        
        // If no such directed edge exists, then there's no such face in the mesh.
        // The edge must be a boundary edge.
        // In this case, the reverse orientation edge must have a face.
        if (it == de2fi.end()) {
            assert(de2fi.find(std::make_pair(vertex_j, vertex_i)) != de2fi.end());
            return -1;
        }
        
        return it->second;
    }

    void trimesh_t::build(const unsigned long num_vertices, const unsigned long num_triangles, const triangle_t* triangles, const unsigned long num_edges, const edge_t* edges) {
        /*
        Generates all half edge data structures for the mesh given by its vertices 'self.vs'
        and faces 'self.faces'.
        
        Python version used heavily
        */
        
        assert(triangles);
        assert(edges);
        
        directed_edge2index_map_t de2fi;
        for (int fi = 0; fi <num_triangles; ++fi) {
            const triangle_t& tri = triangles[fi];
            de2fi[std::make_pair(tri.v[0], tri.v[1])] = fi;
            de2fi[std::make_pair(tri.v[1], tri.v[2])] = fi;
            de2fi[std::make_pair(tri.v[2], tri.v[0])] = fi;
        }
        
        clear();
        _vertex_halfedges.resize(num_vertices, -1);
        _face_halfedges.resize(num_triangles, -1);
        _edge_halfedges.resize(num_edges, -1);
        _halfedges.reserve(num_edges*2);
        
        for (int ei = 0; ei <num_edges; ++ei) {
            const edge_t& edge = edges[ei];
            
            // Add the halfedge_t structures to the end of the list.
            const index_t he0index = _halfedges.size();
            _halfedges.push_back(halfedge_t());
            halfedge_t& he0 = _halfedges.back();
            
            const index_t he1index = _halfedges.size();
            _halfedges.push_back(halfedge_t());
            halfedge_t& he1 = _halfedges.back();
            
            // The face will be -1 if it is a boundary half-edge.
            he0.face = directed_edge2face_index(de2fi, edge.v[0], edge.v[1]);
            he0.to_vertex = edge.v[1];
            he0.edge = ei;
            
            // The face will be -1 if it is a boundary half-edge.
            he1.face = directed_edge2face_index(de2fi, edge.v[1], edge.v[0]);
            he1.to_vertex = edge.v[0];
            he1.edge = ei;
            
            // Store the opposite half-edge index.
            he0.opposite_he = he1index;
            he1.opposite_he = he0index;
            
            // Also store the index in our _directed_edge2he_index map.
            assert(_directed_edge2he_index.find(std::make_pair(edge.v[0], edge.v[1])) == _directed_edge2he_index.end());
            assert(_directed_edge2he_index.find(std::make_pair(edge.v[1], edge.v[0])) == _directed_edge2he_index.end());
            _directed_edge2he_index[std::make_pair(edge.v[0], edge.v[1])] = he0index;
            _directed_edge2he_index[std::make_pair(edge.v[1], edge.v[0])] = he1index;
            
            // If the vertex pointed to by a half-edge doesn't yet have an out-going
            // halfedge, store the opposite halfedge.
            // Also, if the vertex is a boundary vertex, make sure its
            // out-going halfedge is a boundary halfedge.
            // NOTE: Halfedge data structure can't properly handle butterfly vertices.
            //       If the mesh has butterfly vertices, there will be multiple outgoing
            //       boundary halfedges.  Because we have to pick one as the vertex's outgoing
            //       halfedge, we can't iterate over all neighbors, only a single wing of the
            //       butterfly.
            if (_vertex_halfedges[he0.to_vertex] == -1 || -1 == he1.face) {
                _vertex_halfedges[he0.to_vertex] = he0.opposite_he;
            }
            if (_vertex_halfedges[he1.to_vertex] == -1 || -1 == he0.face) {
                _vertex_halfedges[he1.to_vertex] = he1.opposite_he;
            }
            
            // If the face pointed to by a half-edge doesn't yet have a
            // halfedge pointing to it, store the halfedge.
            if (-1 != he0.face && _face_halfedges[he0.face] == -1) {
                _face_halfedges[he0.face] = he0index;
            }
            if (-1 != he1.face && _face_halfedges[he1.face] == -1) {
                _face_halfedges[he1.face] = he1index;
            }
            
            // Store one of the half-edges for  the edge.
            assert(_edge_halfedges[ei] == -1);
            _edge_halfedges[ei] = he0index;
        }
        
        // Now that all the half-edges are created, set the remaining next_he field.
        // We can't yet handle boundary halfedges, so store them for  later.
        std::vector<index_t> boundary_heis;
        for (int hei = 0; hei <_halfedges.size(); ++hei) {
            halfedge_t& he = _halfedges.at(hei);
            // Store boundary halfedges for  later.
            if (-1 == he.face) {
                boundary_heis.push_back(hei);
                continue;
            }
            
            const triangle_t& face = triangles[he.face];
            const index_t i = he.to_vertex;
            
            index_t j = -1;
            if (face.v[0] == i) j = face.v[1];
            else if (face.v[1] == i) j = face.v[2];
            else if (face.v[2] == i) j = face.v[0];
            assert(-1 != j);
            
            he.next_he = _directed_edge2he_index[std::make_pair(i,j)];
        }
        
        // Make a map from vertices to boundary halfedges (indices) originating from them.
        // NOTE: There will only be multiple originating boundary halfedges at butterfly vertices.
        std::map<index_t, std::set<index_t>> vertex2outgoing_boundary_hei;
        for (std::vector<index_t>::const_iterator hei = boundary_heis.begin(); hei != boundary_heis.end(); ++hei) {
            const index_t originating_vertex = _halfedges[_halfedges[*hei].opposite_he].to_vertex;
            vertex2outgoing_boundary_hei[originating_vertex].insert(*hei);

            if (vertex2outgoing_boundary_hei[originating_vertex].size()> 1) {
                std::cerr <<"Butterfly vertex encountered.\n";
            }
        }
        
        // for  each boundary halfedge, make its next_he one of the boundary halfedges
        // originating at its to_vertex.
        for (std::vector<index_t>::const_iterator hei = boundary_heis.begin(); hei != boundary_heis.end(); ++hei) {
            halfedge_t& he = _halfedges[*hei];
            
            std::set<index_t>& outgoing = vertex2outgoing_boundary_hei[he.to_vertex];

            if (!outgoing.empty()) {
                std::set<index_t>::iterator outgoing_hei = outgoing.begin();
                he.next_he = *outgoing_hei;
                
                outgoing.erase(outgoing_hei);
            }
        }
    }

    std::vector<trimesh_t::index_t> trimesh_t::boundary_vertices() const {
        // Returns a list of the vertex indices on the boundary.

        std::set<index_t> result;
        for (int hei = 0; hei <_halfedges.size(); ++hei)
        {
            const halfedge_t& he = _halfedges[hei];
            
            if (-1 == he.face)
            {
                // result.extend(self.he_index2directed_edge(hei))
                result.insert(he.to_vertex);
                result.insert(_halfedges[he.opposite_he].to_vertex);
            }
        }
        
        return std::vector<index_t>(result.begin(), result.end());
    }

    std::vector<std::pair<trimesh_t::index_t, trimesh_t::index_t>> trimesh_t::boundary_edges() const {  
        // Returns a list of undirected boundary edges (i,j).  If (i,j) is in the result, (j,i) will not be.
        std::vector<std::pair<index_t, index_t>> result;

        for (int hei = 0; hei <_halfedges.size(); ++hei) {
            const halfedge_t& he = _halfedges[hei];
            
            if (-1 == he.face)
            {
                result.push_back(he_index2directed_edge(hei));
            }
        }
        
        return result;
    }

    void unordered_edges_from_triangles(const unsigned long num_triangles, const triangle_t* triangles, std::vector<edge_t>& edges_out) {
        typedef std::set<std::pair<index_t, index_t>> edge_set_t;
        edge_set_t edges;
        for (int t = 0; t <num_triangles; ++t) {
            edges.insert(std::make_pair(std::min(triangles[t].i(), triangles[t].j()), std::max(triangles[t].i(), triangles[t].j())));
            edges.insert(std::make_pair(std::min(triangles[t].j(), triangles[t].k()), std::max(triangles[t].j(), triangles[t].k())));
            edges.insert(std::make_pair(std::min(triangles[t].k(), triangles[t].i()), std::max(triangles[t].k(), triangles[t].i())));
        }
        
        edges_out.resize(edges.size());
        int e = 0;
        for (edge_set_t::const_iterator it = edges.begin(); it != edges.end(); ++it, ++e) {
            edges_out.at(e).start() = it->first;
            edges_out.at(e).end() = it->second;
        }
    }
}