STEP BY STEP

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

很久没骑这么远了。第一次骑三天都不知道是怎么坚持过来的。总的感觉，身体比去年好太多了。虽然脚累，但是呼吸完全不乱，完全不会心累。中午吃错东西了，一开始就来个闹肚子、过敏。我对某种食物强过敏，当时就在为自己生命安全思考了。本打算骑到附近的梅溪湖算了，但是吹了下风，觉得过敏消散了，就坚持走走看吧。

雷锋大道真的好远，从河西走完整个望城。经过去年的洗心禅寺，发现好近额。为啥去年我走这么点距离会觉得好累了。有些改变确实在不知不觉。不管是身，还是心。永远健康充满活力的身体，是追寻幸福的必要条件。

最后，附图一张。

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

排序的重要性就不用说了。stl里面的很多容器都是自序的，比如set和map，从begin到end就是一个顺序，毕竟最基本的概念就是一颗搜索树。另外，stack、queue、priority_queue的顺序也是内部实现决定的，不用多说了。list则不支持随机迭代器，内部也提供sort函数。所以，stl里面的sort泛型算法针对的是vector，deque，原生数组之类的数据。正常情况下，用的最多的也就是vector了。

插入排序和快速排序大家都知道是什么。但是，stl里面的sort并不仅仅是一个快速排序，而是一个复杂的组合体，有人把它叫做intro_sort。下面说一下introsort的几点要点。

1> 递归层次的限制，超过一定层次直接堆排序剩余数据。

2> 数据长度小于一定值，直接返回，不继续排序。

3>递归结束后，得到多段相互之间有序的数据（段内部无序），再进行一次优化的插入排序。

4>轴元素选取采用三点取中法，以避免分割不均。

sort的代码如下：

1
2
3
4
5
6
7
template <class RandomAccessIterator>
inline void sort(RandomAccessIterator first, RandomAccessIterator last) {
  if (first != last) {
    __introsort_loop(first, last, value_type(first), __lg(last - first) * 2);
    __final_insertion_sort(first, last);
  }
}

从代码中可以看到，sort内部调用的是introsort_loop排序得到多段互相之间有序的小段，再进行插入排序。其中，lg是计算递归层次，定义如下：

1
2
3
4
5
6
7
// 找出 2^k <= n 的最大值k。例，n=7，得k=2，n=20，得k=4，n=8，得k=3。
template <class Size>
inline Size __lg(Size n) {
  Size k;
  for (k = 0; n > 1; n >>= 1) ++k;
  return k;
}

因此，可以看出长度为64的数据，其递归层次不能超过16。如果超过，则进行堆排序。__introsort_loop的代码如下，是整个算法的关键。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
template <class RandomAccessIterator, class T, class Size>
void __introsort_loop(RandomAccessIterator first,
                      RandomAccessIterator last, T*,
                      Size depth_limit) {
  // 以下，__stl_threshold 是個全域常數，稍早定義為 const int 16。
  while (last - first > __stl_threshold) {
    if (depth_limit == 0) {             // 至此，切割惡化
      partial_sort(first, last, last);  // 改用 heapsort
      return;
    }
    --depth_limit;
    // 以下是 median-of-three partition，選擇一個夠好的樞軸並決定切割點。
    // 切割點將落在迭代器 cut 身上。
    RandomAccessIterator cut = __unguarded_partition
      (first, last, T(__median(*first, *(first + (last - first)/2),
                               *(last - 1))));
    // 對右半段遞迴進行 sort.
    __introsort_loop(cut, last, value_type(first), depth_limit);
    last = cut;
    // 現在回到while 迴圈，準備對左半段遞迴進行 sort.
    // 這種寫法可讀性較差，效率並沒有比較好。
    // RW STL 採用一般教科書寫法（直觀地對左半段和右半段遞迴），較易閱讀。
  }
}

从代码里面可以看出，stl里面定义了个常数stl_threshold决定要不要继续循环，如果数据的长度比该值小，那么就不继续递归排序了，除非是调用了堆排，长度小于stl_threshold的段是没有经过排序的。当层次耗尽，depth_limit为0，直接调用partial_sort堆排序返回。另外的部分，就是采用三点取中法决定轴元素。还有这个代码写成了类似尾递归的形式，与教科书上不一致，较难阅读和理解。
__median的代码如下，容易理解。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
template <class T>
inline const T __median(const T a, const T b, const T c) {
  if (a < b)
    if (b < c)      // a < b < c
      return b;
    else if (a < c) // a < b, b >= c, a < c
      return c;
    else
      return a;
  else if (a < c)   // c > a >= b
    return a;
  else if (b < c)       // a >= b, a >= c, b < c
    return c;
  else
    return b;
}

