STEP BY STEP

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

这是一个很老的话题了。所谓的反射机制，其实就是MFC里面的对象动态创建。比如说，给一个字符串代表类名，比如”CMesh”，根据这个字符串创建出一个CMesh对象。有不同的实现方法，在MFC里面是用一张巨大的链表网络把所有的类的信息链接起来，里面存储了类名和构造函数，动态创建的时候按照一定的顺序去里面搜索创建函数（构造函数）。

与其用一张链表网，不如用一个字典保存起来，键是类名字符串，值是对应的构造函数。看起来很简单的样子，确实也很简单额。只是实现起来有些trick，比如字典一般是类的静态成员等。。。假设只有无参数的构造函数，首先定义上面所述的内容，类的字典。

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
//头文件
    typedef void* (*CreateFuntion)(void);

    class CClassFactory
    {
    public:
        CClassFactory() {}
        ~CClassFactory() {}

    public:
        static void* GetClassByName(std::string name)
        {
            std::map<std::string, CreateFuntion>::const_iterator find;
            find = m_clsMap.find(name);  

            if(find==m_clsMap.end())
            {  

                return NULL;
            }
            else
            {
                return find->second();
            }  

        }  

        static void RegistClass(std::string name, CreateFuntion method)
        {
            m_clsMap.insert(std::make_pair(name,method));
        }  

    private:
        static std::map<std::string,CreateFuntion> m_clsMap;
    };
//源文件
std::map<std::string, CreateFuntion> CClassFactory::m_clsMap;

到目前为止，使用CClassFactory即可完成工作了。代码实在过于简单，不用解释了。下面看看一些带trick的东西。

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
    class RegistyClass
    {  

    public:
        RegistyClass(std::string name, CreateFuntion method)
        {
            CClassFactory::RegistClass(name, method);
        }
    };  

    template<class T, const char* name>
    class Register
    {
    public:
        Register()
        {
            const RegistyClass tmp = rc;//这一句不能删除，如果不使用rc，那么rc这个模板成员编译器不会生成代码
        }  

        static void* CreateInstance()
        {
            return new T;
        }  

    public:
        static const RegistyClass rc;
    };  

    template<class T, const char* name>
    const RegistyClass Register<T, name>::rc(name, Register<T, name>::CreateInstance);

#define REGISTRY_CLASS(class_name) \
extern char const array_##class_name[] = #class_name;\
Register<class_name, array_##class_name> register_class_name

RegistyClass不用说了，模板类Register是关键。模板T代表要注册的类类型，常量字符串模板name代表类名。其中的static成员rc用于注册当前类。注意CreateInstance的实现，说明只能调用无参数构造函数。至于Register构造函数里面的那句话也是不能省略的。

最麻烦的是下面的那个REGISTRYCLASS宏，其余代码都来自网络，这个宏确是我修改实现的。class_name代表类，比如CMesh，注意不是类名。这里有几个关键点，一个是extern。在vs2010中，不加extern修饰无法编译通过，具体原因估计是编译器支持不够吧。另一个是#class_name，是把CMesh转换为”CMesh”。最最关键的是`array##class_name`，是生成array_CMesh的变量名。为什么要这么做了，这样注册不同的类就是不同的变量名，避免了变量重定义，否则会链接失败。使用方法是在要注册的类源文件中，包含ClassFactory.h头文件，再添加宏REGISTRY_CLASS(CMesh);即可。

网上提供的代码，还定义了声明类的宏，表示继承Register模板类自动支持反射机制。相关代码经过修改如下所示，这个就没有去具体使用了，因为修改类定义比较麻烦，添加一句宏到现在代码中却是很方便的事情。这个反射机制已经用到我现在的项目中了，实现了一个从配置文件读取信息，选择创建不同的界面的功能。

1
2
3
4
5
6
7
#define DEFINE_CLASS(class_name)\
extern char char array_##class_name[]=#class_name;\
class class_name:public Register<class_name, NameArray>  

#define DEFINE_CLASS_EX(class_name, father_class) \
extern char char array_##class_name[] = #class_name;\
class class_name: public Register<class_name, NameArray>, public father_class

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

我终于把界面换成qt的了，即使我从大二开始就用mfc，一直用到去年上半年。也终于到了，我实在受不了mfc代码混乱的时候了，qt那么方便的东西，我为什么不早点用了。这里就不吐槽mfc的混乱，qt的方便了。在mfc里面建立一个复杂的界面，比如说有dock窗口，或者tab页面的，代码多，而且混乱麻烦，更无语的是，每次我都必须去搜索下，解决各种各样的问题，有的还可以理解，有的就是莫名其妙了。更别说什么代码清晰，跨平台之类的，何止是天方夜谭。

