# 等值线追踪算法

### 等值线追踪算法

//矩形4位2进制对应连接线，连接线对应两个边（索引）上的点。
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}
};

#### 追踪代码

struct Dpos{//迭代中用于返回xy值
Dpos(){l=0;r=0;n=0;}
unsigned int l,r;
bool n;//是否有线
};

//等值线值为isoL时的所有线
struct ClineIsoLevel{
float isoL;
vector<vector<unsigned int>> isoLines;
};

typedef std::vector<ClineIsoLevel> CtLineTrack;

class Contour_line{
public:
......

//追踪生成等值线
bool CreatLineTra(float* field,int n[2],vector<float> threshold,vector<float> &vrts,
CtLineTrack& lns);

//从标量场追踪等值线
void TrackLineV(vector<float> &vrts,CtLineTrack& lns);

private:

......

vector<int> m_bTrick;//判断网格上还有几条线
ClineIsoLevel m_allIsoLines;//记录每一个等值线阈值上的点
//追踪递归
Dpos TrackRecursive(unsigned int x,unsigned int y,unsigned int iso,unsigned int edge);

// c输出追踪出的等值线点信息和线信息
void RenameTrackVerticesAndLines(vector<float> &vrts, unsigned int &nvrts,
CtLineTrack& lns);
}

bool Contour_line::CreatLineTra(float* field,int n[2],vector<float> threshold,vector<float> &vrts,
CtLineTrack& 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;

TrackLineV(vrts,lns);//追踪算法
return true;
}

/* 等值线信息生成（等值线追踪的迭代方式）
* @param [out] vrts等值线顶点
* @param [out] lns等值线信息（顶点连接信息）
*/
void Contour_line::TrackLineV(vector<float> &vrts,CtLineTrack& lns){
if (m_bValidLine)
DeleteLine();

//等值线生成
for (unsigned int iso=0; iso<m_tIsoLevel.size(); iso++) {
ClineIsoLevel currentIsoL;//每次保存该值的所有等值线
currentIsoL.isoL=m_tIsoLevel[iso];
m_bTrick.clear();
m_bTrick.resize(m_Grid[0]*m_Grid[1], 1);//每次更新表格初始化为1
for (unsigned int y=0; y<m_Grid[1]; y++) {
for (unsigned int x=0; x<m_Grid[0]; x++) {
if (m_bTrick[y*m_Grid[0]+x]) {
vector<unsigned int>currentL,currentR;//保存其中一条
unsigned int nx=x,ny=y,lastL,lastR;
Dpos np=TrackRecursive(nx,ny,iso,0),fp=np;
if (np.n) {
lastL=np.l;
lastR=np.r;
while (np.n) {//处理左边线
unsigned int id;
switch (lastL) {//查表判断下一个网格遍历位置
case 1:
id=GetEdgeID(nx, ny, 1);
np=TrackRecursive(nx, ny-=1, iso,posTable[1]);
break;
case 2:
id=GetEdgeID(nx, ny, 2);
np=TrackRecursive(nx+=1, ny, iso,posTable[2]);
break;
case 3:
id=GetEdgeID(nx, ny, 3);
np=TrackRecursive(nx, ny+=1, iso,posTable[3]);
break;
case 4:
id=GetEdgeID(nx, ny, 4);
np=TrackRecursive(nx-=1, ny, iso,posTable[4]);
break;
}
currentL.push_back(id);
lastL=np.r;
}
np=fp; nx=x;ny=y;
while (np.n) {//处理右边线
unsigned int id;
switch (lastR) {
case 1:
id=GetEdgeID(nx, ny, 1);
np=TrackRecursive(nx, ny-=1, iso,posTable[1]);
break;
case 2:
id=GetEdgeID(nx, ny, 2);
np=TrackRecursive(nx+=1, ny, iso,posTable[2]);
break;
case 3:
id=GetEdgeID(nx, ny, 3);
np=TrackRecursive(nx, ny+=1, iso,posTable[3]);
break;
case 4:
id=GetEdgeID(nx, ny, 4);
np=TrackRecursive(nx-=1, ny, iso,posTable[4]);
break;
}
currentR.push_back(id);
lastR=np.r;
}
for (int i=0; i<(currentR.size()+1)/2; i++) {
unsigned int temp=currentR[i];
currentR[i]=currentR[currentR.size()-i-1];
currentR[currentR.size()-i-1]=temp;
}
for (int i=0; i<currentL.size(); i++) {
currentR.push_back(currentL[i]);
}
currentIsoL.isoLines.push_back(currentR);
}
}
}
}
m_allIsoLines=currentIsoL;
RenameTrackVerticesAndLines(vrts, m_nVertices, lns);
}
m_bValidLine=true;
}

