476 lines
11 KiB
C
476 lines
11 KiB
C
#ifndef _DATASTRUCT_H
|
||
#define _DATASTRUCT_H
|
||
#include "sysDefine.h"
|
||
|
||
//直角坐标系下的一个点
|
||
struct cpoint
|
||
{
|
||
double x,y,z;
|
||
cpoint() //初始化坐标值
|
||
{
|
||
x = y = z = MAX_DBL;
|
||
}
|
||
};
|
||
typedef vector<cpoint> cpointArray;
|
||
|
||
//直角坐标点的一些数学运算
|
||
cpoint operator -(cpoint a, cpoint b)
|
||
{
|
||
cpoint m;
|
||
m.x=a.x-b.x;
|
||
m.y=a.y-b.y;
|
||
m.z=a.z-b.z;
|
||
return m;
|
||
}
|
||
|
||
cpoint operator +(cpoint a, cpoint b) //矢量加法
|
||
{
|
||
cpoint m;
|
||
m.x=a.x+b.x;
|
||
m.y=a.y+b.y;
|
||
m.z=a.z+b.z;
|
||
return m;
|
||
}
|
||
|
||
cpoint operator *(double sign,cpoint b) //矢量乘法
|
||
{
|
||
cpoint m;
|
||
m.x=sign*b.x;
|
||
m.y=sign*b.y;
|
||
m.z=sign*b.z;
|
||
return m;
|
||
}
|
||
|
||
//重载逻辑等操作符作用于矢量，判断两个直角点是否相等
|
||
bool operator ==(cpoint a, cpoint b)
|
||
{
|
||
if(fabs(a.x-b.x)<ZERO&&fabs(a.y-b.y)<ZERO&&fabs(a.z-b.z)<ZERO)
|
||
{
|
||
return 1;
|
||
}
|
||
else return 0;
|
||
}
|
||
|
||
double dot(cpoint a, cpoint b) //矢量点乘
|
||
{
|
||
return a.x*b.x+a.y*b.y+a.z*b.z;
|
||
}
|
||
|
||
cpoint cross(cpoint a,cpoint b) //矢量叉乘
|
||
{
|
||
cpoint v;
|
||
v.x = a.y*b.z-a.z*b.y;
|
||
v.y = a.z*b.x-a.x*b.z;
|
||
v.z = a.x*b.y-a.y*b.x;
|
||
return v;
|
||
}
|
||
|
||
//返回两个直角坐标点的中点位置
|
||
cpoint middle_cpoint(cpoint a,cpoint b)
|
||
{
|
||
cpoint c;
|
||
c.x = 0.5*(a.x + b.x);
|
||
c.y = 0.5*(a.y + b.y);
|
||
c.z = 0.5*(a.z + b.z);
|
||
return c;
|
||
}
|
||
|
||
//返回两点之间的一个点 以第一个点为参考点 第三个参数为相对于原线段的比例
|
||
cpoint scale_cpoint(cpoint a,cpoint b,double scale)
|
||
{
|
||
cpoint c;
|
||
c.x = a.x + (b.x - a.x)*scale;
|
||
c.y = a.y + (b.y - a.y)*scale;
|
||
c.z = a.z + (b.z - a.z)*scale;
|
||
return c;
|
||
}
|
||
|
||
cpoint rescale_cpoint(cpoint a,double refr)
|
||
{
|
||
cpoint c;
|
||
double m = sqrt(a.x*a.x+a.y*a.y+a.z*a.z);
|
||
c.x = a.x*refr/m;
|
||
c.y = a.y*refr/m;
|
||
c.z = a.z*refr/m;
|
||
return c;
|
||
}
|
||
|
||
double length_cpoint(cpoint v) //矢量模
|
||
{
|
||
return sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
|
||
}
|
||
|
||
double distance_cpoint(cpoint a, cpoint b)
|
||
{
|
||
cpoint m;
|
||
double d;
|
||
m.x=a.x-b.x;
|
||
m.y=a.y-b.y;
|
||
m.z=a.z-b.z;
|
||
d = sqrt(m.x*m.x + m.y*m.y + m.z*m.z);
|
||
return d;
|
||
}
|
||
|
||
//计算两个向量的夹角
|
||
double cpoint_angle(cpoint a,cpoint b)
|
||
{
|
||
return acos((a.x*b.x+a.y*b.y+a.z*b.z)/(sqrt(a.x*a.x+a.y*a.y+a.z*a.z)*sqrt(b.x*b.x+b.y*b.y+b.z*b.z)))*180.0/pi;
|
||
}
|
||
|
||
//求三角形中心坐标
|
||
cpoint Tri_center(cpoint vec1,cpoint vec2,cpoint vec3)
|
||
{
|
||
cpoint c;
|
||
c.x = (vec1.x+vec2.x+vec3.x)/3.0;
|
||
c.y = (vec1.y+vec2.y+vec3.y)/3.0;
|
||
c.z = (vec1.z+vec2.z+vec3.z)/3.0;
|
||
return c;
|
||
}
|
||
|
||
//球坐标系下的一个点
|
||
struct spoint
|
||
{
|
||
double lon,lat,rad;
|
||
spoint() //初始化坐标值
|
||
{
|
||
lon = lat = rad = MAX_DBL;
|
||
}
|
||
};
|
||
typedef vector<spoint> spointArray;
|
||
|
||
/*直角坐标与球坐标相互转换函数 注意这里使用的球坐标是地理坐标范围 即经度为-180~180 纬度为-90~90*/
|
||
cpoint s2c(spoint s)
|
||
{
|
||
cpoint c;
|
||
c.x = s.rad*sin((0.5 - s.lat/180.0)*pi)*cos((2.0 + s.lon/180.0)*pi);
|
||
c.y = s.rad*sin((0.5 - s.lat/180.0)*pi)*sin((2.0 + s.lon/180.0)*pi);
|
||
c.z = s.rad*cos((0.5 - s.lat/180.0)*pi);
|
||
return c;
|
||
}
|
||
|
||
spoint c2s(cpoint c)
|
||
{
|
||
spoint s;
|
||
s.rad = sqrt(pow(c.x,2)+pow(c.y,2)+pow(c.z,2));
|
||
if (fabs(s.rad)<ZERO) //点距离原点极近 将点置于原点
|
||
{
|
||
s.lat = s.lon = s.rad = 0.0;
|
||
}
|
||
else
|
||
{
|
||
s.lat = 90.0 - acos(c.z/s.rad)*180.0/pi;
|
||
s.lon = atan2(c.y,c.x)*180.0/pi;
|
||
}
|
||
return s;
|
||
}
|
||
|
||
//顶点
|
||
struct vertex
|
||
{
|
||
int id; //索引
|
||
cpoint posic; //直角坐标系位置
|
||
spoint posis; //球坐标系位置
|
||
vertex()
|
||
{
|
||
id = -1; //初始化顶点索引值为-1 这里不需要初始化坐标位置 因为已经由相应的初始化函数完成了初始化
|
||
}
|
||
void set(int i) //设置索引值
|
||
{
|
||
id = i;
|
||
}
|
||
void set(cpoint c) //从直角坐标位置初始化
|
||
{
|
||
posic.x = c.x; posic.y = c.y; posic.z = c.z;
|
||
posis = c2s(posic);
|
||
}
|
||
void set(spoint s) //从球坐标位置初始化
|
||
{
|
||
posis.lon = s.lon; posis.lat = s.lat; posis.rad = s.rad;
|
||
posic = s2c(posis);
|
||
}
|
||
void info() //显示顶点信息
|
||
{
|
||
cout << id << " " << setprecision(16) << posic.x << " " << posic.y << " " << posic.z << " " << posis.lon << " " << posis.lat << " " << posis.rad << endl;
|
||
}
|
||
};
|
||
typedef vector<vertex> vertexArray;
|
||
typedef map<int,vertex> idMap; //顶点索引值映射 用于通过索引值寻找相应顶点
|
||
typedef map<string,vertex> strMap; //顶点位置映射 用于通过顶点位置寻找相应顶点
|
||
typedef map<int,int> outIdMap; //输出msh文件时重新索引三角形顶点集
|
||
|
||
//计算一个顶点向量的中点
|
||
cpoint middle_vertex(vertexArray vert)
|
||
{
|
||
cpoint c;
|
||
c.x = 0; c.y = 0; c.z = 0;
|
||
if (!vert.empty())
|
||
{
|
||
for (int i = 0; i < vert.size(); i++)
|
||
{
|
||
c.x += vert.at(i).posic.x;
|
||
c.y += vert.at(i).posic.y;
|
||
c.z += vert.at(i).posic.z;
|
||
}
|
||
c.x /= vert.size();
|
||
c.y /= vert.size();
|
||
c.z /= vert.size();
|
||
}
|
||
return c;
|
||
}
|
||
|
||
/*旋转顶点的方位
|
||
旋转的方式为首先绕着x轴旋转 然后绕着z轴旋转
|
||
给出参考点的原位置olda与新位置newa以获取旋转参数 然后求取输入点oldb做相同旋转后的新坐标newb
|
||
*/
|
||
|
||
vertex rotate_vertex(vertex olda,vertex newa,vertex oldb)
|
||
{
|
||
vertex newb;
|
||
vertex temp_ref,temp_b;
|
||
double yz_angle = (newa.posis.lat - olda.posis.lat)*pi/180.0;
|
||
//首先绕olda.lon即x轴旋转oldb
|
||
temp_b.posic.x = oldb.posic.x;
|
||
temp_b.posic.y = oldb.posic.y*cos(-1.0*yz_angle)+oldb.posic.z*sin(-1.0*yz_angle);
|
||
temp_b.posic.z = oldb.posic.z*cos(-1.0*yz_angle)-oldb.posic.y*sin(-1.0*yz_angle);
|
||
temp_b.set(temp_b.posic);
|
||
//计算绕x轴旋转后olda的位置 这是后一步旋转需要的参考值
|
||
temp_ref.posic.x = olda.posic.x;
|
||
temp_ref.posic.y = olda.posic.y*cos(-1.0*yz_angle)+olda.posic.z*sin(-1.0*yz_angle);
|
||
temp_ref.posic.z = olda.posic.z*cos(-1.0*yz_angle)-olda.posic.y*sin(-1.0*yz_angle);
|
||
temp_ref.set(temp_ref.posic);
|
||
//注意绕z轴旋转的经度参考位置为olda绕x轴旋转后的经度值
|
||
double xy_angle = (newa.posis.lon - temp_ref.posis.lon)*pi/180.0;
|
||
//绕z轴旋转temp_b z值不变
|
||
newb.id = oldb.id;
|
||
newb.posic.x = temp_b.posic.x*cos(-1.0*xy_angle)+temp_b.posic.y*sin(-1.0*xy_angle);
|
||
newb.posic.y = temp_b.posic.y*cos(-1.0*xy_angle)-temp_b.posic.x*sin(-1.0*xy_angle);
|
||
newb.posic.z = temp_b.posic.z;
|
||
newb.set(newb.posic);
|
||
return newb;
|
||
}
|
||
|
||
//点 点包含索引和一个顶点
|
||
struct point
|
||
{
|
||
int id;
|
||
int maxDepth;
|
||
double minDeg;
|
||
int physic;
|
||
vertex vert;
|
||
point()
|
||
{
|
||
id = -1;
|
||
maxDepth = -1;
|
||
minDeg = -1.0;
|
||
}
|
||
void info()
|
||
{
|
||
cout << id << " " << maxDepth << " " << minDeg << endl;
|
||
vert.info();
|
||
}
|
||
};
|
||
typedef vector<point> pointArray;
|
||
|
||
//折线 折线包含索引和一个顶点向量 顶点从前向后连成折线
|
||
struct line
|
||
{
|
||
int id;
|
||
int maxDepth;
|
||
double minDeg;
|
||
int physic;
|
||
vertexArray vert;
|
||
line()
|
||
{
|
||
id = -1;
|
||
maxDepth = -1;
|
||
minDeg = -1.0;
|
||
}
|
||
void info()
|
||
{
|
||
cout << id << " " << maxDepth << " " << minDeg << endl;
|
||
for (int i = 0; i < vert.size(); i++)
|
||
{
|
||
vert.at(i).info();
|
||
}
|
||
}
|
||
void clear_vert()
|
||
{
|
||
if (!vert.empty()) vert.clear();
|
||
}
|
||
};
|
||
typedef vector<line> lineArray;
|
||
|
||
//多边形 多边形包含索引和一个顶点向量 顶点逆时针连成多边形 注意多边形第一个点和最后一个点应该一致
|
||
struct polygon
|
||
{
|
||
int id;
|
||
int maxDepth;
|
||
double minDeg;
|
||
int physic;
|
||
vertexArray vert;
|
||
polygon()
|
||
{
|
||
id = -1;
|
||
maxDepth = -1;
|
||
minDeg = -1.0;
|
||
}
|
||
void info()
|
||
{
|
||
cout << id << " " << maxDepth << " " << minDeg << endl;
|
||
for (int i = 0; i < vert.size(); i++)
|
||
{
|
||
vert.at(i).info();
|
||
}
|
||
}
|
||
void clear_vert()
|
||
{
|
||
if (!vert.empty()) vert.clear();
|
||
}
|
||
};
|
||
typedef vector<polygon> polygonArray;
|
||
|
||
//圆形
|
||
struct circle
|
||
{
|
||
int id;
|
||
int maxDepth;
|
||
double minDeg;
|
||
double radDeg;
|
||
int physic;
|
||
vertex cen;
|
||
};
|
||
typedef vector<circle> circleArray;
|
||
|
||
//三角形信息结构体，包含三角形的三个顶点索引，逆时针排序
|
||
struct triangle
|
||
{
|
||
int ids[3];//三角形顶点
|
||
int physic; //三角形的物理属性组
|
||
triangle() //初始化顶点索引
|
||
{
|
||
physic = 0; //默认的物理属性组为0
|
||
ids[0] = ids[1] = ids[2] = -1;
|
||
}
|
||
};
|
||
typedef vector<triangle> triangleArray;
|
||
|
||
//平面参数结构体
|
||
struct planePara
|
||
{
|
||
double A,B,C,D;
|
||
planePara()
|
||
{
|
||
A = B = C = D = MAX_DBL;
|
||
}
|
||
};
|
||
|
||
//矢量与平面的交点
|
||
cpoint lineOnPlane(cpoint c,cpoint normal,cpoint p)
|
||
{
|
||
cpoint m;
|
||
m.x = 0; m.y = 0; m.z = 0;
|
||
double t;
|
||
if (dot(normal,p) != 0) //平面与矢量平行
|
||
{
|
||
t = dot(normal,c)/dot(normal,p);
|
||
m.x += p.x*t;
|
||
m.y += p.y*t;
|
||
m.z += p.z*t;
|
||
}
|
||
return m;
|
||
}
|
||
|
||
//正二十面体结构体,包含正二十面体的十二个顶点和二十个面的顶点索引,三角面索引按逆时针排序
|
||
struct Icosa
|
||
{
|
||
vertex vert[12];
|
||
triangle tri[20];
|
||
};
|
||
|
||
//四叉树节点结构体,这是整个算法中最重要的结构体，包含一个指向triangle的指针和指向四个子节点的指针
|
||
struct Qdtree_node
|
||
{
|
||
int id;//节点的编号
|
||
int outId;//节点在文件输出时候的id 这个id会在return_leaf函数中确定
|
||
int depth;//节点深度
|
||
bool outOK; //节点是否可以被输出
|
||
triangle* tri;//节点三角形顶点索引 逆时针
|
||
Qdtree_node* children[4];//四个子节点指针
|
||
Qdtree_node() //初始化变量值
|
||
{
|
||
id = -1;
|
||
depth = -1;
|
||
outOK = true;
|
||
children[0] = children[1] = children[2] = children[3] = NULL;
|
||
tri = new triangle;
|
||
}
|
||
void info()
|
||
{
|
||
cout << id << endl;
|
||
cout << depth << endl;
|
||
cout << outOK << endl;
|
||
cout << tri->ids[0] << " " << tri->ids[1] << " " << tri->ids[2] << endl;
|
||
}
|
||
};
|
||
|
||
//四叉树结构
|
||
struct Qdtree
|
||
{
|
||
Qdtree_node *root;//根节点
|
||
};
|
||
|
||
/*
|
||
// 在现行的代码中 我们不再使用平面投影算法 因此不再需要使用以下的数据类型与函数
|
||
struct point2d
|
||
{
|
||
double x,y;
|
||
point2d()
|
||
{
|
||
x = y = MAX_DBL;
|
||
}
|
||
};
|
||
|
||
//由三个顶点计算平面参数
|
||
planePara get_planePara(cpoint v1,cpoint v2,cpoint v3)
|
||
{
|
||
planePara pl;
|
||
pl.A = (v2.y - v1.y)*(v3.z - v1.z) - (v3.y - v1.y)*(v2.z - v1.z);
|
||
pl.B = (v2.z - v1.z)*(v3.x - v1.x) - (v3.z - v1.z)*(v2.x - v1.x);
|
||
pl.C = (v2.x - v1.x)*(v3.y - v1.y) - (v3.x - v1.x)*(v2.y - v1.y);
|
||
pl.D = -1.0*(pl.A*v1.x + pl.B*v1.y + pl.C*v1.z);
|
||
return pl;
|
||
}
|
||
|
||
//点在平面上的投影
|
||
vertex dotOnPlane(planePara pl,vertex v)
|
||
{
|
||
vertex m;
|
||
double t = (pl.A*v.posic.x + pl.B*v.posic.y + pl.C*v.posic.z + pl.D)/(pl.A*pl.A + pl.B*pl.B + pl.C*pl.C);
|
||
m.posic.x = v.posic.x - pl.A*t;
|
||
m.posic.y = v.posic.y - pl.B*t;
|
||
m.posic.z = v.posic.z - pl.C*t;
|
||
m.set(m.posic);
|
||
return m;
|
||
}
|
||
|
||
//以平面内一条直线为x轴 起点为原点 计算另一个点到新的坐标系内的坐标值
|
||
point2d newXY(vertex v1,vertex v2,vertex v3,vertex p)
|
||
{
|
||
point2d p2d;
|
||
vertex m,ap;
|
||
cpoint dir_map,dir_v123;
|
||
m.posic = v2.posic - v1.posic;
|
||
ap.posic = p.posic - v1.posic;
|
||
//因为cos函数在这种情况下自动可以区分正负情况 所以x值的计算比较简单
|
||
p2d.x = dot(ap.posic,m.posic)/length_cpoint(m.posic);
|
||
//下面我们来计算y值 相对比较麻烦 首先计算一下距离
|
||
p2d.y = sqrt(pow(length_cpoint(ap.posic),2) - p2d.x*p2d.x);
|
||
//计算一下三角形的外向法矢量 m与ap的法矢量
|
||
dir_v123 = cross(v2.posic-v1.posic,v3.posic-v1.posic);
|
||
dir_map = cross(v2.posic-v1.posic,ap.posic);
|
||
//如果两个向量同向则y值为正 否则为负
|
||
if (dot(dir_v123,dir_map) < 0) p2d.y *= -1.0;
|
||
return p2d;
|
||
}
|
||
*/
|
||
#endif |