等高线建模核心算法
构造断层面
void TestBuildSurf()
{
//打开线面要素类
GeoClass geoFace,geoLine,geofault;
CSFeatureCls *pFaceCls = NULL, *pLineCls = NULL,*pFaultCls = NULL;
//等高面
geoFace.Open( "gdbp://MapGisLocal/plugintest/sfcls/isoheightsurf" );
geoFace.Get( (void**)&pFaceCls, "CSFeatureCls" );
//用户的线数据
geoLine.Open( "gdbp://MapGisLocal/plugintest/sfcls/line" );
geoLine.Get( (void**)&pLineCls, "CSFeatureCls" );
//断层面
geofault.Open( "gdbp://MapGisLocal/plugintest/sfcls/faultSurf" );
geofault.Get( (void**)&pFaultCls, "CSFeatureCls" );
if( pFaceCls==NULL || pLineCls==NULL || pFaultCls==NULL )
return;
pFaceCls->cls_Clear();
pFaultCls->cls_Clear();
//获取等高线和断层线
vector<CAnyLine*> lines, faultLines;
vector<double> heights;
//矩形范围查找要素集
CSFeatureSet fset;
pLineCls->RectSelect( -1, NULL, &fset );
//遍历要素 获取其中的线要素集
long flag = fset.MoveFirst();
while( flag!=END_OF_SET )
{
TYPE_OBJ_ID oid = 0;
CAnyLine *pLine = new CAnyLine();
CRecord rcd;
char lineType = 0;//线类型
short isNull = 0;
double dValue = 0;//高程值
fset.GetObjID( &oid );
fset.line_Get( MOVE_CURRENT, pLine, &rcd );//id光标,几何信息,属性信息
rcd.GetCharVal( "线类型", lineType, isNull );//取字段值(根据字段名称)
if( lineType==1 )
{
lines.push_back( pLine );//等高线
rcd.GetDoubleVal( "高程值", dValue, isNull );
heights.push_back(dValue/10);
}
else if( lineType==2)
faultLines.push_back( pLine );//断层线
else
delete pLine;
flag = fset.MoveNext();
}
// 统计高程范围
double dMax = -1e+17, dMin = 1e+17;
for(int i=0; i<heights.size(); i++ )
{
if( heights[i]>dMax )
dMax = heights[i];
if( heights[i]<dMin )
dMin = heights[i];
}
dMax += 1.0;
dMin -= 1.0;
//构建断层切割面
vector<CAnySurface*> faultFaces;
for( int i=0; i<faultLines.size(); i++ )
{
CAnySurface *pSurface = new CAnySurface();
vector<D_3DOT> vectorDot3D;
vector<Tri> vectorTri;
int dotNum = faultLines[i]->m_varLin.dotNum()-1;//封闭折线上的点多了一个重复点。所以减1
// 重新定义大小
vectorDot3D.resize(dotNum * 2);//二维点拉伸成三维的面
vectorTri.resize(dotNum * 2);//三角网的个数=2*二维点个数
D_3DOT dotTemp;
Tri triTemp;
for( int dotIndex=0; dotIndex < dotNum; dotIndex++ )
{
//更新三维点的坐标
dotTemp.x = ((D_DOT*)faultLines[i]->m_varLin.ptXY())[dotIndex].x;//二维坐标合集
dotTemp.y = ((D_DOT*)faultLines[i]->m_varLin.ptXY())[dotIndex].y;
dotTemp.z = dMin;
vectorDot3D[dotIndex] = dotTemp;//封闭线最低点坐标
dotTemp.z = dMax;
vectorDot3D[dotIndex + dotNum ] = dotTemp; //封闭线最高点坐标
//更新三角网的坐标(画图,重点)
if( dotIndex==dotNum - 1 )
triTemp.a = dotNum;
else
triTemp.a = dotIndex + 1 + dotNum;
triTemp.b = dotIndex + dotNum;
triTemp.c = dotIndex;
vectorTri[dotIndex] = triTemp;
if( dotIndex==dotNum - 1 )
{
triTemp.a = 0;
triTemp.b = dotNum;
}
else
{
triTemp.a = dotIndex + 1;
triTemp.b = dotIndex + 1 + dotNum ;
}
triTemp.c = dotIndex;
vectorTri[dotIndex + dotNum ] = triTemp;
}
//设置面的点集 三角网集
pSurface->Set( dotNum*2, &vectorDot3D.front(), dotNum*2, (DWORD*)&vectorTri.front() );
//断层面
faultFaces.push_back(pSurface);
//更新面要素
pFaultCls->surface_Append( pSurface );
}
// 获取高程样本点,构建三角网tin
vector<D_3DOT> sampleDots;
for( int i=0; i<lines.size(); i++ )
{
for( int j=0; j<lines[i]->m_varLin.dotNum()-1; j++ )
{
D_3DOT dotTemp;
dotTemp.x = ((D_DOT*)lines[i]->m_varLin.ptXY())[j].x;
dotTemp.y = ((D_DOT*)lines[i]->m_varLin.ptXY())[j].y;
dotTemp.z = heights[i];
sampleDots.push_back( dotTemp );
}
}
// 利用断层面切割tin,并保存
long PntNum = 0, TriNum = 0;
D_3DOT *pTriPnt = NULL;
DWORD *pTri = NULL;
DWORD *pTriTopo = NULL;
tnCreateDelaunayTinByMemPnts( NULL, NULL, &sampleDots.front(), sampleDots.size(), NULL, 0, 0, false,
PntNum, &pTriPnt, TriNum, &pTri, &pTriTopo );
CAnySurface tinS;
tinS.Set( PntNum, pTriPnt, TriNum, pTri, pTriTopo );
pFaceCls->surface_Append( &tinS );
tnFreeTinPntNet( pTriPnt, pTri, pTriTopo );
// 释放内存
for( int i=0; i<lines.size(); i++ )
delete lines[i];
for( int i=0; i<faultLines.size(); i++ )
delete faultLines[i];
for( int i=0; i<faultFaces.size(); i++ )
delete faultFaces[i];
}
切割面
void CutSurf()
{
GeoClass geoFaceLine, geoFaceFault, geoFaceResult;
CSFeatureCls *pFaceLineCls = NULL, *pFaceFaulCls = NULL, *pFaceResultCls = NULL;
CAnySurface faceLine, faceFault;
CMultiSurface resultFaces;
//打开基本要素类
geoFaceLine.Open( "gdbp://MapGisLocal/plugintest/sfcls/isoheightsurf" );
geoFaceLine.Get( (void**)&pFaceLineCls, "CSFeatureCls" );
geoFaceFault.Open( "gdbp://MapGisLocal/plugintest/sfcls/faultSurf" );
geoFaceFault.Get( (void**)&pFaceFaulCls, "CSFeatureCls" );
geoFaceResult.Open( "gdbp://MapGisLocal/plugintest/sfcls/resultSurf" );
geoFaceResult.Get( (void**)&pFaceResultCls, "CSFeatureCls" );
if( pFaceLineCls==NULL || pFaceFaulCls==NULL || pFaceResultCls==NULL )
return;
pFaceResultCls->cls_Clear();
//获取三维面
pFaceLineCls->surface_Get( 1, &faceLine );//高程面
pFaceFaulCls->surface_Get( 1, &faceFault );//断层面
int i=faceFault.GetPointNum();
int j=faceLine.GetPointNum();
//曲面与曲面相互切割重构
G3DXCutSurfaceWithSurface(&faceLine, &faceFault, &resultFaces);
int item=resultFaces.GetNum();
for( long i=0; i<resultFaces.GetNum(); i++ )
{
CAnySurface *pSurf = resultFaces.Get(i);
pFaceResultCls->surface_Append( pSurf );
}
geoFaceLine.Close();
geoFaceFault.Close();
geoFaceResult.Close();
}
构建地址面 二维转三维
1.三维地学建模-》数据导入 -》平面地质图导入