contourtracer.cpp

二维规则网格数据的等值线追踪

// ContourTracer.cpp: implementation of the CContourTracer class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "ContourGenerator.h"
#include "ContourTracer.h"
#include "2DMemAllocator.h"

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////


CContourTracer::CContourTracer()
{
	m_ppGridData = NULL;

	m_pCurveList = NULL;
	m_currentCurveLine = NULL;

	xSide = NULL;
	ySide = NULL;

}

CContourTracer::~CContourTracer()
{
	FreeMemory();	
}



BOOL CContourTracer::ExecuteTracing(float value)
{
	ASSERT( m_gridDataInfo.cols != 0 ); //在调用该函数之前必须调用SetGridDataInfo()函数设置网格数据信息

	if( value < m_gridDataInfo.zMin && value > m_gridDataInfo.zMax )
	{
		return FALSE;
	}


	m_valueTracing = value;


	//1.为xSide和ySide分配内存空间
	AllocateMemory();

	//2.扫描网格纵横边，内插等值点
	//该函数之中，在计算等值点时，发现追踪值与网格点上的数据相同时，
	//在计算前，会都数据做修正(加上一小偏移值)
	InterpolateTracingValue(); 

	//3.先追踪开曲线
	TracingNonClosedContour();

	//4.再追踪闭曲线
	TracingClosedContour();

	//5.释放空间
	//FreeMemory();

	return TRUE;
	



	
}

void CContourTracer::AllocateMemory()
{// 分配xSide，ySide空间
	int cols = m_gridDataInfo.cols;
	int rows = m_gridDataInfo.rows;
	
	if( xSide == NULL )
	{
		//网格中存在rows*(cols-1)条横边，所有需要为xSide分配rows*(cols-1)空间就行

		C2DMemAllocator::AllocMemory2D(xSide,rows,cols-1/*not cols*/);
	}
	
	if( ySide == NULL )
	{
		//网格中存在(rows-1)*cols条纵边，所有需要为ySide分配(rows-1)*cols空间就行

		C2DMemAllocator::AllocMemory2D(ySide,rows-1/*not rows*/,cols);
	}
	
}

void CContourTracer::FreeMemory()
{
	
	if( xSide != NULL )
	{
		C2DMemAllocator::FreeMemory2D(xSide);
		xSide = NULL;
	}
	
	if( ySide != NULL )
	{
		C2DMemAllocator::FreeMemory2D(ySide);
		ySide = NULL;
	}
	
}


//-----消除网格交点奇异Z值------
//在自动网格边上追踪等值点时，有时会遇到网格交点上的Z值与追踪的等值线值相等的情况;
//如果不对这种情况进行处理，网格追踪算法就会出现错误，一般的处理如下：
//将网格交点上的z值加或减一个修正值，来消除奇异点的影响

//void CContourTracer::AmendingData(BOOL  bForTracing  /*是否为追纵算法而修正数据*/, 
//							   float fValueTracing/*所要追踪的等值线值*/)
//{
//
//	ASSERT( m_ppGridData!=NULL );
//
//	int i,j;
//
//	float shift = 0.001f;  //修正值
//
//	if(bForTracing)
//	{
//		//若把与所要追踪的值相同的加修正值
//		for(i=0; i<m_gridDataInfo.rows; i++)
//		{
//			for(j=0; j<m_gridDataInfo.cols; j++)
//			{
//				if( m_ppGridData[i][j] == fValueTracing )  
//					m_ppGridData[i][j] += shift;    //should changing back!!
//			}
//		}
//	}
//	else
//	{
//		float shiftVal = fValueTracing + shift;
//
//		for(i=0; i<m_gridDataInfo.rows; i++)
//		{
//			for(j=0; j<m_gridDataInfo.cols; j++)
//			{
//				if( m_ppGridData[i][j] == shiftVal )  
//					m_ppGridData[i][j] -= shift;    //restore original data
//			}
//		}
//	}
//}

