Tin三角网的建立

#include "StdAfx.h"
#include "Mstruct.h"
#include "math.h"
CMstruct::CMstruct(void)
{
}

CMstruct::~CMstruct(void)
{

}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// LoadFile()
// Parameters:
//
// Returns:
//
// description:
// 加载文件和数据的存储
// -----------------------------------------------------------------------------------------------------------------------
bool CMstruct::LoadFile()
{
//// 清空缓存数
m_Ps.clear();
m_Ts.clear();
m_Effective.clear();
m_Line.clear();

CString filename;
filename = Open();
if (filename != _T(""))
{
CStdioFile file;
if (file.Open(filename,CFile::modeRead))
{
CString str;
int iddata = 0;
while (file.ReadString(str))
{
Ps psdata;
CString array[5];
int n = 0;
int i = 0;

// 遇到分隔符时，返回以一个分隔符和这个之间的字符
for (;n != -1;i++)
{
array[i] = str.Tokenize(_T(","),n);
}
psdata.X = _wtof(array[0]);
psdata.Y = _wtof(array[1]);
if (array[2] == _T(""))
{
array[2] = _T("0.0");
}
psdata.Z = _wtof(array[2]);
psdata.ID = iddata;
m_Ps.push_back(psdata);
iddata++;
}
}
return true;
}
else
{
return false;
}
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// Open()
// Parameters:
//
// Returns:
//
// description:
// 通过通用对话框选择文件
// -----------------------------------------------------------------------------------------------------------------------
CString CMstruct::Open()
{
CString filename;
CString fileext;//文件扩展名
CString filepathname=_T("");
//下边为弹出对话框，选择要读取的文件
CFileDialog fpdlg(TRUE, NULL, NULL,OFN_HIDEREADONLY | OFN_OVERWRITEPROMPT,_T("(文件类型)(*.txt)|*.txt|(文件类型)(*.txt)|*.txt|所有文件(*.*)|*.*||"),NULL);
if(fpdlg.DoModal() == IDOK)
{
filename = fpdlg.GetFileName();
fileext = fpdlg.GetFileExt();//文件扩展名
filepathname = fpdlg.GetPathName();
return filepathname;
}
return filepathname;
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// ClosePoint(Ps psdata)
// Parameters:
// Ps psdata 点的数据对象
// Returns:
// id 返回最近点的id
// description:
// 找到该点的最近点
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::ClosePoint(Ps psdata)
{
float d,dmin=0;
int idmin;
int m = 0;
vector::iterator it = m_Ps.begin();
for (;it!=m_Ps.end();it++)
{

// 当遇到本身的点后，直接跳过

去执行下一个循环
if ((it->X == psdata.X) && (it->Y == psdata.Y)&&(it->Z == psdata.Z))
{
continue;
}
d = sqrt(pow((it->X-psdata.X),2)+pow((it->Y-psdata.Y),2)+pow((it->Z-psdata.Z),2));
if (m == 0)
{
idmin = it->ID;
dmin = d;
m = 1;
}
if (d{
dmin = d;
idmin = it->ID;
}
}
return idmin; // 返回最近点的id
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// FirstClosePoint
// Parameters:
// Ps psdata
// Returns:
// 最近点的ID
// description:
// 当插入点不在任何三角形内时找离他最近的数据点
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::FirstClosePoint(Ps psdata)
{
float d,dmin=0;
int idmin;
int m = 0;
vector::iterator it = m_Ps.begin();
for (;it!=m_Ps.end();it++)
{

// 当遇到本身的点后，直接跳过去执行下一个循环
if ((it->X == psdata.X) && (it->Y == psdata.Y)&&(it->Z == psdata.Z))
{
continue;
}
d = sqrt(pow((it->X-psdata.X),2)+pow((it->Y-psdata.Y),2)/*+pow((it->Z-psdata.Z),2)*/);
if (m == 0)
{
idmin = it->ID;
dmin = d;
m = 1;
}
if (d{
dmin = d;
idmin = it->ID;
}
}
return idmin; // 返回最近点的id
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// CreateTin
// Parameters:
//
// Returns:
//
// description:
// 创建三维网格
// -----------------------------------------------------------------------------------------------------------------------
void CMstruct::CreateTin()
{
int L = 1;
int K;
Ps psdata;
Ts tsdata;
tsdata.IB1 = m_Ps.front().ID;
tsdata.IB2 = ClosePoint(m_Ps.front()); //找到最近点并返回id
psdata = FindPt(tsdata.IB2);

tsdata.IB3 = EntClosePt(m_Ps.front(),psdata);
m_Ts.push_back(tsdata);
for (K=0;K!=L;K++)/*K为种子三角形的指针，L为当前三角形的指针，通过递归生成三角网，在K=L时结束。*/
{
Ts dts;
int i = 0;
vector ::iterator it = m_Ts.begin();
for (;it!=m_Ts.end();it++)
{
if (K == i)
{
dts.IB1 = it->IB1;
dts.IB2 = it->IB2;
dts.IB3 = it->IB3;
}
i++;
}
if (m_Ts.size() == 1)
{
IsEffective(dts,1);
if (m_Effective.size()>0) /*在1-2边新生成三角形，排除没有有效点的情况。*/
{
Ts newts; // 新的三角形数据
newts.IB1 = dts.IB1;
newts.IB2 = dts.IB2;
Ps pt1,pt2;
pt1 = FindPt(dts.IB1);
pt2 = FindPt(dts.IB2);
newts.IB3 = EntClosePt(pt1,pt2);
if (IsRepeat(newts))
{
m_Ts.push_back(newts);
L++;
}
}
}
IsEffective(dts,2);
if (m_Effective.size()>0) /*在2-3边新生成三角形，排除没有有效

点的情况。*/
{
Ts newts; // 新的三角形数据
newts.IB1 = dts.IB2;
newts.IB2 = dts.IB3;
Ps pt1,pt2;
pt1 = FindPt(dts.IB2);
pt2 = FindPt(dts.IB3);
newts.IB3 = EntClosePt(pt1,pt2);
if (IsRepeat(newts))
{
m_Ts.push_back(newts);
L++;
}
}
IsEffective(dts,3);
if (m_Effective.size()>0) /*在1-3边新生成三角形，排除没有有效点的情况。*/
{
Ts newts; // 新的三角形数据
newts.IB1 = dts.IB1;
newts.IB2 = dts.IB3;
Ps pt1,pt2;
pt1 = FindPt(dts.IB1);
pt2 = FindPt(dts.IB3);
newts.IB3 = EntClosePt(pt1,pt2);
if (IsRepeat(newts))
{
m_Ts.push_back(newts);
L++;
}
}
}
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// EntClosePt(Ps pt1,Ps pt2)
// Parameters:
// Ps pt1,Ps pt2 一边的起点和端点
// Returns:
// id 返回最近点的id
// description:
// 找到三角形一边的最近点
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::EntClosePt(Ps pt1,Ps pt2)
{
int idmax;
float d1,d2,d3,s,C,Cmax=0;

// 判断是否已经构建出三角形
if (m_Ts.empty())
{
vector ::iterator it = m_Ps.begin();
for (;it!=m_Ps.end();it++)
{
if(((it->X == pt1.X) && (it->Y == pt1.Y)&&(it->Z == pt1.Z)) || ((it->X == pt2.X) && (it->Y == pt2.Y)&&(it->Z == pt2.Z)))
{
continue;
}
/*求三点连接的三边的长度*/
d1=sqrt(pow((it->X-pt1.X),2)+pow((it->Y-pt1.Y),2)+pow((it->Z-pt1.Z),2));
d2=sqrt(pow((it->X-pt2.X),2)+pow((it->Y-pt2.Y),2)+pow((it->Z-pt2.Z),2));
d3=sqrt(pow((pt1.X-pt2.X),2)+pow((pt1.Y-pt2.Y),2)+pow((pt1.Z-pt2.Z),2));
s=(d1*d1+d2*d2-d3*d3)/(2*d1*d2);
if(s<1)C=acos(s);else C=0; /*用arccos求角度*/
if(C>Cmax) /*找最大角*/
{
Cmax=C;
idmax = it->ID;
}
}
return idmax;
}
vector ::iterator it = m_Effective.begin();
for (;it!=m_Effective.end();it++)
{
if(((it->X == pt1.X) && (it->Y == pt1.Y)&&(it->Z == pt1.Z)) || ((it->X == pt2.X) && (it->Y == pt2.Y)&&(it->Z == pt1.Z)))
{
continue;
}
/*求三点连接的三边的长度*/
d1=sqrt(pow((it->X-pt1.X),2)+pow((it->Y-pt1.Y),2)+pow((it->Z-pt1.Z),2));
d2=sqrt(pow((it->X-pt2.X),2)+pow((it->Y-pt2.Y),2)+pow((it->Z-pt2.Z),2));
d3=sqrt(pow((pt1.X-pt2.X),2)+pow((pt1.Y-pt2.Y),2)+pow((pt1.Z-pt2.Z),2));
s=(d1*d1+d2*d2-d3*d3)/(2*d1*d2);
if(s<1)C=acos(s);else C=0; /*用arccos求角度*/
if(C>Cmax) /*找最大角*/
{
Cmax=C;
idmax = it->ID;
}
}
return idmax;
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// Draw()
// Parameters:
//
// Returns:
//
// description:
// 在屏幕上画出三角形网格
//

-----------------------------------------------------------------------------------------------------------------------
void CMstruct::Draw()
{
CreateTin();
vector::iterator it = m_Ts.begin();
for (;it!=m_Ts.end();it++)
{
Ps pt1,pt2,pt3;
bool f;

pt1 = FindPt(it->IB1);
pt2 = FindPt(it->IB2);
pt3 =FindPt(it->IB3);
AcGePoint3d point1(pt1.X,pt1.Y,pt1.Z);
AcGePoint3d point2(pt2.X,pt2.Y,pt2.Z);
AcGePoint3d point3(pt3.X,pt3.Y,pt3.Z);
AcDbLine *pline1 =new AcDbLine;
AcDbLine *pline2 =new AcDbLine;
AcDbLine *pline3 =new AcDbLine;
pline1->setStartPoint(point1);
pline1->setEndPoint(point2);
f = IsDrawed(it->IB1,it->IB2);
if (f)
{
line lidata;
lidata.b1 = it->IB1;
lidata.b2 = it->IB2;
m_Line.push_back(lidata);
HP::AppendAndClose(pline1);

}
else
{
delete pline1;
pline1 = NULL;
}


pline2->setStartPoint(point2);
pline2->setEndPoint(point3);

f = IsDrawed(it->IB2,it->IB3);
if (f)
{
line lidata;
lidata.b1 = it->IB2;
lidata.b2 = it->IB3;
m_Line.push_back(lidata);
HP::AppendAndClose(pline2);

}
else
{
delete pline2;
pline2 = NULL;
}


pline3->setStartPoint(point1);
pline3->setEndPoint(point3);
f = IsDrawed(it->IB1,it->IB3);
if (f)
{
line lidata;
lidata.b1 = it->IB1;
lidata.b2 = it->IB3;
m_Line.push_back(lidata);
HP::AppendAndClose(pline3);

}
else
{
delete pline3;
pline3 = NULL;
}

}

}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// IsEffective()
// Parameters:
// int n 选择条件
// Returns:
//
// description:
// 找出各边有效点
// -----------------------------------------------------------------------------------------------------------------------
void CMstruct::IsEffective(Ts tsdata,int n)
{
m_Effective.clear();// 清空数据
int p1,p2,p3;
int m1,m2;
switch(n)
{
case 1:
{
p1=tsdata.IB1; /*以1-2边为界*/
p2=tsdata.IB2;
p3=tsdata.IB3;
break;
}
case 2:
{
p1=tsdata.IB2; /*以2-3边为界*/
p2=tsdata.IB3;
p3=tsdata.IB1;
break;
}
case 3:
{
p1=tsdata.IB1; /*以1-3边为界*/
p2=tsdata.IB3;
p3=tsdata.IB2;
}
}
Ps pt1,pt2,pt3;
pt1 = FindPt(p1);
pt2 = FindPt(p2);
pt3 = FindPt(p3);

if(pt3.Y-(pt2.Y-pt1.Y)/(pt2.X-pt1.X)*pt3.X-(pt1.Y*pt2.X-pt2.Y*pt1.X)/(pt2.X-pt1.X)<0)
{
m1 = -1;
}
else
{
m1 = 1; /*求第3点所在区间，正为1，负为-1*/
}
vector ::iterator it = m_Ps.begin();
for (;it!=m_Ps.end();it++)
{
if(((it->X == pt1.X) && (it->Y == pt1.Y)) || ((it->X == pt2.X) && (it->Y == pt2.Y)))
{
continue;
}
if (it->Y-(pt2.Y-pt1.Y)/(pt2.X-pt1.X)*it->X-(pt1.Y*pt2.X-pt2.Y*pt1.X)/(pt2.X-pt1.X)<0)
{
m2 = -1;
}
else
{
m2 = 1;
}
if (m1 != m2)
{
Ps psda

ta;
psdata.ID = it->ID;
psdata.X = it->X;
psdata.Y = it->Y;
m_Effective.push_back(psdata);
}
}
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// IsRepeat()
// Parameters:
// Ts tsdata 三角形的数据对象
// Returns:
//
// description:
// 判断三角形是否重复
// -----------------------------------------------------------------------------------------------------------------------
bool CMstruct::IsRepeat(Ts tsdata)
{
vector ::iterator it = m_Ts.begin();
for (;it!=m_Ts.end();it++)
{
/*重复的三角形三个点雷同有6种排列组合*/
if((it->IB1==tsdata.IB1&&it->IB2==tsdata.IB2&&it->IB3==tsdata.IB3)||
(it->IB1==tsdata.IB1&&it->IB2==tsdata.IB3&&it->IB3==tsdata.IB2)||
(it->IB1==tsdata.IB2&&it->IB2==tsdata.IB1&&it->IB3==tsdata.IB3)||
(it->IB1==tsdata.IB2&&it->IB2==tsdata.IB3&&it->IB3==tsdata.IB1)||
(it->IB1==tsdata.IB3&&it->IB2==tsdata.IB1&&it->IB3==tsdata.IB2)||
(it->IB1==tsdata.IB3&&it->IB2==tsdata.IB2&&it->IB3==tsdata.IB1))
return false; /*重复返回false*/
}
return true;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// FindPt(int id)
// Parameters:
// int id 点的id
// Returns:
// 返回点的数据
// description:
// 通过id找到该点的具体数据
// -----------------------------------------------------------------------------------------------------------------------
Ps CMstruct::FindPt(int id)
{
Ps psdata;
vector::iterator it = m_Ps.begin();
for (;it!=m_Ps.end();it++)
{
if (it->ID == id)
{
psdata.ID = it->ID;
psdata.X = it->X;
psdata.Y = it->Y;
psdata.Z = it->Z;
break;
}
}
return psdata;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// IsDrawed
// Parameters:
// int id1,int id2 两点的id
// Returns:
//
// description:
// 判断经过这两点的直线是否已经画出
// -----------------------------------------------------------------------------------------------------------------------
bool CMstruct::IsDrawed(int id1,int id2)
{
if (m_Line.size() == 0)
{
return true;
}
vector::iterator it = m_Line.begin();
for (;it!=m_Line.end();it++)
{
// 已经画出直线
if ((id1==it->b1 && id2==it->b2)||(id2==it->b1 && id1==it->b2))
{
return false;
}
}
return true;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// PointZ()
// Parameters:
//
// Returns:
// 返回一个点的z坐标
// description:
// 计算出高程值
// ------------------------------------------------------------------------------

-----------------------------------------
float CMstruct::PointZ(double px,double py)
{
float A,B,C,D;
float z;
A = Ptdata[0].Y*(Ptdata[1].Z-Ptdata[2].Z)+Ptdata[1].Y*(Ptdata[2].Z-Ptdata[0].Z)+Ptdata[2].Y*(Ptdata[0].Z-Ptdata[1].Z);
B = Ptdata[0].Z*(Ptdata[1].X-Ptdata[2].X)+Ptdata[1].Z*(Ptdata[2].X-Ptdata[0].X)+Ptdata[2].Z*(Ptdata[0].X-Ptdata[1].X);
C = Ptdata[0].X*(Ptdata[1].Y-Ptdata[2].Y)+Ptdata[1].X*(Ptdata[2].Y-Ptdata[0].Y)+Ptdata[2].X*(Ptdata[0].Y-Ptdata[1].Y);
D = -(A*Ptdata[0].Y+B*Ptdata[0].Y+C*Ptdata[0].Z);
z = -(A*px+B*py+D)/C;
return z;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// SearchTriangle
// Parameters:
// pt 插入点
// Returns:
// 0 or 1 若找到三角形 返回0 否则返回1
// description:
// 该函数是为了求出插入点所在的三角形的平面的三个顶点
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::SearchTriangle(AcGePoint3d pt)
{
int i;
i = FindCloseLine(pt);
vector::iterator it = m_Ts.begin();
for (;it!=m_Ts.end();it++)
{
Ps p1,p2,p3;
p1 = FindPt(it->IB1);
p2 = FindPt(it->IB2);
p3 = FindPt(it->IB3);
// 判断边所在的三角形，有六种情况
if (m_Line[i].b1 == it->IB1&&m_Line[i].b2 == it->IB2)
{
int f = IsIntersect(it->IB1,it->IB2,it->IB3,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}
if (m_Line[i].b1 == it->IB1&&m_Line[i].b2 == it->IB3)
{
int f = IsIntersect(it->IB1,it->IB3,it->IB2,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}
if (m_Line[i].b1 == it->IB2&&m_Line[i].b2 == it->IB1)
{
int f = IsIntersect(it->IB2,it->IB1,it->IB3,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}
if (m_Line[i].b1 == it->IB2&&m_Line[i].b2 == it->IB3)
{
int f = IsIntersect(it->IB2,it->IB3,it->IB1,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}
if (m_Line[i].b1 == it->IB3&&m_Line[i].b2 == it->IB2)
{
int f = IsIntersect(it->IB3,it->IB2,it->IB1,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}
if (m_Line[i].b1 == it->IB3&&m_Line[i].b2 == it->IB1)
{
int f = IsIntersect(it->IB3,it->IB1,it->IB2,pt);
if (f == 0)
{
Ptdata[0] = p1;
Ptdata[1] = p2;
Ptdata[2] = p3;
return 1;
}
}

}
return 0;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// FindCloseLine
// Parameters:
// AcGePoint3d pt 插入点
// Returns:
// 返回最近的直线在数

组里的索引
// description:
// 找到最近的一条边
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::FindCloseLine(AcGePoint3d pt)
{
int i = 0;
int f; //标识是哪一条直线最近
vector ::iterator it = m_Line.begin();
for (;it!=m_Line.end();it++)
{
Ps pt1,pt2;
float distance,dmin=0;
pt1 = FindPt(it->b1);
pt2 = FindPt(it->b2);
AcGePoint3d point1(pt1.X,pt1.Y,0);
AcGePoint3d point2(pt2.X,pt2.Y,0);
AcDbLine *pline1 =new AcDbLine(point1,point2);
//////////////////////////////////////////////////////////////////////////
//块表记录
AcDbBlockTable *pBlockTable;
acdbHostApplicationServices()->workingDatabase()->getSymbolTable(pBlockTable,AcDb::kForRead);
AcDbBlockTableRecord *pBlockRecord;
pBlockTable->getAt(ACDB_MODEL_SPACE,pBlockRecord,AcDb::kForWrite);
AcDbObjectId lineID;
pBlockRecord->appendAcDbEntity(lineID,pline1);
pBlockRecord->close();
pBlockTable->close();
//////////////////////////////////////////////////////////////////////////
//找到直线上到插入点的最近点
AcGePoint3d ptClose;
pline1->getClosestPointTo(pt,ptClose);
distance = acutDistance(asDblArray(pt),asDblArray(ptClose));
pline1->close();
//////////////////////////////////////////////////////////////////////////
DleteLine(lineID);
if (dmin == 0)
{
dmin = distance;
f = i;
}
if (distance < dmin)
{
dmin = distance;
f = i;
}
i++;
}
return f;
}

// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// IsIntersect
// Parameters:
// int id1,int id2,int id3,AcGePoint3d pt
// Returns:
// 1 or 0 有交点返回1，没有返回0.
// description:
// 判断插入点与第三点在平面上的投影的点的连线是否与插入点的最近边相交
// -----------------------------------------------------------------------------------------------------------------------
int CMstruct::IsIntersect(int id1,int id2,int id3,AcGePoint3d pt)
{
Ps pt1,pt2,pt3;
pt1 = FindPt(id1);
pt2 = FindPt(id2);
pt3 = FindPt(id3);
//////////////////////////////////////////////////////////////////////////
AcGePoint3dArray parray;
AcGePoint3d point1(pt1.X,pt1.Y,0);
AcGePoint3d point2(pt2.X,pt2.Y,0);
AcGePoint3d point3(pt3.X,pt3.Y,0);
AcDbLine *pline1 = new AcDbLine(point1,point2);
AcDbLine *pline2 = new AcDbLine(point3,pt);
AcDbObjectId ObjId1,ObjId2;
AcDbBlockTable *pTable;
acdbHostApplicationServices()->workingDatabase()->getSymbolTable(pTable,AcDb::kForRead);
AcDbBlockTableRecord *pTableRec;
pTable->getAt(ACDB_MODEL_SPACE,pTableRec,AcDb::kForWrite);
pTableRec->appendAcDbEntity(ObjId1,pline1);
pTableRec->appendAcDbEntity(ObjId2,pline2);
//////////////////////////////////////////////////////////////////////////
////看直线是

否有交点
((AcDbEntity *&)pline1)->intersectWith((AcDbEntity *&)pline2,AcDb::kOnBothOperands,parray);
pline1->close();
pline2->close();
pTableRec->close();
pTable->close();
DleteLine(ObjId1);
DleteLine(ObjId2);
if (parray.length()>0)
{
return 1;
}
return 0;
}


// -----------------------------------------------------------------------------------------------------------------------
// Summary:
// DleteLine
// Parameters:
// AcDbObjectId &id
// Returns:
// void
// description:
// 删除已画的直线
// -----------------------------------------------------------------------------------------------------------------------
void CMstruct::DleteLine(AcDbObjectId &id)
{
AcDbLine *pline = new AcDbLine;
acdbOpenAcDbEntity((AcDbEntity *&)pline,id,AcDb::kForWrite);
pline->erase();
pline->close();
}