STEP BY STEP

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

上一篇文章说了提取拉普拉斯矩阵的第二个（或者其它）特征向量能得到一个有规律分布的标量场。现在讲的是如何在这样一个规律分布的标量场上提取出等值线。注意，这是在一个三角网格上面提取等值线，并不是真正的曲面。提取出等值线后，可以做一些相关的应用。

那么，如何提取等值线了？

首先，如果只是绘制出等值线，可以简单的设置纹理贴图值（比如间隔多少行，设置纹理值为白色），标量场对应一维纹理坐标。更详细的说明，可以参考我的一篇博文：在三维模型上可视化标量场及注意方面。

现在考虑如何计算等值线上的间隔点。先看单个的三角面，无论如何可以得到一条等值线的。最简单的思想，就是把单个三角面上的等值线连接起来。

给定一个三角面，如何计算从该三角面出发的一条等值线了？

第一步，求出三个顶点的标量值顺序，那么标量范围在[fMin,fMax]之间。可以任选区间内的一个值，比如中间值作为该等值线的标量值。

第二步，得到等值线标量值在[fMin,fMax]所在边pEdge的位置ptBeg，作为起始点，加入点集合isoLine。

第三步，在pEdge所在的三角面上，找pEdge相邻的两条边的等值点(ptBeg的标量值相等的点)ptNext。如果存在这样的点，那么修改ptBeg = ptNext，加入点集合isoLine。同时，修改pEdge为对应的相邻边pEdge->next()或者是pEdge->next()->next()。如果不存在ptNext，进入第四步，同时设置不成环标志。否则，pEdge = pEdge->opposite()，(效果是进入相邻面，再循环第三步），如果发现ptNext等于isoLine的第一个点，那么说明已经成环，进入第四步，同时设置成环标志。

第四步，返回等值线，以及等值线是否成环状。

相关代码如下，里面使用了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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
bool CSEMesh::ComputeIsoLine(std::vector<double> scalarField, Enriched_Model::Facet_iterator pBegFacet,
    std::vector<bool> bHasIsoLine, IsoLine isoline)
{
    //该面已经有通过的isoline
    if (bHasIsoLine[pBegFacet->tag()] == true)
    {
        return false;
    }

    isoline.pts.clear();
    std::vector<int> visitFaces;//这次访问过的face

    Enriched_Model::Halfedge_around_facet_circulator pHalfedge = pBegFacet->facet_begin();

    int vertexs[3];
    int i = 0;
    do
    {
        vertexs[i++] = pHalfedge->vertex()->tag();
    } while (++pHalfedge != pBegFacet->facet_begin());

    struct ScalarInfor
    {
        int index;
        double value;
        bool operator < (const ScalarInfor si) const
        {
            return value < si.value;
        }
        ScalarInfor(int i, double v) : index(i), value(v) {}
    };
    std::vector<ScalarInfor> vecSi;
    for (int i = 0; i < 3; ++i)
    {
        vecSi.push_back(ScalarInfor(vertexs[i], scalarField[vertexs[i]]));
    }
    std::sort(vecSi.begin(), vecSi.end());

    int nMinIndex = vecSi[0].index;
    int nMaxIndex = vecSi[2].index;

    Enriched_Model::Edge_iterator pEdge;
    pHalfedge = pBegFacet->facet_begin();
    do
    {
        if (nMinIndex == pHalfedge->vertex()->tag()  nMaxIndex == pHalfedge->opposite()->vertex()->tag()
        ||  nMaxIndex == pHalfedge->vertex()->tag()  nMinIndex == pHalfedge->opposite()->vertex()->tag())
        {
            pEdge = pHalfedge;
            break;
        }
    } while (++pHalfedge != pBegFacet->facet_begin());

    Enriched_Model::Point pt;
    isoline.isoValue = (scalarField[nMinIndex] + scalarField[nMaxIndex]) / 2.0;//当前等值线的值

    Enriched_Model::Point ptBeg = pEdge->opposite()->vertex()->point();
    int iBeg = pEdge->opposite()->vertex()->tag();

    Enriched_Model::Point ptEnd = pEdge->vertex()->point();
    int iEnd = pEdge->vertex()->tag();
    GetInterpolationPt(ptBeg, ptEnd, scalarField[iBeg], scalarField[iEnd], isoline.isoValue, pt);//得到对边的插值点作为终点
    isoline.pts.push_back(pt);

    //下一个三角面是pEdge所在的相邻三角面
    bool bCircle = true;
    while (1)
    {
        if (pEdge->is_border())
        {
            bCircle = false;
            break;
        }

        int iBeg = pEdge->opposite()->vertex()->tag();
        int iEnd = pEdge->vertex()->tag();
        int iNext = pEdge->next()->vertex()->tag();
        if (IsValueBetween(isoline.isoValue, scalarField[iNext], scalarField[iBeg]))
        {
            Enriched_Model::Point ptNext = pEdge->next()->vertex()->point();
            Enriched_Model::Point ptBeg = pEdge->opposite()->vertex()->point();
            GetInterpolationPt(ptNext, ptBeg, scalarField[iNext], scalarField[iBeg], isoline.isoValue, pt);//得到对边的插值点作为终点
            isoline.pts.push_back(pt);
            pEdge = pEdge->next()->next();

            if (IsSamePoint(pt, isoline.pts.front()))
            {
                isoline.pts.back() = isoline.pts.front();
                break;//已经成环状
            }
        }
        else if (IsValueBetween(isoline.isoValue, scalarField[iEnd], scalarField[iNext]))
        {
            Enriched_Model::Point ptNext = pEdge->next()->vertex()->point();
            Enriched_Model::Point ptEnd = pEdge->vertex()->point();
            GetInterpolationPt(ptEnd, ptNext, scalarField[iEnd], scalarField[iNext], isoline.isoValue, pt);//得到对边的插值点作为终点
            isoline.pts.push_back(pt);
            pEdge = pEdge->next();

            if (IsSamePoint(pt, isoline.pts.front()))
            {
                isoline.pts.back() = isoline.pts.front();
                break;//已经成环状
            }
        }
        else
        {
            bCircle = false;
            break;
        }

        visitFaces.push_back(pEdge->facet()->tag());
        pEdge = pEdge->opposite();
    }

    if (bCircle)
    {
        for (int i = 0; i < visitFaces.size(); ++i)
        {
            bHasIsoLine[visitFaces[i]] = true;
        }
    }

    return bCircle;
}

最后要讲的是如何遍历整个模型的三角面，使得计算出来的等值线有一定的规律？诚然，任意的顺序都能计算出等值线来。但是，会丢失一些信息。我用的广度优先搜索的顺序计算。我先找到最小的标量值，从这个顶点所在的一个面开始进行广度优先搜索遍历所有的面。同时，在队列里面保存进入队列的顺序，以进行下一步的等值线排序。最终，我得到的是一些根据标量值和进入队列顺序排序过的一堆等值线。相关代码如下：

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
67
68
69
70
//计算所有等值线
    void CSEMesh::ComputeIsoLines(std::vector<double> scalarField)
    {
        IsoLine isoline;
        std::vector<bool> bHasIsoLine(FaceNum(), false);

        m_isoLines.clear();

        //找到值最小的顶点
        Enriched_Model::Vertex_iterator pVertex = m_pModel->vertices_begin();
        double fMinValue = 1e8;
        Enriched_Model::Vertex_iterator pMinVertex;
        while (pVertex != m_pModel->vertices_end())
        {
            if (scalarField[pVertex->tag()] < fMinValue)
            {
                fMinValue = scalarField[pVertex->tag()];
                pMinVertex = pVertex;
            }
            pVertex++;
        }

        vector<bool> bVisited(FaceNum(), false);

        Enriched_Model::Facet_iterator pSourceFace = pMinVertex->vertex_begin()->facet();
        {
            struct FaceInfor
            {
                Enriched_Model::Face_iterator pFace;
                int nOrder;
                FaceInfor(Enriched_Model::Face_iterator face, int order) : pFace(face), nOrder(order){}
            };

            std::queue<FaceInfor> queFaces;
            bVisited[pSourceFace->tag()] = true;
            queFaces.push(FaceInfor(pSourceFace, 0));

            while (queFaces.empty() == false)
            {
                FaceInfor faceInfor = queFaces.front();
                queFaces.pop();

                if (ComputeIsoLine(scalarField, faceInfor.pFace, bHasIsoLine, isoline))
                {
                    isoline.order = faceInfor.nOrder;
                    m_isoLines.push_back(isoline);
                }

                Enriched_Model::Halfedge_around_facet_circulator pHalfedge = faceInfor.pFace->facet_begin();
                do
                {
                    if (pHalfedge->opposite()->is_border())
                    {
                        continue;
                    }

                    Enriched_Model::Face_iterator pFace = pHalfedge->opposite()->facet();

                    if (bVisited[pFace->tag()] == false)
                    {
                        queFaces.push(FaceInfor(pFace, faceInfor.nOrder + 1));
                        bVisited[pFace->tag()] = true;
                    }
                } while (++pHalfedge != faceInfor.pFace->facet_begin());
            }

        }

        sort(m_isoLines.begin(), m_isoLines.end());
    }

等值线效果图：

发表于 分类于 阅读次数： 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稀疏矩阵类，可以参考我其它的两篇文章。下面是第二个特征向量的分布效果，能够观察到是一个非常有规律的标量场。

发表于 分类于 ， 阅读次数： Valine：
本文字数： 1.2k 阅读时长 ≈ 1 分钟

最新版的cgal4.5，提供了个Triangulated Surface Mesh Segmentation包，利用该部分代码可以计算三维模型每个面或者顶点的sdf（Shape Diameter Function）属性。sdf的具体定义可以参考论文：L. Shapira, A. Shamir, and D. Cohen-Or. Consistent mesh partitioning and skeletonisation using the shape diameter function. The Visual Computer, 24(4):249–259, 2008.

那么sdf可以做什么了？sdf最终能表示三维模型上面的一个标量场，该标量场代表的是模型厚度的分布。因此，利用该标量场可以对模型进行分割，不管是自动化的还是交互的都可以。如何计算sdf了？可以利用定义自己去实现算法，或者利用库，比如cgal。

关于如何使用sdf标量场进行模型分割？cgal的文档说的很清楚，先用k个gmm进行软聚类（每个cluster是k-means初始化的），再利用软聚类的结果进行硬聚类，这一步实际上利用graph-cut求能量最小化。更清楚的解释，参考CGAL的相关文档。刘利刚老师用sdf进行交互式分割的文章，也是利用gmm和graph-cut作为工具实现的，原理和cgal的实现思路基本一致。

下面将一下，我使用该代码碰到的问题。计算sdf的函数是CGAL::sdf_values，这是一个复杂的模板函数，该函数的第一个参数是一个模板类Polyhedron的引用。由于cgal使用了非常复杂的模板语法，所以经常会碰到一些恶心的语法问题，编译无法通过。这一次也是。由于默认情况下，使用cgal都会继承CGAL::Polyhedron_3以实现自定义的多面体类。所以，使用这个函数的时候就会出现模板参数不匹配的情况。怎么解决了。首先，没必要去修改自己的多面体类，也不可能去修改cgal库，再重新构造个cgal默认的Polyhedron_3对象也很浪费。

其实，这里可以考虑到转型。类型转换有两种考虑，转值和转指针。如果将自定义多面体类转换为默认的多面体类，肯定会生成个具大的临时对象，多浪费啊。所以了，可以先取地址转换为基类指针再解引用就获得了cgal默认的Polyhedron_3对象了。我的代码如下：std::pair min_max_sdf = CGAL::sdf_values((Polyhedron)m_pModel, m_sdf_property_map);

至于其它部分的使用，参考cgal手册即可，非常简单方便。下面再贴几张图，展示下sdf标量场的分布，以及利用sdf的分割结果。

sdf标量场，白色的是等值线。

sdf分割模型结果：

至于分割边界不是很整齐，这只是跟模型密度有关系，模型越密，能够得到越平滑的边界，也有相关的算法优化边界。

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

上一篇文章讲了如何用MATLAB计算稀疏矩阵的特征向量，但是我们最终的目的是使用C++达到这个需求。大致搜索了下，发现使用C++调用MATLAB计算引擎的方式最方便，就采用了，目前的计算速度和方便性都能满足要求。下面讲一下如何操作。

1.首先是安装MATLAB，并且保证版本正确，比如32位程序只能调用32位的MATLAB，如我的机器是64位win8.1，所以必须强制安装32位的MATLAB，才能配合我用vs2010编写的32位程序。如果，我用vs2010开放64位程序，那么当然是可以使用64MATLAB的，但是配置64位开源库是个大工程，有些库并没有针对64位版本测试过。

2.假设我的MATLAB安装在C:\Program Files (x86)\MATLAB下面，那么include目录是C:\Program Files (x86)\MATLAB\R2013a\extern\include，lib目录是C:\Program Files (x86)\MATLAB\R2013a\extern\lib\win32\microsoft，bin目录是C:\Program Files (x86)\MATLAB\R2013a\bin\win32。include和lib在vs里面设置好就行了，至于bin目录需要加入到path环境变量中，最后重启电脑生效，当然你一个个去找需要的dll也是可以的。

3.配置好之后，就可以使用MATLAB引擎了。大致需要1>打开引擎，2>用C++创建变量，3>把变量对应到MATLAB命令中，4>执行MATLAB命令，5>将MATLAB变量传回C++这样的一些操作。比如，engOpen打开引擎，mxCreateSparse用于创建稀疏矩阵，engPutVariable将变量传入MATLAB，engEvalString执行命令，engGetVariable返回结果。具体怎么操作还是得参考相应教程。另外，可以用C++创建不同类型的变量，稀疏矩阵是比较特殊的一种。

4.至于如何用MATLAB计算稀疏矩阵的特征向量，参阅上一篇文章。

5.关键的一步，如何设置mxCreateSparse返回的稀疏矩阵的值。先看下面的代码吧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
mxArray* sm = mxCreateSparse(m_nRow, m_nColumn, m_nNum, mxREAL);
double* pr = mxGetPr(sm);
mwIndex* ir = mxGetIr(sm);
mwIndex* jc = mxGetJc(sm); 

int k = 0;
int nLastC = -1;
for (int i = 0; i < m_nNum; ++i)
{
    pr[i] = m_pSI[i].v;
    ir[i] = m_pSI[i].r;

    if (m_pSI[i].c != nLastC)
    {
        jc[k] = i;
        nLastC = m_pSI[i].c;
        ++k;
    }
}
jc[k] = m_nNum;

pr对应的是值，ir对应的是行，关键的是jc，jc并不是列，而是到当前列为止出现的值数目，具体可以查阅MATLAB文档，或者理解下以上代码，正确初始化稀疏矩阵是成功的关键。
下面随附我的计算拉普拉斯矩阵FiedlerVector的函数。

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
67
void CSparseMatrix::GetFiedlerVector(double* pVector, int nChoose)
    {
        assert(nChoose >= 0  nChoose <= 5);

        Engine *ep = NULL;
        if (!(ep = engOpen(NULL)))//打开matlab引擎
        {
            fprintf(stderr, "\nCan't start MATLAB engine\n");
            return;
        }

        mxArray* sm = mxCreateSparse(m_nRow, m_nColumn, m_nNum, mxREAL);
        double* pr = mxGetPr(sm);
        mwIndex* ir = mxGetIr(sm);
        mwIndex* jc = mxGetJc(sm); 

        int k = 0;
        int nLastC = -1;
        for (int i = 0; i < m_nNum; ++i)
        {
            pr[i] = m_pSI[i].v;
            ir[i] = m_pSI[i].r;

            if (m_pSI[i].c != nLastC)
            {
                jc[k] = i;
                nLastC = m_pSI[i].c;
                ++k;
            }
        }
        jc[k] = m_nNum;

        engPutVariable(ep, "sm", sm);
        engEvalString(ep, "max_d = eigs(sm, 1);");//获得最大绝对值的特征值

        double fMaxE;
        mxArray* max_d = engGetVariable(ep, "max_d");
        memcpy(fMaxE, (void *)mxGetPr(max_d), sizeof(double));//获得最大的特征值(假设特征值是从最大排列到0)

        //更新对角线元素，相当于sm = sm - fMaxE，这样会使计算出来的特征值都减去fMaxE
        pr = mxGetPr(sm);
        for (int i = 0; i < m_nNum; ++i)
        {
            pr[i] = m_pSI[i].v;
            if (m_pSI[i].r == m_pSI[i].c)//对角线元素
            {
                pr[i] -= fMaxE;
            }
        }

        engPutVariable(ep, "sm", sm);
        engEvalString(ep, "[V,D] = eigs(sm);");
        engEvalString(ep, "[newD,IX] = sort(diag(D));");
        engEvalString(ep, "newV=V(:,IX);");

        const int EVAl_STR_LEN = 100;
        char szEvalString[EVAl_STR_LEN];
        sprintf(szEvalString, "Fiedler = newV(:, %d);", nChoose);
        engEvalString(ep, szEvalString);

        mxArray* fiedler = engGetVariable(ep, "Fiedler"); //返回计算结果
        memcpy(pVector, (void *)mxGetPr(fiedler), sizeof(double) * m_nRow);

        mxDestroyArray(sm);
        mxDestroyArray(fiedler);
        engClose(ep);
    }

参照以上函数基本可以知道本文所讲的内容了，一切尽在代码中。

发表于 分类于 阅读次数： Valine：
本文字数： 825 阅读时长 ≈ 1 分钟

首先，说明下为什么要用MATLAB。原因是，支持稀疏矩阵的矩阵库比较少，即使支持去学习使用时间也更多。所以，就用MATLAB了。目的是使用MATLAB计算稀疏矩阵的特征向量，再用C++调用MATLAB代码。

问题1：如何定义稀疏矩阵。MATLAB使用sparse函数创建稀疏矩阵。如下调用：S = sparse(i,j,s,m,n)，根据参数创建稀疏矩阵。i和j代表不为0元素的行和列值，s则是对应的值，m和n是矩阵的大小。当然还有其它接口，具体查阅MATLAB文档。

问题2：怎么计算稀疏矩阵的特征向量？MATLAB提供了eigs函数。其实，还有个eig函数，经本人实验，在使用C++调用MATLAB引擎的过程中，该函数失败了。那么，就只能使用eigs函数了。eigs函数有几种形式。下面说下我用到的几种形式。

形式1：d = eigs(A)，返回的是稀疏矩阵A的特征值（所有的特征值构成一个向量d，本人的理解是没有保证d是有序的）。

形式2：d = eigs(A,k)，返回的是最大的k个特征值，文档也没有说是有序的，如果设置k为1，那么返回的是最大的特征值。

形式3：[V,D] = eigs(A)，V是特征向量构成的矩阵，每列是一个特征矩阵，D则是一个对角矩阵，对角线上的值是和V对应的特征值。不过，V和D的大小都固定为6。

形式4：[V,D] = eigs(A，k)，则可以返回最大的k个特征向量和对应的特征值，但是也没有排序，需要根据D排序V，以方便使用。

问题3：如何实现返回k个根据特征值排序后的特征向量？这是最有价值的问题了。直接提供代码吧。

1
2
3
[V,D] = eigs(sm,k);
[newD,IX] = sort(diag(D)));;
newV=V(:,IX);

newV就是根据特征值排序后的特征向量，sm当然需要先建立了。
写在最后的话，之所以使用MATLAB计算稀疏矩阵的特征向量，是因为可以用C++调用。所以，关于这篇文章这些知识就够了，也算是我为自己做的整理吧。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 4.9k 阅读时长 ≈ 4 分钟

前面提到priority_queue的内部使用了stl的泛型heap算法。现在让我们看看怎么来实现heap。至于heap，就不作仔细介绍了，本质上是用一个数组表示的完全二叉树，并且父节点总是大于（或者小于）子节点的值。我们先来看push_heap。

下图是push_heap的示意图。

从图中可以看到算法的过程是将新加入堆的值（50），层层上挪，直到正确的位置。下面来看，摘录出来的代码。

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
template <class RandomAccessIterator>
inline void push_heap(RandomAccessIterator first, RandomAccessIterator last) {
    // 注意，调用该函数时候，新元素位于最后一个位置(last-1)。
    __push_heap_aux(first, last, distance_type(first), value_type(first));
}

template <class RandomAccessIterator, class Distance, class T>
inline void __push_heap_aux(RandomAccessIterator first,
    RandomAccessIterator last, Distance*, T*) {
    __push_heap(first, Distance((last - first) - 1), Distance(0),
        T(*(last - 1)));
    // (last-first)–1代表新元素的索引，0是堆首的索引，*(last - 1)是新加入的值
}