//扫描网格的纵、横边，并线性插值计算等值点的情况
//将各边上的等值点情况存储于xSide和ySide数组中，
// xSide存储所有横边上的等值线情况, ySide存储所有纵边上的等值点情况
//在插值计算时，对『与追踪值相等的数据』要进行修正处理后才计算，但在做修正处理时不要改变原来的数据
void CContourTracer::InterpolateTracingValue()
{
/*      网格点标识如下:
        
            (i+1,j)·--------·(i+1,j+1)
                    |        |
                    |        |
                    |        |
	                |        |
	         (i,j) ·--------·(i,j+1)

              i:表示行号(向上增加)
			  j:表示列号(向右增加)
			  标识一个网格交点时，行号在前，列号在右，如：(i,j)
*/
/*        xSide,ySide中存储r值，(w为追踪值)
          对于网格横边，r = (w - pData[i][j]) / (pData[i][j+1]-pData[i][j]);
		  对于网格纵边，r = (w - pData[i][j]) / (pData[i+1][j]-pData[i][j]);

		  由于浮点运算的误差，xSide[i][j],ySide[i][j]有可能等于1.0或0.0 
		  考虑如下情况：
		  1。当追踪值与网格点上的值很接近(但不相等)时，由于运算误差，就会等于1.0
		     比如追踪0值时，遇到如下边:
			       20 ·--------·-0.00000016   此边上有0值，但计算 (0-20)/(-0.00000016-20) == 1.0 

       
		  2。当网格边上两端点上的值相差很悬殊时。
		     比如追踪2值，遇到如下边：
		     1.70141E+038 ·--------·1   此边上有2值，计算(2-1.70141E+038) / (1-1.70141E+038) == 1.0 

        网格边上有等值点时，理论上比例值不会等于0或1；
		 但由于计算误差，我们在算法中判断时，认为0或1也是有等值点的

        所以xSide,ySide中存储的值是[0,1]的闭区间，不是(0,1)的开区间
*/
	ASSERT( m_ppGridData!=NULL );
	ASSERT( xSide != NULL );
	ASSERT( ySide != NULL );
	
	int i,j;
	
	int   rows = m_gridDataInfo.rows;
	int   cols = m_gridDataInfo.cols;
	float  w   = m_valueTracing; 
	float** pData = m_ppGridData;
	
	float H1,H2; //分别记录一条边的两个点上的数据值
	float flag; 
	float shift = 0.001f;  //修正值


	/* 扫描并计算横边上的等值点,有rows*(cols-1)条横边需要扫描*/
	for(i=0; i<rows ; i++)
	{
		for(j=0;j<cols-1;j++)
		{				
			/*考查横边(i,j)上的左交点(i,j)上的值pData[i][j]
			                 和右交点(i,j+1)的值pData[i][j+1]*/
//			if( pData[i][j] == pData[i][j+1] )
//			{
//				// -2表示此边无等值点
//				xSide[i][j] = -2.0f;
//			}
//			else
//			{
//				xSide[i][j] = (w - pData[i][j]) / ( pData[i][j+1] - pData[i][j] );
//				
//				if( ( xSide[i][j] <= 0 )  || ( xSide[i][j] >= 1 ) )
//					xSide[i][j] = -2.0f;
//			}

			H1 = pData[i][j]; 			H2 = pData[i][j+1];		

			if( H1 == H2 )
			{
				xSide[i][j] = -2.0f;
			}
			else
			{
				flag = (w-H1) * (w-H2);
				if( flag > 0 )
				{
					xSide[i][j] = -2.0f;
				}
				else if( flag < 0 )
				{
					xSide[i][j] = (w-H1) / (H2-H1) ; 
					ASSERT(xSide[i][j]>=0 && xSide[i][j]<=1.0f);
				}
				else if( flag == 0)
				{//其中有一值与追踪值w相等，则修正之(加上一小值偏移量)
					if( H1 == w )
					{
						H1 += shift;
					}
					else
					{
						H2 += shift;
					}
					
					xSide[i][j] = (w-H1) / (H2-H1) ; 
				//	ASSERT(xSide[i][j]>=0 && xSide[i][j]<=1.0f);

				}
				else
				{
					ASSERT(FALSE);
				}

			}
			
			


		}
	}

	/* 扫描并计算纵边上等值点,有(rows-1)*cols条纵边需要扫描*/
	for(i=0; i<rows-1;i++)
	{
		for(j=0; j<cols; j++)
		{			
			/*考查纵边(i,j)上的下交点(i,j)上的值pData[i][j]
			                 和上交点(i+1,j)的值pData[i+1][j]*/

//			if( pData[i][j] == pData[i+1][j] )
//			{
//				/* -2表示此边无等值点，或已追踪过，以后不再考虑*/
//				ySide[i][j] = -2.0f;
//			}
//			else
//			{
//				ySide[i][j] = (w - pData[i][j]) / ( pData[i+1][j] - pData[i][j] );
//				
//				if( ( ySide[i][j] <= 0 )  || ( ySide[i][j] >= 1 ) )
//					ySide[i][j] = -2.0f;
//			}
                                             
			H1 = pData[i][j];            
			H2 = pData[i+1][j];		     
                                         
			if( H1 == H2 )               
			{
				ySide[i][j] = -2.0f;
			}
			else
			{
				flag = (w-H1) * (w-H2);
				if( flag > 0 )
				{
					ySide[i][j] = -2.0f;
				}
				else if( flag < 0 )
				{/*
				 网格边上有等值点时，存储的值∈[0,1] (闭区间)
                 */
					ySide[i][j] = (w-H1) / (H2-H1) ; 
					ASSERT(ySide[i][j]>=0 && ySide[i][j]<=1.0f);
				}
				else if( flag == 0 )
				{//如果其中有一值与追踪值w相等，则修正之(加上一小值偏移量)
					if( H1 == w )
					{
						H1 += shift;
					}
					else
					{
						H2 += shift;
					}
					
					ySide[i][j] = (w-H1) / (H2-H1) ; 
				//	ASSERT(ySide[i][j]>=0 && ySide[i][j]<=1.0f);
				}
				else
				{
					ASSERT(FALSE);
				}


			}

		}
	}
}

