morphometry.cpp

这是一个GPS相关的程序

		Set_Parameters_NoData(x, y, true);
	}
	else
	{
		z				= pDTM->asDouble(x, y);
		SubMatrix[4]	= 0.0;

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

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

				if( pDTM->is_InGrid(ix, iy) )
				{
					SubMatrix[iSub[i]]	= z - pDTM->asDouble(ix, iy);
				}
				else
				{
					SubMatrix[iSub[i]]	= 0.0;
				}
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
inline bool CMorphometry::Get_SubMatrix5x5(int x, int y, double SubMatrix[25])
{
	int		i, ix, iy, jx, jy;
	double	z;

	if( pDTM->is_NoData(x, y) )
	{
		Set_Parameters_NoData(x, y, true);
	}
	else
	{
		z	= pDTM->asDouble(x,y);

		for(i=0, iy=y-2; iy<=y+2; iy++)
		{
			jy	= iy < 0 ? 0 : ( iy >= Get_NY() ? Get_NY() - 1 : iy );

			for(ix=x-2; ix<=x+2; ix++, i++)
			{
				jx	= ix < 0 ? 0 : ( ix >= Get_NX() ? Get_NX() - 1 : ix );

				SubMatrix[i]	= pDTM->is_InGrid(jx, jy) ? pDTM->asDouble(jx, jy) - z : 0.0;
			}
		}

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//					The Methods							 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
// Maximum Slope (Travis et al., 1975, Peucker & Douglas, 1975))
//
// Travis, M.R., Elsner, G.H., Iverson, W.D., and Johnson, C.G. 1975:
//		VIEWIT: computation of seen areas, slope, and aspect for land-use planning.
//		USDA F.S. Gen. Tech. Rep. PSW-11/1975, 70p. Berkeley, California, U.S.A.
//
//---------------------------------------------------------
void CMorphometry::Do_MaximumSlope(int x, int y)
{
	int		i, ix, iy, j, Aspect;

	double	z, zm[8], Slope, Curv, hCurv, a, b;

	//-----------------------------------------------------
	if( pDTM->is_NoData(x, y) )
	{
		Set_Parameters_NoData(x, y, true);
	}
	else
	{
		//-------------------------------------------------
		z		= pDTM->asDouble(x, y);
        Slope	= Curv	= 0.0;

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

			if( !pDTM->is_InGrid(ix, iy) )
			{
				zm[i]	= 0.0;
			}
			else
			{
				zm[i]	= atan((z - pDTM->asDouble(ix, iy)) / Get_Length(i));
				Curv	+= zm[i];

				if( zm[i] > Slope )
				{
					Aspect	= i;
					Slope	= zm[i];
				}
			}
		}

		//-------------------------------------------------
		if( Aspect < 0.0 )
		{
			Set_Parameters_NoData(x, y);
		}
		else
		{
			//---------------------------------------------
			// Let's now estimate the plan curvature...

			for(i=Aspect+1, j=0, a=0.0; i<Aspect+8; i++, j++)
			{
				if( zm[i % 8] < 0.0 )
				{
					a	= j + zm[(i - 1) % 8] / (zm[(i - 1) % 8] - zm[i % 8]);
					break;
				}
			}

			if( a != 0.0 )
			{
				for(i=Aspect+7, j=0, b=0.0; i>Aspect; i--, j++)
				{
					if( zm[i % 8] < 0.0 )
					{
						b	= j + zm[(i + 1) % 8] / (zm[(i + 1) % 8] - zm[i % 8]);
						break;
					}
				}

				hCurv	=  45.0 * (a + b) - 180.0;
			}
			else
			{
				hCurv	=  180.0;
			}

			//---------------------------------------------
			Set_Parameters(x, y,
				Slope,
				Aspect * M_PI_045,
				Curv,
				zm[Aspect] + zm[(Aspect + 4) % 8],
				hCurv
			);
		}
    }
}

//---------------------------------------------------------
// Maximum Triangle Slope
//
// Tarboton, D.G. (1997):
//		'A new method for the determination of flow directions and upslope areas in grid digital elevation models',
//		Water Ressources Research, Vol.33, No.2, p.309-319
//
//---------------------------------------------------------
void CMorphometry::Do_Tarboton(int x, int y)
{
	int		i, ix, iy, j;
	double	z, zm[8], iSlope, iAspect, Slope, Aspect, G, H;

	//-----------------------------------------------------
	if( pDTM->is_NoData(x, y) )
	{
		Set_Parameters_NoData(x, y, true);
	}
	else
	{
		z		= pDTM->asDouble(x, y);

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

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

				if( pDTM->is_InGrid(ix, iy) )
				{
					zm[i]	=  z - (pDTM->asDouble(ix, iy) - z);
				}
				else
				{
					zm[i]	=  z;
				}
			}
		}

		//---------------------------------------------
        Slope	=  0.0;
		Aspect	= -1.0;

		for(i=0, j=1; i<8; i++, j=(j+1)%8)
		{
			if( i % 2 )	// i => diagonal
			{
				G		= (z		- zm[j]) / Get_Cellsize();
				H		= (zm[j]	- zm[i]) / Get_Cellsize();
			}
			else		// i => orthogonal
			{
				G		= (z		- zm[i]) / Get_Cellsize();
				H		= (zm[i]	- zm[j]) / Get_Cellsize();
			}

			if( H < 0.0 )
			{
				iAspect	= 0.0;
				iSlope	= G;
			}
			else if( H > G )
			{
				iAspect	= M_PI_045;
				iSlope	= (z - zm[i % 2 ? i : j]) / (sqrt(2.0) * Get_Cellsize());
			}
			else
			{
				iAspect	= atan(H / G);
				iSlope	= sqrt(G*G + H*H);
			}

			if( iSlope > Slope )
			{
				Aspect	= i * M_PI_045 + (i % 2 ? M_PI_045 - iAspect : iAspect);
				Slope	= iSlope;
			}
		}

		//---------------------------------------------
		if( Aspect < 0.0 )
		{
			Set_Parameters_NoData(x, y);
		}
		else
		{
			Set_Parameters(x, y, atan(Slope), Aspect);
		}
	}
}

