最近导师提了一个实现等值线算法的要求。在思考要怎么写的时候，突然想到表面重建就是重建的等值面，那来个降维打击，不就是等值线了嘛！
网上以及许多文章里的等值线算法都是基于三角形网格的算法，由于我项目的数据性质，我就使用正方形的网格算法，主要是方便编写哈哈。
等值线算法主要步骤分为以下几步：
构建网格，并为每个网格分配索引。根据每个网格索引去索引边表，获取其所在边位置（并赋予边编号），并记录点坐标和边编号。通过每个网格索引去索引线表，索引出点的连接顺序（连接顺序形式为边编号），用边编号从记录的信息中索引出点的坐标位置。最后根据所有点的坐标绘制连线。如该图所示，我们构建一个N*M的网格，网格上每个黑点都有一个数据值（降水量）。我们令黑点上的值大于阈值（等值线值）时为1，小于阈值时为0，我们可以保存在一个4位2进制内。我们设定一个顶点顺序(白色数字)和边顺序(绿色数字)我们：
根据该顺序，我们保存4位2进制根据如下图：
顶点和边的4位二进组都根据上图保存，得到的数最后转化为10进制即可。我们根据此便可以画出所有的情况：
如图所示，我们会发现所有的情况为以上16种。比如我们的1，2，4号顶点标记为0。3号顶点标记为1。我们能获得4位2进制数0010(10进制为2)，我们发现该情况边会连接在第三号边和第四号边上，即：
可以得到新的4位2进制0011(10进制为3)，因为连接线处于三号边和4号边上。所以我们会得到edgeTable[2]=3。我们列出所有的16种情况就可以建立边表：
//矩形4位2进制对应边索引。如0010——>0011，然后转为10进制 const unsigned char edgeTable[16]={ 0, 9, 3, 10, 6, 15, 5, 12, 12, 5, 15, 6, 10, 3, 9, 0 };我们还可以同样的建立一个线表，表示该网格连接的线为哪两条边(由于有可能出现2条线的情况，所以我们第二维用大小为5的数组，便于之后的遍历)。
//矩形4位2进制对应连接线，连接线对应两个边（索引）上的点。255指边上无点 const unsigned char lineTable[16][5]={ {255,255,255,255,255}, {4,1,255,255,255}, {2,1,255,255,255}, {4,2,255,255,255}, {3,2,255,255,255}, {4,3,2,1,255}, {3,1,255,255,255}, {4,3,255,255,255}, {4,3,255,255,255}, {3,1,255,255,255}, {3,2,1,4,255}, {3,2,255,255,255}, {4,2,255,255,255}, {2,1,255,255,255}, {4,1,255,255,255}, {255,255,255,255,255} };网格构建完毕之后，分配索引就很简单了，我们只要遍历网格，判断其是否大于等值线阈值然后赋值即可。代码如下：
//等值线生成 for (int iso=0; iso<m_tIsoLevel.size(); iso++) {//多条等值线情况 for (unsigned int y=0; y<m_Grid[1]; y++) { for (unsigned int x=0; x<m_Grid[0]; x++) { //计算网格内的顶点放置信息索引 unsigned int tableIndex=0; int stx=m_stride*x,sty=m_stride*y;//stride为绘制网格的缩放值，为1就是1:1 if (m_ptScalarField[6*(sty*m_sliceX+stx)+2]<m_tIsoLevel[iso]) tableIndex|=1; if (m_ptScalarField[6*(sty*m_sliceX+stx+m_stride)+2]<m_tIsoLevel[iso]) tableIndex|=2; if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx+m_stride)+2]<m_tIsoLevel[iso]) tableIndex|=4; if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx)+2]<m_tIsoLevel[iso]) tableIndex|=8; ... } } }
我们已经获得了所有的网格索引，可以通过网格索引判断顶点落在哪个边上了。我们为所有边编号：
通过网格索引获取边信息，然后转化为所有边的编号，然后获取边的两个顶点，进行插值计算点的位置，我们就保存该点的位置信息和边编号。
通过网格索引，获取连线的两个顶点信息，通过在边索引表阶段保存的点信息和编号，保存连线的最终点信息。绘制效果如下图：
若我们选择等值线阈值为4，则可生成图示效果。
下面贴出源码，.h文件：
#include <stdio.h> #include <glm/glm.hpp> #include <glm/gtc/matrix_transform.hpp> #include <glm/gtc/type_ptr.hpp> #include <map> #include <string> #include <vector> #include <iterator> #include "Ct_table.hpp" using namespace std; struct CtVertexID{ //点的信息 unsigned int newID;//索引 double x,y;//位置 }; struct CtLine{ //线的信息 unsigned int pointID[2];//存有每个点的索引 }; typedef std::map<unsigned int,CtVertexID> ID2VertexID;//边序号to点信息 typedef std::vector<CtLine> CtLineV;//保存所有的线的vector class Contour_line{ public: Contour_line(); ~Contour_line(); //从标量场生成等值线(vrts保存所有点，lns保存所有的线上点点索引) bool CreatLineV(float* field,int n[2],vector<float> threshold,vector<float>& vrts, vector<int> &lns); //从标量场生成直线段 void GenerateLineV(vector<float>& vrts, vector<int> &lns); //等值线创建成功则返回true bool IsLineValid() const { return m_bValidLine; } //删除等值线 void DeleteLine(); //获取等值线有关的信息 unsigned int GetNumVertices() const { return m_nVertices; } unsigned int GetNumLines() const { return m_nLines; } void SetStride(int s){m_stride=s;} private: unsigned int m_nVertices,m_nLines;//顶点个数与线段个数 ID2VertexID m_i2v;//形成等值线的顶点列表(以边索引为key，等值面的点为value) CtLineV m_lineVertex;//形成等值线的顶点列表 int m_Grid[2];//网格大小x*y,x:[0],y:[1] const float * m_ptScalarField; //保存标量场的样本值 vector<float> m_tIsoLevel; //阈值 bool m_bValidLine; //表面是否生成成功 int m_stride;//网格偏移量 int m_sliceX;//数据的列数 // 边id unsigned int GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nEdgeNo); // 顶点ID unsigned int GetVertexID(unsigned int nX, unsigned int nY); // 计算边上的等值点 CtVertexID CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nEdgeNo,int isoI); // 通过边两端的隐式函数值的线性插值计算等值点 CtVertexID Interpolate(double fX1, double fY1, double fX2, double fY2, float tVal1, float tVal2,int isoI); // 以输出形式存储顶点和等值线几何信息 void RenameVerticesAndLines(vector<float> &vrts, unsigned int &nvrts, vector<int> &lns, unsigned int &nlns); };
.cpp文件：
#include "contour.hpp" Contour_line::Contour_line(){ m_Grid[0]=0; m_Grid[1]=0; m_stride=1; m_nVertices=0; m_nLines=0; m_ptScalarField=nullptr; m_bValidLine=false; } Contour_line::~Contour_line(){ DeleteLine(); } void Contour_line::DeleteLine(){ m_Grid[0]=0; m_Grid[1]=0; m_nVertices=0; m_nLines=0; m_ptScalarField=nullptr; m_tIsoLevel.clear(); m_bValidLine=false; } /* 等值线生成（标量场创建） * @param [in] field标量场信息 * @param [in] threshold等值线阈值 * @param [out] vrts等值线顶点 * @param [out] lns等值线信息（顶点连接信息） */ bool Contour_line::CreatLineV(float* field,int n[2],vector<float> threshold,vector<float>& vrts, vector<int>&lns){ if (field==nullptr) return false; m_sliceX=n[0]; m_Grid[0]=(n[0]-1)/m_stride; m_Grid[1]=(n[1]-1)/m_stride; m_tIsoLevel=threshold; m_ptScalarField=field; //生成等值线上的点和线 GenerateLineV(vrts, lns); return true; } /* 等值线信息生成（标量场创建） * @param [in] field标量场信息 * @param [in] threshold等值线阈值 * @param [out] vrts等值线顶点 * @param [out] lns等值线信息（顶点连接信息） */ void Contour_line::GenerateLineV(vector<float> &vrts, vector<int> &lns){ if (m_bValidLine) { DeleteLine(); } //等值线生成 for (int iso=0; iso<m_tIsoLevel.size(); iso++) { for (unsigned int y=0; y<m_Grid[1]; y++) { for (unsigned int x=0; x<m_Grid[0]; x++) { //计算网格内的顶点放置信息索引 unsigned int tableIndex=0; int stx=m_stride*x,sty=m_stride*y; int nanCount=0; if (m_ptScalarField[6*(sty*m_sliceX+stx)+2]==-1.0){ nanCount++; } else if (m_ptScalarField[6*(sty*m_sliceX+stx)+2]<m_tIsoLevel[iso]) tableIndex|=1; if (m_ptScalarField[6*(sty*m_sliceX+stx+m_stride)+2]==-1.0){ nanCount++; } else if (m_ptScalarField[6*(sty*m_sliceX+stx+m_stride)+2]<m_tIsoLevel[iso]) tableIndex|=2; if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx+m_stride)+2]==-1.0){ nanCount++; } else if(m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx+m_stride)+2]<m_tIsoLevel[iso]) tableIndex|=4; if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx)+2]==-1.0){ nanCount++; } else if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx)+2]<m_tIsoLevel[iso]) tableIndex|=8; if (edgeTable[tableIndex]!=0) { //计算边上的顶点 if (edgeTable[tableIndex]&1) { CtVertexID pt=CalculateIntersection(stx, sty, 1,iso); unsigned int id=GetEdgeID(x, y, 1); m_i2v.insert(ID2VertexID::value_type(id,pt)); } if (edgeTable[tableIndex]&8) { CtVertexID pt=CalculateIntersection(stx, sty, 4,iso); unsigned int id=GetEdgeID(x, y, 4); m_i2v.insert(ID2VertexID::value_type(id,pt)); } if(x==m_Grid[0]-1){ if (edgeTable[tableIndex]&2) { CtVertexID pt=CalculateIntersection(stx, sty, 2,iso); unsigned int id=GetEdgeID(x, y, 2); m_i2v.insert(ID2VertexID::value_type(id,pt)); } } if(y==m_Grid[1]-1){ if (edgeTable[tableIndex]&4) { CtVertexID pt=CalculateIntersection(stx, sty, 3,iso); unsigned int id=GetEdgeID(x, y, 3); m_i2v.insert(ID2VertexID::value_type(id,pt)); } } //生成等值线,获取线上的点在哪条边上 for(unsigned int i=0;lineTable[tableIndex][i]!=255;i+=2){ CtLine line; unsigned int pointID1,pointID2; pointID1=GetEdgeID(x, y, lineTable[tableIndex][i]); pointID2=GetEdgeID(x, y, lineTable[tableIndex][i+1]); line.pointID[0]=pointID1; line.pointID[1]=pointID2; m_lineVertex.push_back(line); } } } } RenameVerticesAndLines(vrts, m_nVertices, lns, m_nLines); } m_bValidLine=true; } /* 获取边索引 * @param [in] nX当前网格x值 * @param [in] nY当前网格y值 * @param [in] nEdgeNo网格索引（1-4） */ unsigned int Contour_line::GetEdgeID(unsigned int nX, unsigned int nY, unsigned int nEdgeNo){ switch (nEdgeNo) { case 1: return GetVertexID(nX, nY); break; case 2: return GetVertexID(nX+1, nY)+1; break; case 3: return GetVertexID(nX, nY+1); break; case 4: return GetVertexID(nX, nY)+1; break; default: //无效值 return -1; break; } } /*！ *获取顶点所在边ID * @param [in] nX，nY网格位置 * @return顶点所在边ID */ unsigned int Contour_line::GetVertexID(unsigned int nX, unsigned int nY){ return 2*(nY*(m_Grid[0]+1)+nX); } /*！ * 通过插值计算边缘上的等值点（从样本量） * @param [in] nX，nY网格位置 * @param [in] nEdgeNo边数 * @return网格顶点信息 */ CtVertexID Contour_line::CalculateIntersection(unsigned int nX, unsigned int nY, unsigned int nEdgeNo,int isoI){ float x1, y1, x2, y2; unsigned int v1x = nX, v1y = nY; unsigned int v2x = nX, v2y = nY; switch (nEdgeNo) { case 1: v2x+=m_stride; break; case 2: v1x+=m_stride; v2x+=m_stride; v2y+=m_stride; break; case 3: v2x+=m_stride; v1y+=m_stride; v2y+=m_stride; break; case 4: v2y+=m_stride; break; default: break; } //获取边的两点坐标 x1=m_ptScalarField[6*(v1y*m_sliceX+v1x)]; y1=m_ptScalarField[6*(v1y*m_sliceX+v1x)+1]; x2=m_ptScalarField[6*(v2y*m_sliceX+v2x)]; y2=m_ptScalarField[6*(v2y*m_sliceX+v2x)+1]; float val1=m_ptScalarField[6*(v1y*m_sliceX+v1x)+2]; float val2=m_ptScalarField[6*(v2y*m_sliceX+v2x)+2]; CtVertexID intersection=Interpolate(x1, y1, x2, y2, val1, val2, isoI); return intersection; } /*！ *通过网格边缘两端隐式函数值的线性插值计算等值点 * @param [in] fX 1，fY 1终点坐标1 * @param [in] fX 2，fY 2终点坐标2 * @param [in] tVal 1在终点坐标1处的标量值 * @param [in] tVal 2在终点坐标2处的标量值 * @返回顶点信息 */ CtVertexID Contour_line::Interpolate(double fX1, double fY1, double fX2, double fY2, float tVal1, float tVal2,int isoI){ CtVertexID interp; float mu=float(m_tIsoLevel[isoI]-tVal1)/(tVal2-tVal1); interp.x=fX1+mu*(fX2-fX1); interp.y=fY1+mu*(fY2-fY1); return interp; } /*！ * 复制信息到顶点和线缓存中 * @param [out] vrts顶点坐标 * @param [out] nvrts顶点数 * @param [out] lns等值线的几何信息 * @param [out] nlns等值线的数量 */ void Contour_line::RenameVerticesAndLines(vector<float> &vrts, unsigned int &nvrts, vector<int> &lns, unsigned int &nlns){ static unsigned int nextID=0; auto mapIt=m_i2v.begin(); auto vecIt=m_lineVertex.begin(); //为每个顶点标记在顶点缓存中的编号 while (mapIt!=m_i2v.end()) { (*mapIt).second.newID=nextID; nextID++; mapIt++; } //将等值线vector内的顶点编号更新为顶点缓存中的编号 while (vecIt!=m_lineVertex.end()) { for (unsigned int i=0; i<2; i++) { unsigned int newID=m_i2v[(*vecIt).pointID[i]].newID; (*vecIt).pointID[i]=newID; } vecIt++; } //复制所有的顶点到顶点缓存vrts中 mapIt=m_i2v.begin(); nvrts=int(m_i2v.size()); for (unsigned int i=0; i<nvrts; i++,mapIt++) { vrts.push_back((*mapIt).second.x); vrts.push_back((*mapIt).second.y); } //复制所有的线索引 vecIt=m_lineVertex.begin(); nlns=int(m_lineVertex.size()); for (unsigned int i=0; i<nlns; i++,vecIt++) { lns.push_back((*vecIt).pointID[0]); lns.push_back((*vecIt).pointID[1]); } //释放空间 m_i2v.clear(); m_lineVertex.clear(); }最后我根据我的数据绘制了总体效果图为：
由于数据中有许多nan（空）值，所以这里我就直接在空值不去绘制线条。