__unguarded_partition与教科书上的partion实现不同，实现简单而且高效。其代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
template <class RandomAccessIterator, class T>
RandomAccessIterator __unguarded_partition(RandomAccessIterator first,
                                           RandomAccessIterator last,
                                           T pivot) {
  while (true) {

    while (*first < pivot) ++first; // first 找到 >= pivot 的元素，就停下來
    --last;                     // 調整
    while (pivot < *last) --last;   // last 找到 <= pivot 的元素，就停下來
    // 注意，以下first < last 判斷動作，只適用於random iterator
    if (!(first < last)) return first;  // 交錯，結束迴圈。
    iter_swap(first, last);             // 大小值交換
    ++first;                            // 調整
  }
}

现在回到sort函数中，最后的插入排序__final_insertion_sort。该函数再进行一次插入排序，由于小段之间都是有序的，因此该函数很快就能执行完。其代码如下：

1
2
3
4
5
6
7
8
9
10
template <class RandomAccessIterator>
void __final_insertion_sort(RandomAccessIterator first,
                            RandomAccessIterator last) {
  if (last - first > __stl_threshold) {
    __insertion_sort(first, first + __stl_threshold);
    __unguarded_insertion_sort(first + __stl_threshold, last);
  }
  else
    __insertion_sort(first, last);
}

该函数判断当前排序范围是不是小于等于于stl_threshold（16），如果是，直接进行优化的插入排序insertion_sort。否则，排序前面的16个元素，再调用unguarded_insertion_sort排序剩下的元素。unguarded_insertion_sort的相关定义代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
template <class RandomAccessIterator>
inline void __unguarded_insertion_sort(RandomAccessIterator first,
                                RandomAccessIterator last) {
  __unguarded_insertion_sort_aux(first, last, value_type(first));
}

template <class RandomAccessIterator, class T>
void __unguarded_insertion_sort_aux(RandomAccessIterator first,
                                    RandomAccessIterator last, T*) {
  for (RandomAccessIterator i = first; i != last; ++i)
    __unguarded_linear_insert(i, T(*i));
}

__insertion_sort的相关代码如下：

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
template <class RandomAccessIterator>
void __insertion_sort(RandomAccessIterator first, RandomAccessIterator last) {
  if (first == last) return;
  for (RandomAccessIterator i = first + 1; i != last; ++i)  // 外迴圈
    __linear_insert(first, i, value_type(first));   // first,i形成一個子範圍
}

template <class RandomAccessIterator, class T>
inline void __linear_insert(RandomAccessIterator first,
                                  RandomAccessIterator last, T*) {
  T value = *last;      // 記錄尾元素
  if (value < *first) { // 尾比頭還小（那就別一個個比較了，一次做完…）
    copy_backward(first, last, last + 1); // 將整個範圍向右遞移一個位置
    *first = value;     // 令頭元素等於原先的尾元素值
  }
  else
    __unguarded_linear_insert(last, value);
}

template <class RandomAccessIterator, class T>
void __unguarded_linear_insert(RandomAccessIterator last, T value) {
  RandomAccessIterator next = last;
  --next;
  // insertion sort 的內迴圈
  // 注意，一旦不出現逆轉對（inversion），迴圈就可以結束了。
  while (value < *next) {   // 逆轉對（inversion）存在
    *last = *next;      // 轉正
    last = next;            // 調整迭代器
    --next;             // 前進一個位置
  }
  *last = value;
}

从上面的代码可以看出，关于插入排序的代码比sort其它部分的代码都要多，可见stl对插入排序进行了特别的优化，以使得其在处理接近有序数据的时候，速度非常快。从代码中可以看到，这些函数最终都调用了__unguarded_linear_insert函数，该函数优化了插入排序的实现。该函数的思路是从后往前找是否存在比value小的元素，如果还存在，就往后移动数据一个位置。最终，得到的是value需要插入的位置。根据stl源码剖析作者的说法，该函数不需要判断越界，因此提高了速度。在大数据量的排序中，优势很大。毕竟，该函数会被sort频繁调用。

关于插入排序代码实现的具体思路，参考stl源码剖析一书。本文的代码也取自该书提供的源码。

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

这里只讲解稀疏矩阵表示的线性系统。最新版本的eigen已经支持稀疏矩阵，以及相关操作了。eigen相对其它库的最大优势，默认情况下，只需要eigen的头文件即可，不需要其它依赖库，除非用对应的加速解法（需要更快的库支持）。所以，eigen就是个牛逼的模板矩阵库。虽然在使用的时候需要套用不少模板参数，有点麻烦。