//---------------------------------------------------------
// Least Squares or Best Fit Plane (Beasley & Huggins 1982, Costa-Cabral & Burgess 1994)
//
// Beasley, D.B. and Huggins, L.F. 1982:
//		ANSWERS: User抯 manual.
//		U.S. EPA-905/9-82-001, Chicago, IL. 54pp.
//
// Costa-Cabral, M., and Burges, S.J., 1994:
//		Digital Elevation Model Networks (DEMON): a model of flow over hillslopes for computation of contributing and dispersal areas
//		Water Resources Research, v. 30, no. 6, p. 1681-1692.
//
//---------------------------------------------------------
void CMorphometry::Do_LeastSquare(int x, int y)
{
	double	zm[9], a, b;

	if( Get_SubMatrix3x3(x, y, zm) )
	{
		a = ((zm[2] + 2 * zm[5] + zm[8]) - (zm[0] + 2 * zm[3] + zm[6])) / (8 * Get_Cellsize());
		b = ((zm[6] + 2 * zm[7] + zm[8]) - (zm[0] + 2 * zm[1] + zm[2])) / (8 * Get_Cellsize());

		if( a != 0.0 )
		{
			Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_180 + atan2(b, a));
		}
		else if( b > 0.0 )
		{
			Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_270);
		}
		else if( b < 0.0 )
		{
			Set_Parameters(x, y, atan( sqrt(a*a + b*b) ), M_PI_090);
		}
		else
		{
			Set_Parameters_NoData(x, y);
		}
	}
}

//---------------------------------------------------------
// Quadratic Function Approximation (Bauer, Rohdenburg & Bork, 1985)
//
// Bauer, J. / Rohdenburg, H. / Bork, H.-R., (1985):
//		'Ein Digitales Reliefmodell als Vorraussetzung fuer ein deterministisches Modell der Wasser- und Stoff-Fluesse',
//		Landschaftsgenese und Landschaftsoekologie, H.10, Parameteraufbereitung fuer deterministische Gebiets-Wassermodelle,
//		Grundlagenarbeiten zu Analyse von Agrar-Oekosystemen, (Eds.: Bork, H.-R. / Rohdenburg, H.), p.1-15
//
//---------------------------------------------------------
void CMorphometry::Do_FD_BRM(int x, int y)
{
	double	zm[9], D, E, F, G, H;

	if( Get_SubMatrix3x3(x, y, zm) )
	{
		D	= ( (zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8]) - 2 * (zm[1] + zm[4] + zm[7]) )	/ _DX_2;
		E	= ( (zm[0] + zm[6] + zm[1] + zm[7] + zm[2] + zm[8]) - 2 * (zm[3] + zm[4] + zm[5]) )	/ _DX_2;
        F	= (  zm[8] + zm[0] - zm[7] )														/ _4DX_2;
		G	= ( (zm[2] - zm[0]) + (zm[5] - zm[3]) + (zm[8]-zm[6]) )								/ _6DX;
		H	= ( (zm[6] - zm[0]) + (zm[7] - zm[1]) + (zm[8]-zm[2]) )								/ _6DX;

		Set_Parameters_Derive(x, y, D, E, F, G, H);
	}
}