/* 追踪算法
* @param [out] vrts等值线顶点
* @param [out] lns等值线信息（顶点连接信息）
*/
Dpos Contour_line::TrackRecursive(unsigned int x,unsigned int y,unsigned int iso,unsigned int edge){
Dpos apos;
if (x<m_Grid[0]&&y<m_Grid[1]&&m_bTrick[y*m_Grid[0]+x]) {
//计算网格内的顶点放置信息索引
unsigned int tableIndex=0;
int stx=m_stride*x,sty=m_stride*y;

//存在空值情况
int nanCount=0,nanNum=0;
if (m_ptScalarField[6*(sty*m_sliceX+stx)+2]==-1.0){
nanCount++;
nanNum=1;
}
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++;
nanNum=2;
}
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++;
nanNum=3;
}
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++;
nanNum=4;
}
else if (m_ptScalarField[6*((sty+m_stride)*m_sliceX+stx)+2]<m_tIsoLevel[iso])
tableIndex|=8;

//处理单nan特殊情况
if (nanCount==1) {
int mid=pow(2, nanNum-1);
int left=(nanNum-1)<1?4:(nanNum-1);
left=pow(2, left-1);
int right=(nanNum+1)>4?1:(nanNum+1);
right=pow(2, right-1);
if ((tableIndex&left)==(tableIndex&right)) {
tableIndex|=(tableIndex&left)&mid;
nanCount=0;
}
}

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

if (!nanCount) {
//生成等值线,获取线上的点在哪条边上
apos.n=1;
if(edge!=0){//当非第一次遍历时在r中存放接下来的点
for (unsigned int i=0; lineTable[tableIndex][i]!=255; i+=2) {
unsigned int pointID1,pointID2;
pointID1=lineTable[tableIndex][i];
pointID2=lineTable[tableIndex][i+1];
if (pointID1==edge)
apos.r=pointID2;//用r表示下个点
else if(pointID2==edge)
apos.r=pointID1;
}
}else{
apos.l=lineTable[tableIndex][0];
apos.r=lineTable[tableIndex][1];
}
if(lineTable[tableIndex][2]!=255){//处理一个网格有两条线情况，0，1，2三个值
if (m_bTrick[y*m_Grid[0]+x]==1)
m_bTrick[y*m_Grid[0]+x]=2;
else
m_bTrick[y*m_Grid[0]+x]=0;
}
else
m_bTrick[y*m_Grid[0]+x]=0;
return apos;
}
}
m_bTrick[y*m_Grid[0]+x]=0;
}
apos.n=0;
return  apos;
}

/*！
* 复制信息到顶点和线缓存中
* @param [out] vrts顶点坐标
* @param [out] nvrts顶点数
* @param [out] lns等值线的几何信息
*/
void Contour_line::RenameTrackVerticesAndLines(vector<float> &vrts, unsigned int &nvrts,
CtLineTrack& lns){
static unsigned int nextID=0;
auto mapIt=m_i2v.begin();

//为每个顶点标记在顶点缓存中的编号
while (mapIt!=m_i2v.end()) {
(*mapIt).second.newID=nextID;
nextID++;
mapIt++;
}

//将等值线m_allIsoLines内的顶点编号更新为顶点缓存中的编号
for (int j=0; j<m_allIsoLines.isoLines.size(); j++) {//每个阈值下的每条等值线
for (int k=0; k<m_allIsoLines.isoLines[j].size(); k++) {//每条等值线上的点
unsigned int newID=m_i2v[m_allIsoLines.isoLines[j][k]].newID;
m_allIsoLines.isoLines[j][k]=newID;
}
}

//复制所有的顶点到顶点缓存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);
}

lns.push_back(m_allIsoLines);

//释放空间
m_i2v.clear();
m_lineVertex.clear();
}


• 4
点赞
• 25
收藏
觉得还不错? 一键收藏
• 1
评论
03-16 1772
04-08 5508
05-26
04-24 302
07-19 695
02-28 749
01-07
08-22 573
10-14
02-21
04-14

### “相关推荐”对你有帮助么？

• 非常没帮助
• 没帮助
• 一般
• 有帮助
• 非常有帮助

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。