eigen里面定义稀疏矩阵，然后加入数据的一种方式如下：

1
2
3
4
5
6
7
8
Eigen::SparseMatrix<float> smA(m_nRow, m_nColumn);

std::sort(m_pSI, m_pSI + m_nNum);//排序，列主序
for (int i = 0; i < m_nNum; ++i)// 初始化非零元素
{
    smA.insert(m_pSI[i].r, m_pSI[i].c) = m_pSI[i].v;
}
smA.makeCompressed();

代码非常简单，makeCompressed是可选的，具体原因不怎么清楚，参考官方文档。
知道定义了稀疏矩阵之后，就可以解稀疏线性系统Ax=b了。A作为稀疏系数矩阵输入。eigen里面有多种解线性系统的方法，官方文档有说明。其中，QR分解适合于任何大小的稀疏矩阵A，并且被推荐用来解最小二乘法问题，适用范围很广。下面的是QR分解求线性系统的代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
// 一个QR分解的实例
Eigen::SparseQR<Eigen::SparseMatrix<float>, Eigen::COLAMDOrdering<int>> linearSolver;
// 计算分解
linearSolver.compute(smA);

Eigen::VectorXf vecB(m_nRow);
for (int i = 0; i < m_nRow; ++i)
{
    vecB[i] = pColumnB[i];
}
Eigen::VectorXf vecX = linearSolver.solve(vecB);// 求一个A*x = b

for (int i = 0; i < m_nColumn; ++i)
{
    pX[i] = vecX[i];
}

代码的思路是先定义一个QR分解的线性solver，然后调用compute分解。再传入b，再调用solve函数求解，再传出数据。逻辑很简单。值得注意的是，这里用了float来计算，是因为qr分解比较耗内存，如果矩阵过大或者非0值过多就会计算不过来。因此，用float来减少一半内存，不过精度比double下降了，具体看实际需要选择float还是double。至于QR分解的原理，请参考维基百科等。
如果使用迭代的办法解线性方程组，该怎么用eigen实现了。eigen里面也提供了迭代的solver，不过要求稀疏矩阵A是方阵。那么，把A转换为方阵即可。Ax=b的两边同时乘以A的转置，得到A’Ax=A’b(A’是A的转置)，再用迭代solver求解。如果得到的答案误差较大，可以在两边乘以个（-1）。如下代码所示。

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
//把A*X = b，转换为A'*A*X = A'*b，其中A'是A的转置，貌似计算结果不是很对，故没有使用该函数
    int CSparseMatrix::SolveLinearSystemByEigenIterative(double* pColumnB, double* pX)
    {
        Eigen::SparseMatrix<float> smA(m_nRow, m_nColumn);

        std::sort(m_pSI, m_pSI + m_nNum);//排序，列主序
        for (int i = 0; i < m_nNum; ++i)//初始化非零元素
        {
            smA.insert(m_pSI[i].r, m_pSI[i].c) = m_pSI[i].v;
        }

        Eigen::SparseMatrix<float> smATranspose = smA.transpose();

        Eigen::VectorXf vecB(m_nRow);
        for (int i = 0; i < m_nRow; ++i)
        {
            vecB[i] = pColumnB[i];
        }
        vecB = smATranspose * vecB;//b = A'*b
        vecB = (-1) * vecB;

        smA = smATranspose * smA;//A = A'*A
        smA = (-1) * smA;
        smA.makeCompressed();
        //std::cout << smA << std::endl;

        Eigen::BiCGSTAB<Eigen::SparseMatrix<float>> linearSolver;
        linearSolver.compute(smA);
        Eigen::VectorXf vecX = linearSolver.solve(vecB);// 求一个A*x = b

        std::cout << "#iterations:     " << linearSolver.iterations() << std::endl;
        std::cout << "estimated error: " << linearSolver.error()      << std::endl;

        for (int i = 0; i < m_nColumn; ++i)
        {
            pX[i] = vecX[i];
        }

        return 1;
    }