template <class RandomAccessIterator, class Distance, class T>
void __push_heap(RandomAccessIterator first, Distance holeIndex,
    Distance topIndex, T value) {
    Distance parent = (holeIndex - 1) / 2;  // 找出父節點
    while (holeIndex > topIndex  *(first + parent) < value) {
        // 尚未到达顶端，且父节点小于新值
        // 由于以上使用 operator<，可知 STL heap 是max-heap
        *(first + holeIndex) = *(first + parent);   // 令洞值为父值
        holeIndex = parent; // percolate up：调整洞号，向上提升至父节点。
        parent = (holeIndex - 1) / 2;   // 新洞的父节点
    }    // 持续至顶端，或满足 heap 的次序特性为止。
    *(first + holeIndex) = value;   // 令洞值为新值。
}

push_heap的用法是输入迭代器对，并且保证[first,last-1)是最大堆，*(last-1)是新加入的元素。push_heap调用辅助函数push_heap_aux。至于为什么需要这个辅助函数了？应该是为了提取出distance_type和value_type吧，这两个内联函数的定义，可以参考stl源码剖析迭代器的那章。下面来思考真正的实现函数push_heap。这个函数需要新加入元素位置holeIndex和堆首位置topIndex，另外还有保存好的新加入值。算法的过程很简单，就是上溯holeIndex，找到真正满足条件的位置（无法继续上回溯），然后把value放入该位置即可。
pop_heap实际上是一个相反的过程。实现思路是将堆大小加一后，再找出最后一个元素应该放入的位置holeIndex，最后再加入这个值。示意图如下：