inline void CContourTracer::CalcAndSaveOnePointCoord(int i, int j, BOOL bHorizon,float &x, float &y)
{
	/*static*/ float deltX = (m_gridDataInfo.xMax - m_gridDataInfo.xMin) / ( m_gridDataInfo.cols - 1 );
	/*static*/ float deltY = (m_gridDataInfo.yMax - m_gridDataInfo.yMin) / ( m_gridDataInfo.rows - 1 );
	// should not be static, because diffirent input data has diffirent xyMin Max
	
	
	if( bHorizon )
	{//在横边上
		x = m_gridDataInfo.xMin + ( j + xSide[i][j] ) * deltX;
		y = m_gridDataInfo.yMin + i * deltY;		
	}
	else
	{//在纵边上
		x = m_gridDataInfo.xMin + j * deltX;
		y = m_gridDataInfo.yMin + ( i + ySide[i][j] ) * deltY;
	}

	//Saving Coord
	CGeoPoint point(x,y);
	m_currentCurveLine->Add(point);

}

//当下一个等值点找到后做相应的处理
void CContourTracer::HandlingAfterNextPointFounded(int i, int j, BOOL bHorizon)
{//参数说明：i,j分别是等值点所在边的编号，bHorizon指明所在边是横边还是纵边

	//当下一个等值点找到后做相应的处理,如下：
	//1.记录该等值点的i,j
	//2.计算并保存该等值点的坐标
	//3.标志该等值点所在边的已经搜索过

	//验证i∈[0,rows-1], j∈[0,cols-1]
	ASSERT( i>=0 && i<=m_gridDataInfo.rows-1 && j>=0 && j<=m_gridDataInfo.cols-1 );

	//1.
	NextPoint.i = i;
	NextPoint.j = j;
    NextPoint.bHorV = bHorizon;

	//2.
	CalcAndSaveOnePointCoord(i,j,bHorizon,NextPoint.x,NextPoint.y);

	//3.
	if( NextPoint.bHorV )
	{
		xSide[i][j] = -2.0f; //已经追踪过
	}
	else
	{	
		ySide[i][j] = -2.0f; //已经追踪过
	}

}

void CContourTracer::TracingNextPoint()
{
/*  
  1.先确定出等值线的前进方向(自下向上、由左向右、自上向下、由右向左，其中之一)
  2.再追踪下一个等值点

  前进方向可以如下判定：
    if( 当前点.行号 > 前一点.行号 )
    {
	    下---->上
	}
	else if( 当前点.列号 > 前一点.列号 )
    {
		左---->右
	}
    else if( 当前点在横边上 )
	{
		上---->下
	}
	else
	{
	     右---->左