?? 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 + -