flow_parallel.cpp

这是一个GPS相关的程序

///////////////////////////////////////////////////////////
//														 //
//				Deterministic Infinity					 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_DInf(	int x, int y )
{
	int		Direction, x1, y1, x2, y2;

	double	Slope, Aspect, z;

	Get_Gradient(x, y, Slope, Aspect);

	if( Aspect >= 0.0 )
	{
		Direction	= (int)(Aspect / M_PI_045);
		Aspect		=  fmod(Aspect , M_PI_045) / M_PI_045;

		z			= pDTM->asDouble(x, y);
		x1			= Get_xTo(Direction + 0, x);
		y1			= Get_yTo(Direction + 0, y);
		x2			= Get_xTo(Direction + 1, x);
		y2			= Get_yTo(Direction + 1, y);

		if( (!is_InGrid(x1, y1) || z > pDTM->asDouble(x1, y1))
		&&	(!is_InGrid(x2, y2) || z > pDTM->asDouble(x2, y2)) )
		{
			Add_Fraction(x, y,  Direction			, 1.0 - Aspect);
			Add_Fraction(x, y, (Direction + 1) % 8	,       Aspect);
		}
		else
		{
			Add_Fraction(x, y, pDTM->Get_Gradient_NeighborDir(x, y));
		}
	}
}


///////////////////////////////////////////////////////////
//														 //
//				Multiple Flow Direction					 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CFlow_Parallel::Set_MFD(	int x, int y )
{
	int		i, ix, iy;

	double	z, d, dzSum, dz[8];

	z		= pDTM->asDouble(x, y);
	dzSum	= 0.0;

	for(i=0; i<8; i++)
	{
		ix		= Get_xTo(i, x);
		iy		= Get_yTo(i, y);

		if( pDTM->is_InGrid(ix, iy) )
		{
			d		= z - pDTM->asDouble(ix, iy);
		}
		else
		{
			ix		= Get_xTo(i + 4, x);
			iy		= Get_yTo(i + 4, y);

			if( pDTM->is_InGrid(ix, iy) )
			{
				d		= pDTM->asDouble(ix, iy) - z;
			}
			else
			{
				d		= 0.0;
			}
		}

		if( d > 0.0 )
		{
			dzSum	+= (dz[i]	= pow(d / Get_Length(i), MFD_Converge));
		}
		else
		{
			dz[i]	= 0.0;
		}
	}

	if( dzSum )
	{
		for(i=0; i<8; i++)
		{
			if( dz[i] )
			{
				Add_Fraction(x, y, i, dz[i] / dzSum );
			}
		}
	}
}


///////////////////////////////////////////////////////////
//														 //
//				Braunschweiger Reliefmodell				 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Flow_BRM.h"

#define BRM_nint(x)		((int)(x >= 0 ? x + 0.5 : x - 0.5))


