等高线建模核心算法

构造断层面 

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.三维地学建模-》数据导入 -》平面地质图导入

 

 

 

posted @ 2018-07-19 09:45  lightmare  阅读(482)  评论(0编辑  收藏  举报