下面，说说我自己逐渐搭建的这个框架吧。其实，这个框架是从我上一个项目里，使用mfc的单文档单视图，主窗口里面有几个dockpane，视图里面放置了渲染OpenGL的子窗口，过渡来的。在qt的框架里面，mfc所有的这些恶心的东西都没有了，只剩下一个MainWindow类。在这个类里面，创建菜单，工具条，状态条，主窗口，dock窗口，tab窗口。代码简洁方便。渲染opengl的窗口设置为MainWindow的主widget。渲染窗口只是一个普通的widget而已，可以任意放置。

这样的做法能够把界面和渲染窗口独立出来，渲染窗口对应一个独立的类(COglWidget)，所以方便创建不同类型的复杂界面。我觉得我应该先用viso画个类图，才能说清楚。下图是我的界面。中间是COglWidget。右侧是一个DockWidget，里面放置了个设置Mesh属性的widget。

通过主窗口的构造函数，可以看看如何创建这个界面，虽然在主窗口里面做这么的事情不是个好的习惯。

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
CMainWindow::CMainWindow(QWidget *parent)
        : QMainWindow(parent)
    {
        //setWindowFlags(this->windowFlags()  ~Qt::WindowMaximizeButtonHint);//隐藏最大化按钮
        setWindowTitle(tr("Skeleton Extract"));

        m_pSEMesh = new CSEMesh();
        m_pOglForSE = new COglWidget(this, m_pSEMesh);
        m_pSEMesh->SetOglWidget(m_pOglForSE);
        setCentralWidget((QWidget*)m_pOglForSE);//设置中间是渲染OpenGL的widget      

        m_pPCMesh = new CPointCloudMesh();
        m_pOglForPC = new COglWidget(this, m_pPCMesh);
        //m_pOglForPC->SetEye(0.0, 0.0, 1.0);
        m_pOglForPC->setFixedSize(960, 720);
        m_pOglForPC->setWindowTitle(tr("Point Cloud For Geodesic Matrix Eigen Vector"));
        m_pOglForPC->move((QApplication::desktop()->width() - m_pOglForPC->width()) / 2,
            (QApplication::desktop()->height() - m_pOglForPC->height())/2);
        m_pOglForPC->hide();
        m_pOglForPC->setWindowFlags(Qt::Window);

        CreateActions();
        CreateMenus();
        CreateToolBars();
        CreateStatusBar();
        CreateDockWidgets();

        setMinimumSize(960, 720);
        setWindowState(Qt::WindowMaximized);

        std::string strName = "MeshData\\bunny_3k.m";
        m_pSEMesh->LoadModel(strName.c_str());
        m_pMeshNameLabel->setText(QString(strName.c_str()));
    }

再来看一个使用tabwidget的界面，记得在mfc中使用tab页面应该是比较麻烦，我也没用过。不过，在qt里面，唯一需要改变的是把tabwidget设置为central的widget，把其它widget加入到tab页面里面就行了。效果如图：

这个界面的创建代码如下，这里我把菜单，工具条，状态栏都去掉了。

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
CMainWindow::CMainWindow(QWidget *parent)
        : QMainWindow(parent)
    {
        //setWindowFlags(this->windowFlags()  ~Qt::WindowMaximizeButtonHint);//隐藏最大化按钮
        setWindowTitle(tr("Virtual Box"));

        m_pMesh = new CCubeMesh();
        m_pOglWidget = new COglWidget(this, m_pMesh);

        m_pMainViewMesh = new CCubeMesh();
        m_pMainViewMesh->DrawStyle(CCubeMesh::MAIN_VIEW);
        m_pMainViewWidget = new CViewWidget(this, m_pMainViewMesh);

        m_pTopViewMesh = new CCubeMesh();
        m_pTopViewMesh->DrawStyle(CCubeMesh::TOP_VIEW);
        m_pTopViewWidget = new CViewWidget(this, m_pTopViewMesh);

        m_pSideViewMesh = new CCubeMesh();
        m_pSideViewMesh->DrawStyle(CCubeMesh::SIDE_VIEW);
        m_pSideViewWidget = new CViewWidget(this, m_pSideViewMesh);

        //CreateActions();
        //CreateMenus();
        //CreateToolBars();
        //CreateStatusBar();
        //CreateDockWidgets();

        setMinimumSize(960, 720);
        setWindowState(Qt::WindowMaximized);

        //std::string strName = "cube.obj";
        //m_pMesh->LoadModel(strName.c_str());
        //m_pMeshNameLabel->setText(QString(strName.c_str()));
        m_pTabWidget = new QTabWidget(this);
        m_pTabWidget->addTab(m_pOglWidget, "Cube View");
        m_pTabWidget->addTab(m_pMainViewWidget, "Main View");
        m_pTabWidget->addTab(m_pTopViewWidget, "Top View");
        m_pTabWidget->addTab(m_pSideViewWidget, "Side View");
        setCentralWidget((QWidget*)m_pTabWidget);
    }

