C++ SurfaceMesh::get_simplex_up方法代码示例

Vector getNormal(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<3> faceID)
{
    Vector norm;
    auto   name = mesh.get_name(faceID);

    // Get the three vertices
    auto a = *mesh.get_simplex_up({name[0]});
    auto b = *mesh.get_simplex_up({name[1]});
    auto c = *mesh.get_simplex_up({name[2]});

    if ((*faceID).orientation == 1)
    {
        norm = cross(c-b, a-b);
    }
    else if ((*faceID).orientation == -1)
    {
        norm = cross(a-b, c-b);
    }
    else
    {
        // std::cerr << "ERROR(getNormal): Orientation undefined, cannot compute "
                  // << "normal. Did you call compute_orientation()?" << std::endl;
        throw std::runtime_error("ERROR(getNormal): Orientation undefined, cannot compute normal. Did you call compute_orientation()?");
    }
    return norm;
}

double getArea(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<3> faceID)
{
    auto name = mesh.get_name(faceID);
    auto a = *mesh.get_simplex_up({name[0]});
    auto b = *mesh.get_simplex_up({name[1]});
    auto c = *mesh.get_simplex_up({name[2]});
    return getArea(a, b, c);
}

/**
 * http://research.microsoft.com/en-us/um/people/chazhang/publications/icip01_ChaZhang.pdf
 * http://chenlab.ece.cornell.edu/Publication/Cha/icip01_Cha.pdf
 */
double getVolume(const SurfaceMesh &mesh)
{
    bool orientError = false;
    double volume = 0;
    for (auto faceID : mesh.get_level_id<3>())
    {
        auto   name = mesh.get_name(faceID);
        auto   a = (*mesh.get_simplex_up({name[0]})).position;
        auto   b = (*mesh.get_simplex_up({name[1]})).position;
        auto   c = (*mesh.get_simplex_up({name[2]})).position;

        Vector norm;
        double tmp;
        if ((*faceID).orientation == 1)
        {
            // a->b->c
            norm = cross(b, c);
            tmp  = dot(a, norm);

        }
        else if ((*faceID).orientation == -1)
        {
            // c->b->a
            norm = cross(b, a);
            tmp  = dot(c, norm);
        }
        else
        {
            orientError = true;
        }

        // * Probably less efficient way #1
        // tmp = dot(a, getNormal(mesh, faceID));

        // * Less efficient way #2
        // norm = getNormalFromTangent(getTangent(mesh, faceID));
        // auto wedge = a^b^c;
        // tmp = std::sqrt(wedge|wedge);
        // auto sgn = dot((a+b+c)/3, norm);
        // if(sgn <= 0) {
        //     tmp = -1*tmp;
        // }
        volume += tmp;
    }
    if(orientError){
        std::cerr << "ERROR getVolume(): Orientation undefined for one or more "
                  << "simplices. Did you call compute_orientation()?" << std::endl;
    }
    return volume/6;
}

tensor<double, 3, 2> getTangent(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<3> faceID)
{
    auto cover = mesh.get_name(faceID);
    auto vertexID = mesh.get_simplex_up({cover[0]});
    std::set<SurfaceMesh::KeyType> next(std::begin(cover)+1, std::end(cover));
    return surfacemesh_detail::getTangentF(mesh, (*vertexID).position, vertexID, next);
}

std::tuple<double, double, int, int> getMinMaxAngles(
    const SurfaceMesh &mesh,
    int maxMinAngle,
    int minMaxAngle)
{
    double minAngle = 360;
    double maxAngle = 0;
    int small = 0;
    int large = 0;

    // for each triangle
    for(auto nid : mesh.get_level_id<3>()){
        auto name = mesh.get_name(nid);
        auto a = *mesh.get_simplex_up({name[0]});
        auto b = *mesh.get_simplex_up({name[1]});
        auto c = *mesh.get_simplex_up({name[2]});
        std::array<double, 3> angles;
        try{
            angles[0] = angle(a,b,c);
            angles[1] = angle(b,a,c);
            angles[2] = angle(a,c,b);
        }
        catch (std::runtime_error& e){
            throw std::runtime_error("ERROR(getMinMaxAngles): Cannot compute angles of face with zero area.");
        }

        for(double angle : angles){
            if (angle < minAngle){
                minAngle = angle;
            }
            if (angle > maxAngle){
                maxAngle = angle;
            }
            if (angle < maxMinAngle){
                ++small;
            }
            if (angle > minMaxAngle){
                ++large;
            }
        }
    }
    return std::make_tuple(minAngle, maxAngle, small, large);
}

