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());
    }

等值线效果图：

相关文章