迭代解法对应的头文件是#include “Eigen/IterativeLinearSolvers”。QR分解求线性系统的完整代码如下，相关定义很容易看懂。

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
//适合于解决最小二乘法问题
    int CSparseMatrix::SolveLinearSystemByEigenQR(double* pColumnB, double* pX)
    {
        Eigen::SparseMatrix<float> smA(m_nRow, m_nColumn);

        std::sort(m_pSI, m_pSI + m_nNum);//排序，列主序
        for (int i = 0; i < m_nNum; ++i)// 初始化非零元素
        {
            smA.insert(m_pSI[i].r, m_pSI[i].c) = m_pSI[i].v;
        }
        smA.makeCompressed();

        //std::cout << smA << std::endl;

        // 一个QR分解的实例
        Eigen::SparseQR<Eigen::SparseMatrix<float>, Eigen::COLAMDOrdering<int>> linearSolver;
        // 计算分解
        linearSolver.compute(smA);

        Eigen::VectorXf vecB(m_nRow);
        for (int i = 0; i < m_nRow; ++i)
        {
            vecB[i] = pColumnB[i];
        }
        Eigen::VectorXf vecX = linearSolver.solve(vecB);// 求一个A*x = b

        for (int i = 0; i < m_nColumn; ++i)
        {
            pX[i] = vecX[i];
        }

        return 1;
    }

QR分解对应的头文件是#include “Eigen/SparseQR”。</pre>

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

以前写过一篇C++调用MATLAB引擎计算特征向量，里面讲了如何配置环境等等。现在又有一个新的需求，求解线性系统。又想起了MATLAB这个工具，于是又写了个小类。这里要求解的是AX=B。其中，A是m*n的矩阵，X是n维行向量，B是m维列向量。MATLAB里面求解很简单X=A\B。

该类的头文件如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
//求解线性系统AX = B.
    //A是Row*Column的矩阵,B是Row*1的列向量
    class CLinearSystem
    {
    public:
        CLinearSystem() {}
        CLinearSystem(double* pMatrixA, int nRow, int nColumn, double* pColumnB);
        ~CLinearSystem() {}

    public:
        void GetResult(double* pX);//pX保证能输出m_nColumn个元素

    private:
        double* m_pMatrixA;
        double* m_pColumnB;
        int m_nRow;
        int m_nColumn;
    };

源文件如下：

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
#include "stdafx.h"
#include "LinearSystem.h"
#include <cassert>
#include "engine.h"
#pragma comment(lib, "libeng.lib")
#pragma comment(lib, "libmx.lib")
#pragma comment(lib, "libmat.lib")

    CLinearSystem::CLinearSystem(double* pMatrixA, int nRow, int nColumn, double* pColumnB)
        : m_pMatrixA(pMatrixA), m_nRow(nRow), m_nColumn(nColumn), m_pColumnB(pColumnB)
    {

    }

    void CLinearSystem::GetResult(double* pX)
    {
        assert(pX != NULL);

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

        mxArray* matrixA = mxCreateDoubleMatrix(m_nRow, m_nColumn, mxREAL);
        //matlab里面的矩阵是列主顺序，这里需要变换位置
        double* pA = mxGetPr(matrixA);
        for (int j = 0; j < m_nColumn; ++j)
        {
            for (int i = 0; i < m_nRow; ++i)
            {
                *pA++ = *(m_pMatrixA + i * m_nColumn + j);
            }
        }

        mxArray* columnB = mxCreateDoubleMatrix(m_nRow, 1, mxREAL);
        memcpy(mxGetPr(columnB), m_pColumnB, sizeof(double) * m_nRow * 1);

        engPutVariable(ep, "A", matrixA);
        engPutVariable(ep, "B", columnB);
        engEvalString(ep, "X = A\\B;");

        mxArray* rowX = engGetVariable(ep, "X"); //返回计算结果
        memcpy(pX, mxGetPr(rowX), sizeof(double) * m_nColumn);

        mxDestroyArray(matrixA);
        mxDestroyArray(columnB);
        mxDestroyArray(rowX);
        engClose(ep);
    }

该类里面就一个计算结果的函数。实现的过程也很简单，打开MATLAB引擎，创建了对应的矩阵和向量，计算表达式。注意，”X = A\B;”中必须加转义\。另外，MATLAB中的矩阵是列主序的，也就是按照列存储的。所以，从C++的矩阵变换到MATLAB矩阵需要重新设置值，而不是简单的拷贝内存。
如果，A是巨大的稀疏矩阵，那么就应该创建稀疏矩阵来计算了。测试代码如下：

1
2
3
4
5
6
7
8
//测试线性系统
int nRow = 3, nColumn = 2;
double pA[3][2] = { {4, 5}, {1, 2}, {3, 1} };
double pB[3] = {3, 15, 12};
CLinearSystem ls((double*)pA, nRow, nColumn, pB);
double pX[2];
ls.GetResult(pX);
printf("%f %f\n", pX[0], pX[1]);

输出为：3.000000 -0.600000

