求三维Delaunay四面体剖分算法的各种资料
哪位牛人有关这资料能不能给我一份?
[解决办法]
自己以前随便弄的,你看看把,不是很全了
[code=C/C++][/code]
//多面体体积（凸包）计算（不包括凹多面体）
void CLaserDataEditDlg::OnVolume()
{
// TODO: Add your control notification handler code here
// TODO: Add your control notification handler code here
TriangleList Tlist;//实例一个三角片链
TRIANGLE T;
FPOINT p,p3,p2,pminX,Pn,Pj,Pp,Pq,p_to_p;
//复制一个临时数组，再计算中移除点数据时用
CPointArray p_v_newArray;
p_v_newArray.RemoveAll();
for(int b=0;b<m_pPtArray->GetSize();b++)
{
memcpy(&p_to_p,&m_pPtArray->GetAt(b),sizeof(FPOINT));
p_v_newArray.Add(p_to_p);
}
double DIS = 1000000000;
pminX.x=1000000000;
//获取x坐标最小的坐标,三角形的第一个顶点
for(int i=0;i<p_v_newArray.GetSize();i++)
{
p.x =p_v_newArray.GetAt(i).x;
if(pminX.x >p.x)
{
pminX.x =p.x;
pminX.y =p_v_newArray.GetAt(i).y;
pminX.z =p_v_newArray.GetAt(i).z;
}
else
continue;
}
//获得第二个点，距离最近的点
for(int j=0;j<p_v_newArray.GetSize();j++)
{
Pj =p_v_newArray.GetAt(j);
double Distance =HeightCalculate(&pminX,&Pj);
if(Distance ==0)
{continue;}
if(DIS >Distance)
{
DIS = Distance;
p2 =Pj;
}
if(DIS <Distance)
{
continue;
}
}
//获取第三个顶点
double Dis1 =HeightCalculate(&pminX,&p2);
double cos_value =1;
double cos_temp =1;
for(int n=0;n<p_v_newArray.GetSize();n++)
{
Pn =p_v_newArray.GetAt(n);
if(Pn.x ==pminX.x && Pn.y ==pminX.y && Pn.z ==pminX.z)
{
p_v_newArray.RemoveAt(n);
continue;
}
if(Pn.x ==p2.x && Pn.y ==p2.y && Pn.z ==p2.z)
{
p_v_newArray.RemoveAt(n);
continue;
}
double Dis2 =HeightCalculate(&pminX,&Pn);
double Dis3 =HeightCalculate(&p2,&Pn);
//定义三个角的余弦值，在0-180上是递减函数，值取最小，夹角最大
//余弦定理
cos_temp =(Dis2*Dis2+Dis3*Dis3-Dis1*Dis1)/(2*Dis2*Dis3);
if(cos_value >cos_temp)
{
cos_value =cos_temp;
p3 =Pn;
}
if(cos_value <cos_temp)
{continue;}
}
//移除第三个顶点
for(int m=0;m<p_v_newArray.GetSize();m++)
{
Pp =p_v_newArray.GetAt(m);
if(Pp.x ==p3.x && Pp.y ==p3.y && Pp.z ==p3.z)
{
p_v_newArray.RemoveAt(m);
break;
}
}
//生成第一个三角片
T.line1.PA =pminX;
T.line1.PB =p2;
T.line1.PC =p3;
T.line2.PA =pminX;
T.line2.PB =p3;
T.line2.PC =p2;
T.line3.PA =p2;
T.line3.PB =p3;
T.line3.PC =pminX;
T.mark =1;
//存入三角片链表
Tlist.AddTail(T);
int l=0;
//从其余的点坐标中取出一个点形成第一个三角锥,首先判断是否共面
FPOINT p4,p_centertemp,p_NV;
for (int k=0;k<p_v_newArray.GetSize();k++)
{
p4 =p_v_newArray.GetAt(k);
//判断共面，先要计算三角片的重心，和法向量
p_centertemp =GetCenter(pminX,p2,p3);
p_NV =getNormalVector(pminX,p2,p3);
//判断是否共面
l =Get_The_Four_Point_Pos(p_centertemp,p_NV,p4);
if(l ==0)
{
continue;
}
else
{break;}
}
//四面体新生成三条边，三个三角形,总共六条边，
//新生成的三个三角面{p1,p2,p4}{p1,p3,p4}{p2,p3,p4}
//{p1,p2,p4}
T.line1.PA =pminX;
T.line1.PB =p2;
T.line1.PC =p4;
T.line2.PA =pminX;
T.line2.PB =p4;
T.line2.PC =p2;
T.line3.PA =p2;
T.line3.PB =p4;
T.line3.PC =pminX;
T.mark =2;
Tlist.AddTail(T);
//{p1,p3,p4}
T.line1.PA =pminX;
T.line1.PB =p3;
T.line1.PC =p4;
T.line2.PA =pminX;
T.line2.PB =p4;
T.line2.PC =p3;
T.line3.PA =p3;
T.line3.PB =p4;
T.line3.PC =pminX;
T.mark =3;
Tlist.AddTail(T);
//{p2,p3,p4}
T.line1.PA =p2;
T.line1.PB =p3;
T.line1.PC =p4;
T.line2.PA =p2;
T.line2.PB =p4;
T.line2.PC =p3;
T.line3.PA =p3;
T.line3.PB =p4;
T.line1.PA =p2;
T.mark =4;
Tlist.AddTail(T);
//移除第四个顶点
for(int q=0;q<p_v_newArray.GetSize();q++)
{
Pq =p_v_newArray.GetAt(q);
if(Pq.x ==p3.x && Pq.y ==p3.y && Pq.z ==p3.z)
{
p_v_newArray.RemoveAt(q);
break;
}
}
//实例一个四面体对象与存放四面体的链表
FOUR_TRIANGLE f_t;
F_T_List f_t_list;
POSITION pos1,pos2,pos3,pos4;
pos1 =Tlist.FindIndex(0);
pos2 =Tlist.FindIndex(1);
pos3 =Tlist.FindIndex(2);
pos4 =Tlist.FindIndex(3);
f_t.Ta =Tlist.GetAt(pos1);
f_t.Tb =Tlist.GetAt(pos2);
f_t.Tc =Tlist.GetAt(pos3);
f_t.Td =Tlist.GetAt(pos4);
f_t.T_mark =0;
f_t_list.AddTail(f_t);
//扩展三角锥
FOUR_TRIANGLE f_t_new;
TRIANGLE tn,T_new1,T_new2,T_new3;
FPOINT p_n,p_m,p_1,p_2,p_3,p_center,p_f;
POSITION pos_f_t;
//（循环）由三角锥链表里面取出一个三角锥
for(int u =0;u<f_t_list.GetCount();u++)
{
pos_f_t =f_t_list.FindIndex(u);
f_t =f_t_list.GetAt(pos_f_t);
//（循环）由取出的三角锥中，取出一个面，和面所对的顶点
for(int h =0;h<4;h++)
{
if(h ==0)
{
tn =f_t.Ta;
p_n =GetPointNotInTriangle(tn,f_t);
}
if(h ==1)
{
tn =f_t.Tb;
p_n =GetPointNotInTriangle(tn,f_t);
}
if(h ==2)
{
tn =f_t.Tc;
p_n =GetPointNotInTriangle(tn,f_t);
}
if(h ==3)
{
tn =f_t.Td;
p_n =GetPointNotInTriangle(tn,f_t);
}
//取出面与面所对顶点后，先计算面的中心和法向，计算面的顶点与法向的方向的正负，0为共面，1为同向，2为反向
p_1 =tn.line1.PA;
p_2 =tn.line1.PB;
p_3 =tn.line1.PC;
p_center =GetCenter(p_1,p_2,p_3);
p_f =getNormalVector(p_1,p_2,p_3);
int p_n_d =Get_The_Four_Point_Pos(p_center,p_f,p_n);
//（循环）由数组中剩余数据取出点,计算是否和法相同向
for(int g =0;g<p_v_newArray.GetSize();g++)
{
p_m =p_v_newArray.GetAt(g);
int p_m_d =Get_The_Four_Point_Pos(p_center,p_f,p_m);
//判断是否和面所对的顶点在同一方向
if(p_m_d ==0)//和三角面共面
{ continue; }
if(p_m_d *p_n_d >0)//同向
{ continue; }
if(p_m_d *p_n_d <0)//反向可以存形成新的三角锥，存三角锥链
{
//新三面{p_m,p_1,p_2}{p_m,p_1,p_3}{p_m,p_2,p_3}
T_new1.line1.PA =p_m;
T_new1.line1.PB =p_1;
T.line1.PC =p_2;
T_new1.line2.PA =p_m;
T_new1.line2.PB =p_2;
T_new1.line2.PC =p_1;
T_new1.line3.PA =p_1;
T_new1.line3.PB =p_2;
T_new1.line3.PC =p_m;
Tlist.AddTail(T);
//{p_m,p_1,p_3}
T_new2.line1.PA =p_m;
T_new2.line1.PB =p_1;
T_new2.line1.PC =p_3;
T_new2.line2.PA =p_m;
T_new2.line2.PB =p_3;
T_new2.line2.PC =p_1;
T_new2.line3.PA =p_1;
T_new2.line3.PB =p_3;
T_new2.line3.PC =p_m;
Tlist.AddTail(T);
//{p_m,p_2,p_3}
T_new3.line1.PA =p_m;
T_new3.line1.PB =p_2;
T_new3.line1.PC =p_3;
T_new3.line2.PA =p_m;
T_new3.line2.PB =p_3;
T_new3.line2.PC =p_2;
T_new3.line3.PA =p_3;
T_new3.line3.PB =p_2;
T_new3.line1.PA =p_m;
Tlist.AddTail(T);
f_t_new.Ta =tn;
f_t_new.Tb =T_new1;
f_t_new.Tc =T_new2;
f_t_new.Td =T_new3;
f_t_list.AddTail(f_t);
//添加完后要再离散点数组中移除 刚使用过的 形成新的顶点的 坐标数据
for(int q=0;q<p_v_newArray.GetSize();q++)
{
Pq =p_v_newArray.GetAt(q);
if(Pq.x ==p_m.x && Pq.y ==p_m.y && Pq.z ==p_m.z)
{
p_v_newArray.RemoveAt(q);
break;
}
}
}
}
}
}
//计算体积，由三角锥中取出所生成的三角锥
FOUR_TRIANGLE ft;
TRIANGLE T_v;
FPOINT pa,pb,pc,pd;
POSITION pos_v;
double v =0, V =0;
for(int a =0;a<f_t_list.GetCount();a++)
{
pos_v =f_t_list.FindIndex(a);
ft =f_t_list.GetAt(pos_v);
T_v =ft.Ta;
pa =T_v.line1.PA;
pb =T_v.line1.PB;
pc =T_v.line1.PC;
pd =GetPointNotInTriangle(T_v,ft);
//
//四面体体积
//
v =VOLUME(pa,pb,pc,pd);
V =V +v;
}
V =fabs(V);
CString str;
str.Format("体积是：%f",V);
SetDlgItemText(IDC_RESULT_SHOW,str);
}
