提取三维表面标量场的等值线

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

那么,如何提取等值线了?

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

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

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

第一步,求出三个顶点的标量值顺序,那么标量范围在[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());
}

等值线效果图: