Gauss–Bonnet theorem简介

这里稍作说明：

代码实现

double checkG_B(CMyMesh *pMesh)
{
    double GB = 0;
    for (CMyMesh::MeshVertexIterator viter(pMesh); !viter.end(); ++viter)
    {
        CMyVertex *v = *viter;
        if(v->boundary()) {     //pi
            std::vector<CPoint> vEdge;
            double totalAngle = 0;
            for(CMyMesh::VertexVertexIterator vviter(v); !vviter.end(); ++vviter)
            {
                CMyVertex *pV = *vviter;
                CPoint p = pV->point();
                CPoint e = p - v->point();   //e1
                vEdge.push_back(e);
            }
            for(size_t i = 1; i < vEdge.size(); ++i) {
                CPoint pa = vEdge.at(i - 1),
                       pb = vEdge.at(i);
                double cos_ = pa * pb / mod(pa) / mod(pb);
                totalAngle += acos(cos_);
            }
            GB += (PI - totalAngle);
        }
        if(!v->boundary()) {     //2 * pi
            std::vector<CPoint> vEdge;
            double totalAngle = 0;
            for(CMyMesh::VertexVertexIterator vviter(v); !vviter.end(); ++vviter)
            {
                CMyVertex *pV = *vviter;
                CPoint p = pV->point();
                CPoint e = p - v->point();   //e1
                vEdge.push_back(e);
            }
            for(size_t i = 1; i < vEdge.size(); ++i) {
                CPoint pa = vEdge.at(i - 1),
                       pb = vEdge.at(i);
                double cos_ = pa * pb / mod(pa) / mod(pb);
                totalAngle += acos(cos_);
            }
            //below is the important codes without witch the result will be wrong
            //the angle between the fist vector and the last vector
            CPoint p_first = vEdge.at(0), p_last = vEdge.at(vEdge.size() - 1);
            double cos_ = p_first * p_last / mod(p_first) / mod(p_last);
            GB += (2 * PI - totalAngle);
        }

    }
    return GB;
}