下面看看摘录出来的代码，思想类似于push_heap，只需要求出最终的holeIndex。

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
template <class RandomAccessIterator>
inline void pop_heap(RandomAccessIterator first, RandomAccessIterator last) {
    __pop_heap_aux(first, last, value_type(first));
}

template <class RandomAccessIterator, class T>
inline void __pop_heap_aux(RandomAccessIterator first,
    RandomAccessIterator last, T*) {
    __pop_heap(first, last - 1, last - 1, T(*(last - 1)), distance_type(first));
    // pop动作的結果为底层容器的第一個元素。因此，首先设定欲调整值为尾值，然后將首值調至 
    // 尾节点（所以以上將迭代器result设为last-1）。然后重整 [first, last-1)，
    // 使之重新成一個合格的 heap。
}

template <class RandomAccessIterator, class T, class Distance>
inline void __pop_heap(RandomAccessIterator first, RandomAccessIterator last,
    RandomAccessIterator result, T value, Distance*) {
    *result = *first; // 設定尾值为首值，于是尾值即是結果，
    // 可由调用底层容器之 pop_back() 取出尾值。
    __adjust_heap(first, Distance(0), Distance(last - first), value);
    // 以上欲重新調整 heap，洞号为 0，欲調整值为value。
}

template <class RandomAccessIterator, class Distance, class T>
void __adjust_heap(RandomAccessIterator first, Distance holeIndex,
    Distance len, T value) {
    Distance topIndex = holeIndex;
    Distance secondChild = 2 * holeIndex + 2;   // 洞节点之右子节点
    while (secondChild < len) {
        // 比较洞节点之左右兩个子值，然后以 secondChild 代表较大子节点。
        if (*(first + secondChild) < *(first + (secondChild - 1)))
            secondChild--;
        // Percolate down：令较大大子值为洞值，再令洞号下移至较大子节点处。
        *(first + holeIndex) = *(first + secondChild);
        holeIndex = secondChild;
        // 找出新洞节点的右子节点
        secondChild = 2 * (secondChild + 1);
    }

    if (secondChild == len) { // 沒有右子节点，只有左子节点
        // Percolate down：令左子值为洞值，再令洞号下移至左子节点处。
        *(first + holeIndex) = *(first + (secondChild - 1));
        holeIndex = secondChild - 1;
    }

    // 將欲调整值填入目前的洞号內。注意，此時肯定滿足次序特性。
    // 依侯捷之见，下面直接改為 *(first + holeIndex) = value; 应该可以。
    __push_heap(first, holeIndex, topIndex, value);
}

类似于push_heap，pop_heap也是调用辅助函数pop_heap_aux。pop_heap_aux调用pop_heap。pop_heap调用adjust_heap调整holeIndex，最终在holeIndex处放入value（原最后一个的值）。关键代码是adjust_heap中的循环。循环的主要意思是将holeIndex不断下放，直到最底层。最后的if语句的意思是，如果最底层有左子节点，而没有右子节点，那么最终位置肯定是这个左子节点。最后一句代码的意思是加入value到holeIndex，由于已经调整完毕，所以一个赋值操作也可以达到要求，参见侯捷注释。
sort_heap就比较简单了，不断将极值移动到末尾，不断pop_heap。

1
2
3
4
5
6
7
8
9
// 以下這個 sort_heap() 不允許指定「大小比較標準」
template <class RandomAccessIterator>
void sort_heap(RandomAccessIterator first, RandomAccessIterator last) {
    // 以下，每執行一次 pop_heap()，極值（在STL heap中為極大值）即被放在尾端。
    // 扣除尾端再執行一次 pop_heap()，次極值又被放在新尾端。一直下去，最後即得
    // 排序結果。
    while (last - first > 1)
        pop_heap(first, last--); // 每執行 pop_heap() 一次，操作範圍即退縮一格。
}

最后要将的是make_heap，即将一个迭代器对里面的内容构造为最大堆。代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
// 將 [first,last) 排列為一個 heap。
template <class RandomAccessIterator>
inline void make_heap(RandomAccessIterator first, RandomAccessIterator last) {
    __make_heap(first, last, value_type(first), distance_type(first));
}