从上面的代码中，可以看到创建一个OglWidget需要一个Mesh指针。在我的框架中，COglWidget内中只使用几个固定的Mesh函数。因此，可以把COglWidget实现为依赖于一个IMesh接口，这样就进一步把OpenGL窗口和Mesh分离开来。不管Mesh类的具体实现如何，OpenGL窗口只需要依赖IMesh接口，因此，第一个界面里面的Mesh类实际上组合了一个CGAL的CGAL::Polyhedron_3对象。而第二个界面的Mesh类是一个简单的渲染长方体的类，不过支持渲染多个视图而已。
下面说说，COglWidget类具有哪些功能。首先，渲染Mesh，包括设置视口，eye，投影，光照，背景等等，其次，处理鼠标按键操作模型，我的实现是把所有的鼠标按键消息都丢给了下层的一个trackball类，这个trackball类支持旋转缩放移动模型，另外还写了几个打开文件，以及选取背景色的操作。
现在到了Mesh的部分了。因为跟OpenGL窗口相关的只是IMesh接口，所以，我可以任意实现实际的Mesh类。
IMesh的定义如下，

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
#ifndef IMESH_H
#define IMESH_H
#include "libyx/vertex.h"

namespace LibYX
{
    struct IMesh
    {
        IMesh()
        {
            m_pMatMV = new double[16];
            m_pMatProj = new double[16];
            m_pViewport = new int[4];
        }

        virtual ~IMesh()
        {
            FreeModel();
            delete[] m_pMatMV;
            delete[] m_pMatProj;
            delete[] m_pViewport;
        }

        virtual bool LoadModel(const char* pFilename) { return true; }
        virtual void SaveModel(const char* sFilename) {}
        virtual void UnifyModel() {}

        virtual bool HasModelLoad() const { return true; }
        virtual void FreeModel() {}

        virtual void InitData() {}
        virtual void DestroyData() {}

        virtual void DrawScene()
        {
            glGetDoublev(GL_MODELVIEW_MATRIX, m_pMatMV);
            glGetDoublev(GL_PROJECTION_MATRIX, m_pMatProj);
            glGetIntegerv(GL_VIEWPORT, m_pViewport);
        }

        void WindowWidth(int nWidth) { m_nWinWidth = nWidth; }
        int WindowWidth() const { return m_nWinWidth; }
        void WindowHeight(int nHeight) { m_nWinHeight = nHeight; }
        int WindowHeight() const { return m_nWinHeight; }

    protected:
        double* m_pMatMV;//模型视图矩阵
        double* m_pMatProj;//投影矩阵
        int* m_pViewport;//视口
        int m_nWinWidth;
        int m_nWinHeight;

    protected:
        v3d m_vMax, m_vMin;         // bounding box of the scene
    };
};

#endif

其中，DrawScene()是渲染接口，另外还有加载模型的接口等，里面还有些数据。我用CTriMesh继承IMesh接口，同时组合Enriched_Model *m_pModel的指针，创建一个CGAL的多边形网格模型，实际的操作基本都交给m_pModel。这个类的定义太复杂了，就不贴了。这个类里面我集成了基本操作，比如按照点，线，面方式渲染，渲染包围盒，坐标轴，法线，顶点，边，面的索引等，加载和保存模型等等。因此，CTriMesh作为模型的基类，在实际的项目中，再继承CTriMesh，比如第一个界面中的CSEMesh，在CSEMesh中进行具体的处理。
所以，实际的算法实现都在CSEMesh里面。因为，IMesh和CTriMesh，COglWidget可以固化了。能够变化的只是CMainWindow。如果要添加新的代码，只能继承CTriMesh，在子类里面实现具体的算法，也算是对修改关闭对扩展开放吧。

因此，如果再做类似的东西，我就有一个很方便的框架了。。。鉴于不少人留言要代码，有兴趣的可以参考这个QtRenderMesh。

发表于 分类于 阅读次数： 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关键字，否则会将其视作定义。