//---------------------------------------------------------
void CFlow_Parallel::Set_BRM(	int x, int y )
{
	int		Dir, QBinaer,
			ix[3], iy[3],
			nexp[6];

	double	QLinks, QMitte, QRecht,
			nnei[6];

	//-----------------------------------------------------
	if( is_InGrid(x, y, 1) ) // Rand !!!
	{
		Dir	= BRM_InitRZ(x,y,ix,iy);

		if( Dir >= 0 )
		{
			if( Dir % 2 )
			{
				BRM_GetDiago(Dir,x,y,ix,iy,nnei,nexp);
				BRM_QStreuung(4,1,nnei,nexp,QBinaer,QLinks,QMitte,QRecht);
			}
			else
			{
				BRM_GetOrtho(Dir,x,y,ix,iy,nnei,nexp);
				BRM_QStreuung(6,0,nnei,nexp,QBinaer,QLinks,QMitte,QRecht);
			}

			Add_Fraction(x,y,(Dir+1)%8,BRM_BitMtrx[0][QBinaer] ? QLinks : 0);
			Add_Fraction(x,y,(Dir+0)%8,BRM_BitMtrx[1][QBinaer] ? QMitte : 0);
			Add_Fraction(x,y,(Dir+7)%8,BRM_BitMtrx[2][QBinaer] ? QRecht : 0);
		}
	}
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_Init(void)
{
	int 	i;

	double	DXT	= Get_Cellsize()/2,
			DYT	= Get_Cellsize()/2;

	//-----------------------------------------------------
	BRM_kgexp[0]	= (int)(atan2(DXT     , Get_Cellsize()) * M_RAD_TO_DEG);
	BRM_kgexp[1]	= (int)(atan2(Get_Cellsize(), DYT     ) * M_RAD_TO_DEG) + 1;
	BRM_kgexp[2]	= (int)(atan2(Get_Cellsize(),-DYT     ) * M_RAD_TO_DEG);
	BRM_kgexp[3]	= (int)(atan2(DXT     ,-Get_Cellsize()) * M_RAD_TO_DEG) + 1;

	for(i=0; i<4; i++)
		BRM_kgexp[i+4]	= BRM_kgexp[i] + 180;

	//-----------------------------------------------------
	for(i=0; i<=360; i++)
	{
		BRM_sinus[i]	= -sin(i * M_DEG_TO_RAD);
		BRM_cosin[i]	= -cos(i * M_DEG_TO_RAD);
	}

	//---BRM_idreh---------------------------------------------
	BRM_idreh[0]	= 180;
	BRM_idreh[1]	= 180 -	BRM_nint(atan2(Get_Cellsize(), Get_Cellsize()) * M_RAD_TO_DEG);
	BRM_idreh[2]	=  90;
	BRM_idreh[3]	=		BRM_nint(atan2(Get_Cellsize(), Get_Cellsize()) * M_RAD_TO_DEG);
	BRM_idreh[4]	=   0;

	for(i=1; i<4; i++)
		BRM_idreh[i+4]	= BRM_idreh[i] + 180;
}

//---------------------------------------------------------
int CFlow_Parallel::BRM_InitRZ(int x, int y, int ix[3], int iy[3])
{
	int		i, j, Dir;

	double	Slope, Aspect;

	Get_Gradient(x, y, Slope, Aspect);

	Aspect	*= M_RAD_TO_DEG;

	if( Aspect < 0 )
	{
		return( -1 );
	}

	//---Kategorisierte-Exposition-------------------------
	Dir	= 0;
	while( Aspect > BRM_kgexp[Dir] && Dir < 8 )
		Dir++;
	Dir	%=8;

    //---Finde-Die-3-ZielRasterZellen----------------------
	for(i=0; i<3; i++)	// zxy[]: 0=Recht, 1=Mitte, 2=Links
	{
		j		= (Dir + 7 + i) % 8;
		ix[2-i]	= Get_xTo(j,x);
		iy[2-i]	= Get_yTo(j,y);
	}

	return(Dir);
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_GetOrtho(int Dir, int x, int y, int ix[3], int iy[3], double nnei[6], int nexp[6])
{
	int		jx, jy, i,
			i0	= (Dir + 2) % 8,
			i4	= (Dir + 6) % 8;

	double	Slope, Aspect;

	for(i=0; i<3; i++)
	{
		jx		= ix[i];
		jy		= iy[i];

		Get_Gradient(jx, jy, Slope, Aspect);

		nnei[i]	= M_RAD_TO_DEG * Slope;
		nexp[i]	= (int)(M_RAD_TO_DEG * Aspect);
	}

	jx		= Get_xTo(i0,x);
	jy		= Get_yTo(i0,y);

	Get_Gradient(jx, jy, Slope, Aspect);

	nnei[3]	= M_RAD_TO_DEG * Slope;
	nexp[3]	= (int)(M_RAD_TO_DEG * Aspect);

	jx		= Get_xTo(i4,x);
	jy		= Get_yTo(i4,y);

	Get_Gradient(jx, jy, Slope, Aspect);

	nnei[5]	= M_RAD_TO_DEG * Slope;
	nexp[5]	= (int)(M_RAD_TO_DEG * Aspect);

	Get_Gradient(x, y, Slope, Aspect);

	nnei[4]	= M_RAD_TO_DEG * Slope;
	nexp[4]	= (int)(M_RAD_TO_DEG * Aspect);	//[jy][jx]) ????!!!!...;

	for(i=0; i<6; i++)
		if(nexp[i]<0)
			nexp[i]	= nexp[4];

	for(i=0; i<6; i++)
	{
		nexp[i]	+= BRM_idreh[Dir];

		if(nexp[i]>360)
			nexp[i]	-= 360;
	}
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_GetDiago(int Dir, int x, int y, int ix[3], int iy[3], double nnei[6], int nexp[6])
{
	int		i;

	double	Slope, Aspect;

	Get_Gradient( x   ,  y   , Slope, Aspect);
	nexp[0]	= (int)(M_RAD_TO_DEG * Aspect);
	nnei[0]	= M_RAD_TO_DEG * Slope;

	Get_Gradient(ix[0], iy[0], Slope, Aspect);
	nexp[1]	= (int)(M_RAD_TO_DEG * Aspect);
	nnei[1]	= M_RAD_TO_DEG * Slope;

	Get_Gradient(ix[2], iy[2], Slope, Aspect);
	nexp[2]	= (int)(M_RAD_TO_DEG * Aspect);
	nnei[2]	= M_RAD_TO_DEG * Slope;

	Get_Gradient(ix[1], iy[1], Slope, Aspect);
	nexp[3]	= (int)(M_RAD_TO_DEG * Aspect);
	nnei[3]	= M_RAD_TO_DEG * Slope;

	for(i=1; i<4; i++)
		if(nexp[i]<0)
			nexp[i]	= nexp[0];

	for(i=0; i<4; i++)
	{
		nexp[i]	+= BRM_idreh[Dir];

		if(nexp[i]>360)
			nexp[i]	-= 360;
	}
}

//---------------------------------------------------------
void CFlow_Parallel::BRM_QStreuung(int i64, int g64, double nnei[6], int nexp[6], int &QBinaer, double &QLinks, double &QMitte, double &QRecht)
{
	int		i, j, ix, iy,
			ALinks, AMitte=2, ARecht;

	double	x=1, y=1, sg=0, a,
			s[6], c[6];

	ALinks	= ARecht	= 0;
	QLinks	= QRecht	= 0.0;

	for(i=0; i<i64; i++)
		sg	+= nnei[i];

	sg	= i64 / sg;

	for(i=0; i<i64; i++)
	{
		a		= sg * nnei[i];
		j		= nexp[i];
		s[i]	= a * BRM_sinus[j];
		c[i]	= a * BRM_cosin[j];
	}

    //---QLinks-ermitteln----------------------------------
	for(i=0; i<100; i++)
	{
		ix	= BRM_nint(x) - 1;
		iy	= BRM_nint(y) - 1;

		for(j=0; j<i64; j++)
		{
			a		= BRM_g[g64][j][ix][iy];
			x		+= s[j] * a;
			y		+= c[j] * a;
		}	 

		if(x<1)
		{
			ALinks	= 0;
			QLinks	= 0;
			break;
		}

		if( x>8.99 || y<1 )
		{
			ALinks	= 4;
			QLinks	= 1;			
			break;
		}

		if(y>8.95)
		{
			if(x<1.02)
			{
				ALinks	= 0;
				QLinks	= 0;
			}
			else
			{
				ALinks	= 4;
				QLinks	= (x - 1) / 8;

				if(i64==4)
				{
					if(QLinks<0.5)
						QLinks	= QLinks * (1.67 - QLinks * 1.078);
					else
						QLinks	= QLinks * 0.869 + 0.131;
				}
			}
			break;
		}
	}


    //---QRechts-ermitteln---------------------------------
	x	= 9;
	y	= 1;

	for(i=0; i<100; i++)
	{
		ix	= BRM_nint(x) - 1;
		iy	= BRM_nint(y) - 1;

		for(j=0; j<i64; j++)
		{
			a	= BRM_g[g64][j][ix][iy];
			x	+= s[j] * a;
			y	+= c[j] * a;
		}

		if(x>9)
		{
			ARecht	= 0;
			QRecht	= 0;
			break;
		}

		if( x<1.01 || y<1 )
		{
			ARecht	= 1;
			QRecht	= 1;
			break;
		}

		if(y>8.95)
		{
			if(x>8.98)
			{
				ARecht	= 0;
				QRecht	= 0;
			}
			else
			{
				ARecht	= 1;
				QRecht	= 1 - (x - 1) / 8;

				if(i64==4)					
				{
					if(QRecht<0.5)
						QRecht	= QRecht * (1.67 - QRecht * 1.078);
					else
						QRecht	= QRecht * 0.869 + 0.131;
				}
			}
			break;
		}
	}

	//---QMitte-ist-Rest-/-QBinaer'Bits-setzen-------------
	QMitte	= 1 - QLinks - QRecht;
	QBinaer	= ALinks + AMitte + ARecht;

	if(QMitte<=0.01)
	{
		a		= QLinks + QRecht;
		QLinks	/= a;
		QRecht	/= a;
		QMitte	= 0;
		QBinaer	-= 2;
	}
}