//---------------------------------------------------------
// Quadratic Function Approximation (Heerdegen & Beran, 1984)
//
// Heerdegen, R.G. / Beran, M.A. (1982):
//		'Quantifying source areas through land surface curvature',
//		Journal of Hydrology, Vol.57
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Heerdegen(int x, int y)
{
	double	zm[9], D, E, F, G, H, a, b;

	//-----------------------------------------------------
	if( Get_SubMatrix3x3(x, y, zm) )
	{
		a	=   zm[0] + zm[2] + zm[3] + zm[5] + zm[6] + zm[8];
		b	=   zm[0] + zm[1] + zm[2] + zm[6] + zm[7] + zm[8];
		D	= (0.3 * a - 0.2 * b)								/ _DX_2;
		E	= (0.3 * b - 0.2 * a)								/ _DX_2;
		F	= ( zm[0] - zm[2]                 - zm[6] + zm[8])	/ _4DX_2;
		G	= (-zm[0] + zm[2] - zm[3] + zm[5] - zm[6] + zm[8])	/ _6DX;
        H	= (-zm[0] - zm[1] - zm[2] + zm[6] + zm[7] + zm[8])	/ _6DX;

		Set_Parameters_Derive(x, y, D, E, F, G, H);
	}
}

//---------------------------------------------------------
// Quadratic Function Approximation (Zevenbergen und Thorne, 1986)
//
// Zevenbergen, L.W. and C.R. Thorne. 1987:
//		Quantitative analysis of land surface topography
//		Earth Surface Processes and Landforms, 12: 47-56.
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Zevenbergen(int x, int y)
{
	double	zm[9], D, E, F, G, H;

	//-----------------------------------------------------
	if( Get_SubMatrix3x3(x, y, zm) )
	{
		D	= ((zm[3] + zm[5]) / 2.0 - zm[4])	/ _DX_2;
		E	= ((zm[1] + zm[7]) / 2.0 - zm[4])	/ _DX_2;
		F	=  (zm[0] - zm[2] - zm[6] + zm[8])	/ _4DX_2;
		G	=  (zm[5] - zm[3])					/ _2DX;
        H	=  (zm[7] - zm[1])					/ _2DX;

		Set_Parameters_Derive(x, y, D, E, F, G, H);
	}
}

//---------------------------------------------------------
// Cubic Function Approximation (Haralick, 1991)
//
// R.M. Haralick (1983):
//		'Ridge and Valley Detection on digital images',
//		Computer Vision, Graphics and Image Processing, Vol.22, No.1, p.28-38
//
//---------------------------------------------------------
void CMorphometry::Do_FD_Haralick(int x, int y)
{
	//-----------------------------------------------------
	// Finite Differenzen Methode Matrizen...

	const int 	Mtrx[][5][5]	= {
	{	{ 31,- 5,-17,- 5, 31}, {-44,-62,-68,-62,-44}, {  0,  0,  0,  0,  0}, { 44, 62, 68, 62, 44}, {-31,  5, 17,  5,-31}	},
	{	{ 31,-44,  0, 44,-31}, {- 5,-62,  0, 62,  5}, {-17,-68,  0, 68, 17}, {- 5,-62,  0, 62,  5}, { 31,-44,  0, 44,-31}	},
	{	{  2,  2,  2,  2,  2}, {- 1,- 1,- 1,- 1,- 1}, {- 2,- 2,- 2,- 2,- 2}, {- 1,- 1,- 1,- 1,- 1}, {  2,  2,  2,  2,  2}	},
	{	{  4,  2,  0,- 2,- 4}, {  2,  1,  0,- 1,- 2}, {  0,  0,  0,  0,  0}, {- 2,- 1,  0,  1,  2}, {- 4,- 2,  0,  2,  4}	},
	{	{  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}, {  2,- 1,- 2,- 1,  2}	},	};

	const int	QMtrx[]			= { 4200, 4200, 700, 1000, 700 };

	//-----------------------------------------------------
	int		i, ix, iy, n;
	double	Sum, zm[25], k[5];

	//-----------------------------------------------------
	if( Get_SubMatrix5x5(x, y, zm) )
	{
		for(i=0; i<5; i++)
		{
			for(n=0, Sum=0.0, iy=0; iy<5; iy++)
			{
				for(ix=0; ix<5; ix++, n++)
				{
					Sum	+= zm[n] * Mtrx[i][ix][iy];
				}
			}

			k[i]	= Sum / QMtrx[i];
		}

		Set_Parameters_Derive(x, y, k[4], k[2], k[3], k[1], k[0]);
    }
}