// 以下這組 make_heap() 不允許指定「大小比較標準」。
template <class RandomAccessIterator, class T, class Distance>
void __make_heap(RandomAccessIterator first, RandomAccessIterator last, T*,
    Distance*) {
    if (last - first < 2) return;   // 如果長度為 0 或 1，不必重新排列。
    Distance len = last - first;
    // 找出第一個需要重排的子樹頭部，以 parent 標示出。由於任何葉節點都不需執行 
    // perlocate down，所以有以下計算。parent 命名不佳，名為 holeIndex 更好。
    Distance parent = (len - 2) / 2;

    while (true) {
        // 重排以 parent 為首的子樹。len 是為了讓 __adjust_heap() 判斷操作範圍
        __adjust_heap(first, parent, len, T(*(first + parent)));
        if (parent == 0) return;    // 走完根節點，就結束。
        parent--;                   // （即將重排之子樹的）頭部向前一個節點
    }
}

make_heap中代码的思路也很简单。从原序列的中间位置开始不断调整（调用adjust_heap），每次调整的目的是以当前位置为根的构建一个子堆。至于为什么从中间位置开始就可以了？原因很简单，最底层元素的数目大致就会占了一半了。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 1.3k 阅读时长 ≈ 1 分钟

类似于queue，stack也是个简单的适配器。在实现上，stack与queue非常类似。底层可以使用同样的容器，比如默认都采用deque。因此，stack的源码和queue的非常类似。代码过于简单，都可以自解释了。下面是摘录出来的代码，可以顺便编译通过。

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
#ifndef _STL_STACK_H
#define _STL_STACK_H