前面说了，系数矩阵的事情，现在代码进行了更新，加入了稀疏矩阵，A可以是稀疏矩阵传入，B的话没有当做稀疏矩阵传入了。不过，在调用MATLAB代码时候还是得作为稀疏矩阵。具体代码如下，就不细说了。

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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
#pragma once

#ifndef LIBYX_LINEAR_SYSTEM_H
#define LIBYX_LINEAR_SYSTEM_H

#include "MatlabSparseMatrix.h"

namespace LibYX
{
    //求解线性系统AX = B.
    //A是Row*Column的矩阵,B是Row*1的列向量
    class CLinearSystem
    {
    public:
        CLinearSystem() {}
        CLinearSystem(double* pMatrixA, int nRow, int nColumn, double* pColumnB);
        CLinearSystem(MatlabSparseInfor* pSI, int nNum, int nRow, int nColumn, double* pColumnB);
        ~CLinearSystem() {}

    public:
        void GetResult(double* pX);//pX保证能输出m_nColumn个元素
        void GetResultSparse(double* pX);//中间过程使用稀疏矩阵计算结果

    private:
        double* m_pMatrixA;
        double* m_pColumnB;
        int m_nRow;
        int m_nColumn;

    private://稀疏矩阵形式输入A
        MatlabSparseInfor* m_pSI;
        int m_nNum;
    };
};

#endif

#include "stdafx.h"
#include "LinearSystem.h"
#include <cassert>
#include "engine.h"
#pragma comment(lib, "libeng.lib") 
#pragma comment(lib, "libmx.lib")
#pragma comment(lib, "libmat.lib")

namespace LibYX
{
    CLinearSystem::CLinearSystem(double* pMatrixA, int nRow, int nColumn, double* pColumnB)
        : m_pMatrixA(pMatrixA), m_nRow(nRow), m_nColumn(nColumn), m_pColumnB(pColumnB)
    {

    }

    CLinearSystem::CLinearSystem(MatlabSparseInfor* pSI, int nNum, int nRow, int nColumn, double* pColumnB)
        : m_pSI(pSI), m_nNum(nNum), m_nRow(nRow), m_nColumn(nColumn), m_pColumnB(pColumnB)
    {

    }

    void CLinearSystem::GetResult(double* pX)
    {
        assert(pX != NULL);

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

        mxArray* matrixA = mxCreateDoubleMatrix(m_nRow, m_nColumn, mxREAL);
        //matlab里面的矩阵是列主顺序，这里需要变换位置
        double* pA = mxGetPr(matrixA);
        for (int j = 0; j < m_nColumn; ++j)
        {
            for (int i = 0; i < m_nRow; ++i)
            {
                *pA++ = *(m_pMatrixA + i * m_nColumn + j);
            }
        }

        mxArray* columnB = mxCreateDoubleMatrix(m_nRow, 1, mxREAL);
        memcpy(mxGetPr(columnB), m_pColumnB, sizeof(double) * m_nRow * 1);

        engPutVariable(ep, "A", matrixA);
        engPutVariable(ep, "B", columnB);
        engEvalString(ep, "X = A\\B;");

        mxArray* rowX = engGetVariable(ep, "X"); //返回计算结果
        memcpy(pX, mxGetPr(rowX), sizeof(double) * m_nColumn);

        mxDestroyArray(matrixA);
        mxDestroyArray(columnB);
        mxDestroyArray(rowX);
        engClose(ep);
    }

    void CLinearSystem::GetResultSparse(double* pX)
    {
        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;

        //以下代码构造稀疏列
        mxArray* columnB = mxCreateSparse(m_nRow, 1, m_nRow, mxREAL);
        pr = mxGetPr(columnB);
        ir = mxGetIr(columnB); 
        jc = mxGetJc(columnB); 

        jc[0] = 0;
        for (int i = 0; i < m_nRow; ++i)
        {
            pr[i] = m_pColumnB[i];
            ir[i] = i;
        }
        jc[1] = m_nRow;

        //将构造的稀疏矩阵传入MATLAB命令中
        engPutVariable(ep, "sm", sm);
        engPutVariable(ep, "columnB", columnB);
        engEvalString(ep, "X = sm\\columnB;");

        mxArray* rowX = engGetVariable(ep, "X"); //返回计算结果
        memcpy(pX, mxGetPr(rowX), sizeof(double) * m_nColumn);

        mxDestroyArray(rowX);
        mxDestroyArray(columnB);
        mxDestroyArray(sm);
        engClose(ep);
    }
};

发表于 分类于 ， 阅读次数： 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);
    }

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