/**
 * @brief      Refine the mesh by quadrisection.
 *
 * Note that this function will delete all stored data on edges and faces. But
 * this can be easily fixed.
 *
 * @param      mesh  The mesh
 */
std::unique_ptr<SurfaceMesh> refineMesh(const SurfaceMesh &mesh){
    std::unique_ptr<SurfaceMesh> refinedMesh(new SurfaceMesh);

    // Copy over vertices to refinedMesh
    for (auto vertex : mesh.get_level_id<1>()){
        auto key = mesh.get_name(vertex);
        refinedMesh->insert(key, *vertex);
    }

    // Split edges and generate a map of names before to after
    std::map<std::array<int,2>, int> edgeMap;

    for(auto edge : mesh.get_level_id<2>()){
        auto edgeName = mesh.get_name(edge);
        auto v1 = *mesh.get_simplex_up({edgeName[0]});
        auto v2 = *mesh.get_simplex_up({edgeName[1]});

        auto newVertex = refinedMesh->add_vertex(Vertex(0.5*(v1+v2)));
        edgeMap.emplace(std::make_pair(edgeName, newVertex));
    }

    // Connect edges and face and copy data
    for(auto face : mesh.get_level_id<3>()){
        auto name = mesh.get_name(face);
        int a, b, c;

        // Skip checking if found
        auto it = edgeMap.find({name[0],name[1]});
        a = it->second;

        it = edgeMap.find({name[1],name[2]});
        b = it->second;

        it = edgeMap.find({name[0],name[2]});
        c = it->second;

        refinedMesh->insert({name[0], a});
        refinedMesh->insert({a, name[1]});

        refinedMesh->insert({name[1], b});
        refinedMesh->insert({b, name[2]});

        refinedMesh->insert({name[0], c});
        refinedMesh->insert({c, name[2]});

        refinedMesh->insert({a,b,c});
        refinedMesh->insert({name[0],a,c}, *face);
        refinedMesh->insert({name[1],a,b}, *face);
        refinedMesh->insert({name[2],b,c}, *face);
    }
    return refinedMesh;
}

void printQualityInfo(const std::string &filename, const SurfaceMesh &mesh){

    std::stringstream anglefilename;
    anglefilename << filename << ".angle";

    std::ofstream angleOut(anglefilename.str());
    if(!angleOut.is_open())
    {
        std::stringstream ss;
        ss << "Couldn't open file " << filename << ".angle";
        throw std::runtime_error(ss.str());
    }

    std::stringstream areafilename;
    areafilename << filename << ".area";
    std::ofstream areaOut(areafilename.str());
    if(!areaOut.is_open())
    {
        std::stringstream ss;
        ss << "Couldn't open file " << filename << ".area";
        throw std::runtime_error(ss.str());
    }

    for(auto faceID : mesh.get_level_id<3>()){
        auto name = mesh.get_name(faceID);
        auto a = *mesh.get_simplex_up({name[0]});
        auto b = *mesh.get_simplex_up({name[1]});
        auto c = *mesh.get_simplex_up({name[2]});

        areaOut << getArea(a,b,c) << "\n";
        // Print angles in degrees
        angleOut    << angle(a,b,c) << "\n"
                    << angle(b,a,c) << "\n"
                    << angle(a,c,b) << "\n";
    }
    areaOut.close();
    angleOut.close();
}