namespace std
{
    template <class T, class Sequence = deque<T> >
    class stack {
        friend bool operator==(const stack, const stack);
        friend bool operator<(const stack, const stack);

    public:
        typedef typename Sequence::value_type value_type;
        typedef typename Sequence::size_type size_type;
        typedef typename Sequence::reference reference;
        typedef typename Sequence::const_reference const_reference;
    protected:
        Sequence c; // 底层容器

    public:
        // 以下完全利用 Sequence c 的操作，完成 stack 的操作。
        bool empty() const { return c.empty(); }
        size_type size() const { return c.size(); }
        reference top() { return c.back(); }
        const_reference top() const { return c.back(); }

        // deque 是两头可进出，stack 是末端进，末端出（所以后进者先出）。
        void push(const value_type x) { c.push_back(x); }
        void pop() { c.pop_back(); }
    };

    template <class T, class Sequence>
    bool operator==(const stack<T, Sequence> x, const stack<T, Sequence> y) {
        return x.c == y.c;
    }

    template <class T, class Sequence>
    bool operator<(const stack<T, Sequence> x, const stack<T, Sequence> y) {
        return x.c < y.c;
    }

};

#endif

代码非常类似于queue，除了pop弹出的是末尾元素，top也是返回末尾元素之外。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 1.5k 阅读时长 ≈ 1 分钟

接着上面一篇，这次要讲的是另一个适配器容器queue。queue的操作很简单，先进先出。现在能够使用底层容器来实现，所以也非常简单。总之，这就是个简单的适配器吧。queue需要两个模版参数，类型和底层容器，默认的底层容器是deque。摘录出的代码如下：

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
#ifndef _STL_QUEUE_H
#define _STL_QUEUE_H
#include <deque>

namespace std
{
    template <class T, class Sequence = deque<T> >
    class queue {
        friend bool operator== (const queue x, const queue y);
        friend bool operator< (const queue x, const queue y);

    public:
        typedef typename Sequence::value_type value_type;
        typedef typename Sequence::size_type size_type;
        typedef typename Sequence::reference reference;
        typedef typename Sequence::const_reference const_reference;

    protected:
        Sequence c; // 底层容器

    public:
        // 以下完全利用 Sequence c 的操作，完成 queue 的操作。
        bool empty() const { return c.empty(); }
        size_type size() const { return c.size(); }
        reference front() { return c.front(); }
        const_reference front() const { return c.front(); }
        reference back() { return c.back(); }
        const_reference back() const { return c.back(); }

        // deque 是两头可进出，queue 是末端进，前端出（所以先进者先出）。
        void push(const value_type x) { c.push_back(x); }
        void pop() { c.pop_front(); }
    };

    template <class T, class Sequence>
    bool operator==(const queue<T, Sequence> x, const queue<T, Sequence> y) {
        return x.c == y.c;
    }

    template <class T, class Sequence>
    bool operator<(const queue<T, Sequence> x, const queue<T, Sequence> y) {
        return x.c < y.c;
    }

};

#endif

总之，代码非常简洁。分析下代码，typedef typename Sequence::value_type value_type;这其实，就是声明嵌套模版类型。还有必须注意的是，声明嵌套模版类型，必须加typename关键字，否则会将其视作定义。

发表于 分类于 ， ， 阅读次数： Valine：
本文字数： 2.2k 阅读时长 ≈ 2 分钟

