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

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

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

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

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

第一步，求出三个顶点的标量值顺序，那么标量范围在[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的半边结构，同时该函数标记了当前面是否已经计算过等值线，避免等值线过密。

	bool CSEMesh::ComputeIsoLine(std::vector& scalarField, Enriched_Model::Facet_iterator pBegFacet,
		std::vector& bHasIsoLine, IsoLine& isoline)
	{
		//该面已经有通过的isoline
		if (bHasIsoLine[pBegFacet->tag()] == true)
		{
			return false;
		}

		isoline.pts.clear();
		std::vector 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 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;
	}

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

//计算所有等值线
	void CSEMesh::ComputeIsoLines(std::vector& scalarField)
	{
		IsoLine isoline;
		std::vector 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 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 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());
	}

等值线效果图：

您可能也喜欢：
使用CGAL实现平面分割三维模型	在三维模型上可视化标量场及注意方面	使用CGAL计算SDF分割模型	PaintMeshCutting的实现	hdu 3584 Cube 三维树状数组
无觅