📄 area0.cpp
字号:
/*4月29日
算法：在八个子区域内分别找到最小的四面体，用其外接球半径进行探索
运行结果显示这种算法不对
算出的探索球内包含所有的点
*/
#include"iostream.h"
#include"iomanip.h"
#include"math.h"
#include"fstream.h"
#include"stdlib.h"
double det(double a[3][3]);
const double pi=3.1415926;
double test_r=5;
/**********************************************
***********************************************
类定义：点，球，四面体
***********************************************
**********************************************/
/*****************************************
点的定义，计算两点之间的距离
****************************************/
class point
{
public:
double x,y,z; //坐标
point(double x=0,double y=0,double z=0); //构造函数 默认值为(0,0,0)
double distance(point); //求两点间的距离
point operator -(point p1) ; //重载“-”号，使其为两点坐标对应相减
void operator =(point p); //重载“=”号，使其实现类的赋值
void dis();
};
inline point::point(double a1,double a2,double a3)
//构造函数
{
x=a1;y=a2;z=a3;
}
point point::operator -(point p)
//重载“-”号，使其为两点坐标对应相减
{
return point(x-p.x,y-p.y,z-p.z);
}
void point::operator =(point p)
//重载“=”号，使其实现类的赋值
{
x=p.x;y=p.y;z=p.z;
}
double point::distance(point p2)
//计算两点之间的距离
{
double dis;
dis=sqrt((x-p2.x)*(x-p2.x)+(y-p2.y)*(y-p2.y)+(z-p2.z)*(z-p2.z));
return dis;
}
void point::dis()
//显示点的坐标
{
cout<<"此点的坐标=("<<x<<","<<y<<","<<z<<")"<<endl;
}
point pnt[10000];
int vexcnt;
/*********************************************
球的定义，判断点是否在球内
**********************************************/
class sphere
{
private:
double r; //半径
point centre; //球心坐标
public:
sphere(point a,double rad); //初始化圆心和半径
int judge_sphere(point);
};
sphere::sphere(point a,double rad):centre(a),r(rad)
//初始化圆心和半径
{
}
int sphere::judge_sphere(point p)
//判断点是否在球内
//在球内返回 1，否则返回 0
{
double point::distance(point);
if(centre.distance(p)<r)
return 1;
else return 0;
}
/*****************************************
四面体的定义，内切球，外界球
*******************************************/
class tetrahedron
{
private:
int pn1,pn2,pn3,pn4; //四个节点的编号
//内切球和外接球
double quanlity; //单元的质量
public:
tetrahedron(int pn1=0,int pn2=0,int pn3=0,int pn4=0); //初始化四面体的顶点的编号
double r_in_sph(); //内切球的半径
double r_out_sph(); //外接球的半径
point center_out(); //外接球的球心
point center_in(); //内切球的球心
void reset(int,int,int,int); //重置四面体四个节点的编号
};
tetrahedron::tetrahedron(int p1,int p2,int p3,int p4)
//初始化四面体的顶点的编号
{
pn1=p1;pn2=p2;pn3=p3;pn4=p4;
}
point tetrahedron::center_out()
//球心到任两个点距离相等
{
int i,j;
point pp(0,0,0); //pp返回球心
double a[3][3],b[3],d[3][3];
pp=pnt[pn1]-pnt[pn4];
a[0][0]=pp.x;a[0][1]=pp.y;a[0][2]=pp.z;
pp=pnt[pn2]-pnt[pn3];
a[1][0]=pp.x;a[1][1]=pp.y;a[1][2]=pp.z;
pp=pnt[pn4]-pnt[pn3];
a[2][0]=pp.x;a[2][1]=pp.y;a[2][2]=pp.z;
b[0]=(pnt[pn1].x*pnt[pn1].x+pnt[pn1].y*pnt[pn1].y+pnt[pn1].z*pnt[pn1].z
-pnt[pn4].x*pnt[pn4].x-pnt[pn4].y*pnt[pn4].y-pnt[pn4].z*pnt[pn4].z)/2;
b[1]=(pnt[pn2].x*pnt[pn2].x+pnt[pn2].y*pnt[pn2].y+pnt[pn2].z*pnt[pn2].z
-pnt[pn3].x*pnt[pn3].x-pnt[pn3].y*pnt[pn3].y-pnt[pn3].z*pnt[pn3].z)/2;
b[2]=(pnt[pn4].x*pnt[pn4].x+pnt[pn4].y*pnt[pn4].y+pnt[pn4].z*pnt[pn4].z
-pnt[pn3].x*pnt[pn3].x-pnt[pn3].y*pnt[pn3].y-pnt[pn3].z*pnt[pn3].z)/2;
for(i=0;i<3;i++)
{for(j=1;j<3;j++)
d[i][j]=a[i][j];
d[i][0]=b[i];
}
pp.x=det(d)/det(a);
for(i=0;i<3;i++)
{for(j=0;j<3&&j!=2;j++)
d[i][j]=a[i][j];
d[i][1]=b[i];
}
pp.y=det(d)/det(a);
for(i=0;i<3;i++)
{for(j=0;j<2;j++)
d[i][j]=a[i][j];
d[i][2]=b[i];
}
pp.z=det(d)/det(a);
return pp;
}
double tetrahedron::r_out_sph()
//求四面体外接球的半径
{
double a,b,c,v;
a=pnt[pn1].distance(pnt[pn3])*pnt[pn2].distance(pnt[pn4]);
b=pnt[pn1].distance(pnt[pn4])*pnt[pn2].distance(pnt[pn3]);
c=pnt[pn1].distance(pnt[pn2])*pnt[pn3].distance(pnt[pn4]);
point pp(0,0,0);
double d[3][3];
pp=pnt[pn2]-pnt[pn1];
d[0][0]=pp.x;d[1][0]=pp.y;d[2][0]=pp.z;
pp=pnt[pn3]-pnt[pn1];
d[0][1]=pp.x;d[1][1]=pp.y;d[2][1]=pp.z;
pp=pnt[pn4]-pnt[pn1];
d[0][2]=pp.x;d[1][2]=pp.y;d[2][2]=pp.z;
v=det(d);
v=fabs(v)/6.0;
return sqrt((a+b+c)*(a+b-c)*(b+c-a)*(a-b+c))/(24*v);
}
void tetrahedron::reset(int n1,int n2,int n3,int n4)
//重置四面体四个节点的编号
{
pn1=n1;pn2=n2;pn3=n3;pn4=n4;
}
/*************************************************
**************************************************
管理函数：求最大值
**************************************************
*************************************************/
double doble_max(double a[],int n)
//求数组a[n]的最大值
{
int i;
double max=a[0];
for(i=1;i<n;i++)
if(max<a[i])max=a[i];
return max;
}
double det(double a[3][3])
//求3阶方阵a的行列式
{
return
a[0][0]*a[1][1]*a[2][2]+a[0][1]*a[1][2]*a[2][0]+a[0][2]*a[2][1]*a[1][0]
-a[2][0]*a[1][1]*a[0][2]-a[1][0]*a[0][1]*a[2][2]-a[0][0]*a[1][2]*a[2][1];
}
int judge_conplan(point p1,point p2,point p3,point p4)
//判断空间四点是否共面 共面为1 不共面为0
{
point pp;
double d[3][3];
pp=p2-p1;
d[0][0]=pp.x;d[1][0]=pp.y;d[2][0]=pp.z;
pp=p3-p1;
d[0][1]=pp.x;d[1][1]=pp.y;d[2][1]=pp.z;
pp=p4-p1;
d[0][2]=pp.x;d[1][2]=pp.y;d[2][2]=pp.z;
if(det(d)<1.e-6)
return 1;
else return 0;
}
point chaji(point v1,point v2)
//计算空间两个向量的外积
{
point p;
p.x=v1.y*v2.z-v2.y*v1.z;
p.y=-(v1.x*v2.z-v2.x*v1.z);
p.z=v1.x*v2.y-v2.x*v1.y;
return p;
}
int judge_shareline(point p1,point p2,point p3)
//判断空间三点是否共线 共线为1 不共线为0
{
point pp;
pp=chaji(p3-p1,p2-p1);
if(pp.x<1.e-6&&pp.y<1.e-6&&pp.z<1.e-6)
return 1;
else return 0;
}
double optimize_rad(point a,int a_index,double test_r)
//计算以a为球心的优化探索求半径
{
double opr[8]={0.0}; //记录各个子区域的优化半径
int i=0,j,count[8]={0}; //记录各个象限的节点数
int vex_index[8][10000]; //各行记录一个象限的点的编号
double dis[10000];
sphere explore_sphere(a,test_r);
//for(i=0;i<27;i++)
while(i<vexcnt&&explore_sphere.judge_sphere(pnt[i]))
{
if(pnt[i].x>=a.x&&pnt[i].y>a.y&&pnt[i].z>a.z)
{vex_index[0][count[0]]=i;count[0]++;} //第一象限的点+yoz
else if(pnt[i].x<a.x&&pnt[i].y>=a.y&&pnt[i].z>a.z)
{vex_index[1][count[1]]=i;count[1]++;} //第二象限的点+xoz
else if(pnt[i].x<=a.x&&pnt[i].y<a.y&&pnt[i].z>a.z)
{vex_index[2][count[2]]=i;count[2]++;} //第三象限的点+yoz
else if(pnt[i].x>a.x&&pnt[i].y<=a.y&&pnt[i].z>=a.z)
{vex_index[3][count[3]]=i;count[3]++;} //第四象限的点+xoy+xoz
else if(pnt[i].x>=a.x&&pnt[i].y>a.y&&pnt[i].z<=a.z)
{vex_index[4][count[4]]=i;count[4]++;} //第五象限的点+xoy+yoz
else if(pnt[i].x<a.x&&pnt[i].y>=a.y&&pnt[i].z<=a.z)
{vex_index[5][count[5]]=i;count[5]++;} //第六象限的点+xoy+xoz
else if(pnt[i].x<=a.x&&pnt[i].y<a.y&&pnt[i].z<=a.z)
{vex_index[6][count[6]]=i;count[6]++;} //第七象限的点+xoy+yoz
else if(pnt[i].x!=a.x||pnt[i].y!=a.y||pnt[i].z!=a.z)
{vex_index[7][count[7]]=i;count[7]++;} //第八象限的点+xoz
i++;
} //while{}
for(i=0;i<8;i++)
{for(j=0;j<count[i];j++)
dis[j]=a.distance(pnt[vex_index[i][j]]); //计算子区域内各点距球心a的距离
int k,swap;
double temp;
for(j=0;j<count[i]-1;j++)
for(k=j+1;k<count[i];k++)
if(dis[j]>dis[k])
{temp=dis[j];dis[j]=dis[k];dis[k]=temp;
swap=vex_index[i][j];vex_index[i][j]=vex_index[i][k];vex_index[i][k]=swap;
} //冒泡排序，将子区域内各点的序号按与球心距离从小到大排序
if(count[i]<3) //子区域内节点数小于三个，不能组成四面体
{ opr[i]=dis[count[i]-1];
continue;
}
if(judge_conplan(pnt[a_index],pnt[vex_index[i][0]],pnt[vex_index[i][1]],pnt[vex_index[i][2]])) //如果四点共面
{opr[i]=0;
continue;
}
tetrahedron tet0(a_index,vex_index[i][0],vex_index[i][1],vex_index[i][2]); //距球心距离最近的三个点与其组成一个四面体
sphere out_sph(tet0.center_out(),tet0.r_out_sph()); //四面体的外接球
for(j=3;j<count[i];j++)
if(out_sph.judge_sphere(pnt[vex_index[i][j]])) //有其他点在此四面体的外接球内
opr[i]=dis[j];
if(j==count[i])
opr[i]=2*tet0.r_out_sph(); //不包含其他点时为外接球直径
}
/*for(i=0;i<8;i++) //外层for(i=0;i<8;i++)
cout<<opr[i]<<endl;*/
return doble_max(opr,8); //返回优化探索球半径
} //optimize_rad函数结束
int cnt_in_optsph(int i,int a[])
//将编号为i的点的优化探索球域内包含的点的编号存入数组a中，返回数组的大小
{
int j;
int cnt=0;
sphere sph(pnt[i],optimize_rad(pnt[i],i,test_r));
for(j=0;j<vexcnt;j++)
if(sph.judge_sphere(pnt[j]))
{a[cnt]=j;cnt++;}
return cnt;
}
/*********************************************************
输入 输出函数
*********************************************************/
void read()
//读取节点坐标和节点个数
{
int i;
char ch;
fstream vertex("vertex.txt",ios::in|ios::out);
vexcnt=0;
if(!vertex)
{
cerr<<"vertex.txt can't open"<<endl;
exit(0);
}
while(!vertex.eof())
{
vertex>>pnt[vexcnt].x>>pnt[vexcnt].y>>pnt[vexcnt].z;
vexcnt++;
}
// while(vertex.read((char*)&i,sizeof(double)))
}
void sta_r_pcnt()
//统计节点的探索球半径大小和包含的节点数目 并输出
{
ofstream sta_r("sta_r.txt",ios::out);
ofstream sta_pcnt("sta_pcnt.txt",ios::out);
/* if(!sta)
{cerr<<"文件打开失败"<<endl;
exit(0);
}
*/
int r[3]={0},pcnt[4]={0};
int a[10000];
double opr[10000];
for(int i=0;i<vexcnt;i++)
opr[i]=optimize_rad(pnt[i],i,test_r);
double maxr=doble_max(opr,vexcnt);
for(i=0;i<vexcnt;i++)
{ double or=opr[i]/maxr; //将半径归1化
if(or<=0.05)
r[0]++;
else if(or<=0.1)
r[1]++;
else
r[2]++;
int pnum=cnt_in_optsph(i,a);
if(pnum<=20)
pcnt[0]++;
else if(pnum<=100)
pcnt[1]++;
else if(pnum<=200)
pcnt[2]++;
else
pcnt[3]++;
}
sta_r<<" r的区间"<<setw(15)<<"节点数"<<setw(15)<<"百分比"<<endl;
sta_r<<"(0，0.05]"<<setw(15)<<r[0]<<setw(15)<<double(r[0]*100/vexcnt)<<endl;
sta_r<<"(0.05，0.1]"<<setw(15)<<r[1]<<setw(15)<<double(r[1]*100/vexcnt)<<endl;
sta_r<<"(0.1，0.5]"<<setw(15)<<r[2]<<setw(15)<<double(r[2]*100/vexcnt)<<endl;
sta_pcnt<<"节点数区间"<<setw(15)<<"节点数"<<setw(15)<<"百分比"<<endl;
sta_pcnt<<"(0,20]"<<setw(15)<<pcnt[0]<<setw(15)<<double(pcnt[0]*100/vexcnt)<<endl;
sta_pcnt<<"(20,50]"<<setw(15)<<pcnt[1]<<setw(15)<<double(pcnt[1]*100/vexcnt)<<endl;
sta_pcnt<<"(50,100]"<<setw(15)<<pcnt[2]<<setw(15)<<double(pcnt[2]*100/vexcnt)<<endl;
sta_pcnt<<"(100,inf)"<<setw(15)<<pcnt[3]<<setw(15)<<double(pcnt[3]*100/vexcnt)<<endl;
}
/*************************************************
main函数
**************************************************/
void main()
{
read();
sta_r_pcnt();
}⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -