等高線建模核心演算法

lightmare625發表於2018-07-19

構造斷層面 

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.三維地學建模-》資料匯入 -》平面地質圖匯入

 

 

 

相關文章