暑假时候好好阅读过，stl源码剖析这本书，之后又看了c++标准程序库，effective c++，more
effecttive c++。虽然每次都如醍醐灌顶般，明白了很多知识概念，有种拨开迷雾的感觉，但是终究抵不过
遗忘，所以了还是随便写写吧。
首先从实现代码差不多最短的优先队列开始，虽然简单的代码后面依赖很庞大的东西，比如泛型
迭代器，还有泛型二叉堆算法。我是这样理解优先队列的。首先，数据是存在一个泛型序列容器(线性容器)，
这个可以指定，默认是vector，当然指定list和deque也是可以的。另外还需要给数据类型提供个小于比较仿
函数，有些类型有默认的less泛型仿函数。所以了，对模版比较熟悉的同学，很自然就猜到了，priority_queue
的外观->一个模版（类型，线性容器，less仿函数）。
优先队列支持哪些操作了？什么empty，size，top，直接转交给底层容器就行了，所谓的适配器模式。
只需要思考，如何实现push和pop操作。对算法熟悉的同学，都知道底层数据的组织是用二叉堆的，而且是最大
堆。问题是，我们还需不需要实现个堆操作了。不需要了，stl做了所有的事情。push_heap就是将输入范围的最
后一个元素作为在原堆后面加入的新元素，再重构堆。pop_heap则是将堆首元素交换到堆末，再重构缩小之后的堆。
至于二叉堆的原理，就不仔细介绍了。其实也比较简单，如果不去思考一些复杂度的证明。
到最后，priority_queue的代码可以很清爽的出来了。尤其注意下面的那些typedef，可以说是模版
trick的常见手法，习惯了就好。这种写法还有个用法是为了使内嵌模版声明可见，记得effect c++提到过。以下
源码基本来自stl源码剖析随书奉送代码，做了些修改。

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
#include <algorithm>

template <typename T, typename Sequence = vector<T>,
    typename Compare = less<typename Sequence::value_type >>
class  priority_queue
{
public:
    typedef typename Sequence::value_type value_type;
    typedef typename Sequence::size_type size_type;
    typedef typename Sequence::reference reference;
    typedef typename Sequence::const_reference const_reference; //为了使内嵌模版声明可见

protected:
    Sequence c; // 底层容器
    Compare comp;   // 元素大小比较仿函数

    priority_queue() : c() {}
    explicit priority_queue(const Compare x) : c(), comp(x) {}

    template <class InputIterator>
    priority_queue(InputIterator first, InputIterator last, const Compare x)
        : c(first, last), comp(x) {
        make_heap(c.begin(), c.end(), comp);
    }

    template <class InputIterator>
    priority_queue(InputIterator first, InputIterator last)
        : c(first, last) {
        make_heap(c.begin(), c.end(), comp);
    }

    bool empty() const { return c.empty(); }
    size_type size() const { return c.size(); }
    const_reference top() const { return c.front(); }

    void push(const value_type x)
    {
        try
        {
            c.push_back(x);
            push_heap(c.begin(), c.end(), comp);
        }
        catch (...)
        {
            c.clear();
            throw;
        }
    }

    void pop()
    {
        try
        {
            pop_heap(c.begin(), c.end(), comp);
            c.pop_back();
        }
        catch (...)
        {
            c.clear();
            throw;
        }
    }
};

再啰嗦几句，刚看源码剖析的时候，对某些模版写法很头疼，但是看完书之后，另外看了几本其它的书，
就对这些写法习以为常了，觉得不这么写才怪了。有点艺术的感觉。

发表于 分类于 阅读次数： Valine：
本文字数： 454 阅读时长 ≈ 1 分钟

这里统一考虑二维直线和三维平面的情况。
假设法线（二维直线和三维平面的法线）为n，其上的任意一点为p0。那么，
可以得到方程为：(p - p0) * n = 0。展开得到p * n = p0 * n。令p0 * n = d，由此可以得到
方程可以简单的表示为p * n = d，即(x, y, z) * (nx ,ny, nz) = d，二维情况直线的为(x, y) * (nx, ny) = d。
至此，我们可以方便的用(n,d)代表任意的二维直线或者三维平面了。
下面，开始应用这个表示。
比如，求点到直线或者平面的距离，然后求投影点，再求对称点。
最关键的是如何方便求出距离。假设我们要求q到直线的距离。我们现在的直线方程为
p * n = d。那么经过q与当前直线的平行直线的法线也为n。由于q在经过q的平行直线上面，有q * n = d’。
有符号距离为 d - d’ = d - q * n。通过这个式子，也可以知道d其实就是原点到直线或者平面的有符号距离。
那么投影点q’ = q + (d - d’) = q + (d - q * n)。
对称点则为 q’’ = q’ + (d - d’) = q + 2 * (d - q * n)。
利用这种表示方法，基本上可以直接给出答案。