STEP BY STEP

对模型进行谱分析计算特征向量

发表于 分类于 阅读次数： Valine：
本文字数： 3.3k 阅读时长 ≈ 3 分钟

一个三维模型，可以看做是一个无向图。顶点是图的点，边可以作为图的边，那么可以得到一个邻接矩阵，该矩阵叫做拉普拉斯矩阵，计算这个邻接矩阵的特征向量就是所谓的谱分析。

那么问题来了。如何定义边的权重？额，当然有不同的方法，看你需要分析的特性是什么。一般来说是边的两个对角的cotangle之和。另外，矩阵对角线元素是该行其它元素之和再取负。下一步就是计算该矩阵的特征向量了。如果计算了？我前面写过两篇文章讲述如何使用MATLAB计算稀疏矩阵的特征向量，现在我就是这么做的。首先，拉普拉斯矩阵是一个对称矩阵，而且是一个大型的稀疏矩阵。因此，只能使用稀疏矩阵来存储数据了。而且对称矩阵能保证一些特殊的性质，比如第一个特征向量是k * (1,1,1,1…)。第二个特征向量分布得比较有规律，但是通过实验发现，也许第二个也没有规律，可能是第三或者第四个，因为要利用特征向量的分布，需要综合多个特征向量。

下面是我计算拉普拉斯矩阵的代码，需要考虑角度大于90的情况，以及除0等等。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
//a和b夹角的cot
    double CSEMesh::GetCotangent(double a, double b, double c)
    {
        if (a*b < 1e-8)
        {
            return 0;
        }
        else
        {
            double cosc = (a*a+b*b-c*c)/(2*a*b);

            if (cosc < -1)
                cosc = -1;
            if (cosc > 1)
                cosc = 1;

            double angle = acos(cosc);

            double tan_angle = tan(angle);

            if (fabs(tan_angle) < 1e-8)
                return 1 / 1e-8;
            else
                return 0.5 / tan_angle;
        }
    }

    void CSEMesh::ComputeLaplaceWeight()
    {
        for (Enriched_Model::Edge_iterator pEdge = m_pModel->edges_begin();
            pEdge != m_pModel->edges_end();
            pEdge++)
        {
            double cotA = 1.0, cotB = 1.0;
            double a, b, c;

            Enriched_Model::Edge_iterator pEdgeTmp = pEdge;
            if (!pEdgeTmp->is_border())
            {
                c = pEdgeTmp->length();
                b = pEdgeTmp->next()->length();
                a = pEdgeTmp->next()->next()->length();
                cotA = GetCotangent(a, b, c);
            }

            pEdgeTmp = pEdge->opposite();
            if (!pEdgeTmp->is_border())
            {
                c = pEdgeTmp->length();
                b = pEdgeTmp->next()->length();
                a = pEdgeTmp->next()->next()->length();
                cotB = GetCotangent(a, b, c);
            }

            if(cotA + cotB < 1e-7)
            {
                pEdge->weight() = 1e-7;
                pEdge->opposite()->weight() = 1e-7;
            }
            else
            {
                pEdge->weight() = cotA + cotB;
                pEdge->opposite()->weight() = cotA + cotB;
            }
        }
    }

接下来是计算特征向量的代码，因为混合着cgal等，不是很好分离，但是可以用来参考下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
void CSEMesh::ComputeEigen(int type)
{
        static vector<MatlabSparseInfor> siVec;

        int nRow = m_pModel->size_of_vertices();
        int nColum = nRow;
        m_eigenValues.resize(EIGEN_VECTOR_NUM);
        for (int i = 0; i < EIGEN_VECTOR_NUM; ++i)
        {
            m_eigenVectors[type][i].resize(nRow);
        }

        int nNum = m_pModel->size_of_halfedges();
        siVec.resize(nNum + nRow);

        int i = 0;
        for (Enriched_Model::Edge_iterator pEdge = m_pModel->edges_begin();
            pEdge != m_pModel->edges_end();
            pEdge++)
        {
            int beg = pEdge->opposite()->vertex()->tag();
            int end = pEdge->vertex()->tag();

            siVec[i].r = beg;
            siVec[i].c = end;
            siVec[i].v = -pEdge->weight();
            ++i;

            siVec[i].r = end;
            siVec[i].c = beg;
            siVec[i].v = -pEdge->weight();
            ++i;
        }
        assert(i == nNum);

        for(Enriched_Model::Vertex_iterator pVertex = m_pModel->vertices_begin(); 
            pVertex != m_pModel->vertices_end(); pVertex++)
        {
            Enriched_Model::Vertex::Halfedge_around_vertex_circulator pHalfEdge = pVertex->vertex_begin();
            Enriched_Model::Vertex::Halfedge_around_vertex_circulator d = pHalfEdge;
            double sum = 0.0;

            CGAL_For_all(pHalfEdge, d)
            {
                sum += pHalfEdge->weight();
            }

            if (fabs(sum) < 1e-7)
            {
                sum = 0.0;
            }

            siVec[i].r = pVertex->tag();
            siVec[i].c = pVertex->tag();
            siVec[i].v = sum;
            ++i;
        }
        assert(i == nNum + nRow);
        std::sort(siVec.begin(), siVec.end());

        CMatlabSparseMatrix sm(siVec[0], nNum + nRow, nRow, nColum);
        int nEigenNum = std::min(nRow, (int)EIGEN_VECTOR_NUM);
        sm.GetEigens(m_eigenValues, m_eigenVectors[type], nEigenNum);
}

关于如何实现MATLAB稀疏矩阵类，可以参考我其它的两篇文章。下面是第二个特征向量的分布效果，能够观察到是一个非常有规律的标量场。

相关文章