mrvbf.cpp

这是一个GPS相关的程序

		for(y=0; y<pDEM->Get_NY() && Set_Progress(y, pDEM->Get_NY()); y++)
		{
			for(x=0; x<pDEM->Get_NX(); x++)
			{
				for(iy=-Radius, jy=y-Radius, ky=0, s=d=0.0; iy<=Radius; iy++, jy++, ky++)
				{
					for(ix=-Radius, jx=x-Radius, kx=0; ix<=Radius; ix++, jx++, kx++)
					{
						if( pDEM->is_InGrid(jx, jy) )
						{
							s	+= Kernel.asDouble(kx, ky) * pDEM->asDouble(jx, jy);
							d	+= Kernel.asDouble(kx, ky);
						}
					}
				}

				if( d > 0.0 )
				{
					pSmoothed->Set_Value(x, y, s / d);
				}
				else
				{
					pSmoothed->Set_NoData(x, y);
				}
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CMRVBF::Get_Values(CSG_Grid *pDEM, CSG_Grid *pSlopes, CSG_Grid *pPercentiles, double Resolution)
{
	if( pDEM && pDEM->is_Valid() && pSlopes && pPercentiles )
	{
		int			nx, ny;
		CSG_Grid	Smoothed;

		Get_Smoothed(pDEM, &Smoothed, 5, 3.0);
		Get_Slopes(&Smoothed, pSlopes);

		nx	= 2 + (int)(pDEM->Get_XRange() / Resolution);
		ny	= 2 + (int)(pDEM->Get_YRange() / Resolution);

		pDEM->Create(GRID_TYPE_Float, nx, ny, Resolution, pDEM->Get_XMin(), pDEM->Get_YMin());
		pDEM->Assign(&Smoothed, GRID_INTERPOLATION_NearestNeighbour);

		Get_Percentiles(pDEM, pPercentiles, 6);

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CMRVBF::Get_Percentile(CSG_Grid *pDEM, int x, int y, double &Percentile)
{
	if( pDEM && pDEM->is_Valid() && pDEM->is_InGrid(x, y, true) )
	{
		int		iRadius, iPoint, nPoints, nLower, ix, iy;
		double	z	= pDEM->asDouble(x, y);

		for(iRadius=0, nPoints=0, nLower=0; iRadius<m_Radius.Get_Maximum(); iRadius++)
		{
			for(iPoint=0; iPoint<m_Radius.Get_nPoints(iRadius); iPoint++)
			{
				m_Radius.Get_Point(iRadius, iPoint, ix, iy);

				ix	+= x;
				iy	+= y;

				if( pDEM->is_InGrid(ix, iy) )
				{
					nPoints++;

					if( pDEM->asDouble(ix, iy) < z )
					{
						nLower++;
					}
				}
			}
		}

		if( nPoints > 1 )
		{
			Percentile	= (double)nLower / (double)(nPoints - 1.0);

			return( true );
		}
	}

	return( false );
}

//---------------------------------------------------------
bool CMRVBF::Get_Percentiles(CSG_Grid *pDEM, CSG_Grid *pPercentiles, int Radius)
{
	if( pDEM && pDEM->is_Valid() )
	{
		pPercentiles->Create(pDEM->Get_System(), GRID_TYPE_Float);

		m_Radius.Create(Radius);

		for(int y=0; y<pDEM->Get_NY() && Set_Progress(y, pDEM->Get_NY()); y++)
		{
			for(int x=0; x<pDEM->Get_NX(); x++)
			{
				double	Percentile;

				if( !Get_Percentile(pDEM, x, y, Percentile) )
				{
					pPercentiles->Set_NoData(x, y);
				}
				else
				{
					pPercentiles->Set_Value (x, y, Percentile);
				}
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CMRVBF::Get_Slopes(CSG_Grid *pDEM, CSG_Grid *pSlopes)
{
	if( pDEM && pDEM->is_Valid() )
	{
		pSlopes->Create(pDEM->Get_System(), GRID_TYPE_Float);

		for(int y=0; y<pDEM->Get_NY() && Set_Progress(y, pDEM->Get_NY()); y++)
		{
			for(int x=0; x<pDEM->Get_NX(); x++)
			{
				double	Slope, Aspect;

				if( !pDEM->Get_Gradient(x, y, Slope, Aspect) )
				{
					pSlopes->Set_NoData(x, y);
				}
				else
				{
					pSlopes->Set_Value (x, y, 100.0 * tan(Slope));
				}
			}
		}

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CMRVBF::Get_Flatness(CSG_Grid *pSlopes, CSG_Grid *pPercentiles, CSG_Grid *pCF, CSG_Grid *pVF, CSG_Grid *pRF, double T_Slope)
{
//	const int	Interpolation	= GRID_INTERPOLATION_Bilinear;
	const int	Interpolation	= GRID_INTERPOLATION_BSpline;

	if( pSlopes && pSlopes->is_Valid() && pPercentiles && pPercentiles->is_Valid() )
	{
		int		x, y;
		double	xp, yp, Slope, Percentile, cf, vf, rf;

		for(y=0, yp=Get_YMin(); y<Get_NY() && Set_Progress(y); y++, yp+=Get_Cellsize())
		{
			for(x=0, xp=Get_XMin(); x<Get_NX(); x++, xp+=Get_Cellsize())
			{
				if( pSlopes		->Get_Value(xp, yp, Slope     , Interpolation)
				&&	pPercentiles->Get_Value(xp, yp, Percentile, Interpolation) )
				{
					cf	= pCF->asDouble(x, y) * Get_Transformation(Slope, T_Slope, m_P_Slope);
					vf	= cf * Get_Transformation(      Percentile, m_T_Pctl_V, m_P_Pctl);
					rf	= cf * Get_Transformation(1.0 - Percentile, m_T_Pctl_R, m_P_Pctl);

					pCF->Set_Value	(x, y, cf);
					pVF->Set_Value	(x, y, 1.0 - Get_Transformation(vf, 0.3, 4.0));
					pRF->Set_Value	(x, y, 1.0 - Get_Transformation(rf, 0.3, 4.0));
				}
				else
				{
					pVF->Set_NoData	(x, y);
					pRF->Set_NoData	(x, y);
				}
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CMRVBF::Get_MRVBF(int Level, CSG_Grid *pMRVBF, CSG_Grid *pVF, CSG_Grid *pMRRTF, CSG_Grid *pRF)
{
	if( pMRVBF && pVF && pMRRTF && pRF )
	{
		double	d, w, t, p;

		t	= 0.4;
		p	= log((Level - 0.5) / 0.1) / log(1.5);

		for(int y=0; y<Get_NY() && Set_Progress(y); y++)
		{
			for(int x=0; x<Get_NX(); x++)
			{
				if( !pMRVBF->is_NoData(x, y) && !pVF->is_NoData(x, y) )
				{
					d	= pVF->asDouble(x, y);
					w	= 1.0 - Get_Transformation(d, t, p);
					pMRVBF->Set_Value(x, y, w * (Level - 1 + d) + (1.0 - w) * pMRVBF->asDouble(x, y));
				}

				if( !pMRRTF->is_NoData(x, y) && !pRF->is_NoData(x, y) )
				{
					d	= pRF->asDouble(x, y);
					w	= 1.0 - Get_Transformation(d, t, p);
					pMRRTF->Set_Value(x, y, w * (Level - 1 + d) + (1.0 - w) * pMRRTF->asDouble(x, y));
				}
			}
		}

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CMRVBF::Get_Classified(CSG_Grid *pMRF)
{
	if( pMRF && pMRF->is_Valid() )
	{
		for(int y=0; y<Get_NY() && Set_Progress(y); y++)
		{
			for(int x=0; x<Get_NX(); x++)
			{
				if( !pMRF->is_NoData(x, y) )
				{
					double	d	= pMRF->asDouble(x, y);

					if		( d < 0.5 )
						pMRF->Set_Value(x, y, 0.0);
					else if	( d < 1.5 )
						pMRF->Set_Value(x, y, 1.0);
					else if	( d < 2.5 )
						pMRF->Set_Value(x, y, 2.0);
					else if	( d < 3.5 )
						pMRF->Set_Value(x, y, 3.0);
					else if	( d < 4.5 )
						pMRF->Set_Value(x, y, 4.0);
					else if	( d < 5.5 )
						pMRF->Set_Value(x, y, 5.0);
					else
						pMRF->Set_Value(x, y, 6.0);
				}
			}
		}

		return( true );
	}

	return( false );
}


///////////////////////////////////////////////////////////
//														 //
//														 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------