grid_operation.cpp

这是一个GPS相关的程序

CSG_Grid & CSG_Grid::operator /=	(const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
}

CSG_Grid & CSG_Grid::operator /=	(double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Division) );
}

CSG_Grid & CSG_Grid::Divide			(const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
}

CSG_Grid & CSG_Grid::Divide			(double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Division) );
}


///////////////////////////////////////////////////////////
//														 //
//					Operatoren							 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSG_Grid & CSG_Grid::_Operation_Arithmetic(const CSG_Grid &Grid, TSG_Grid_Operation Operation)
{
	if( is_Intersecting(Grid.Get_Extent()) )
	{
		int						x, y;
		double					xWorld, yWorld, Value;
		TSG_Grid_Interpolation	Interpolation;

		Interpolation	=	Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_XMin() - Grid.Get_XMin(), Get_Cellsize()) == 0.0
						&&	Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_YMin() - Grid.Get_YMin(), Get_Cellsize()) == 0.0
						?	GRID_INTERPOLATION_NearestNeighbour
						:	GRID_INTERPOLATION_BSpline;

		for(y=0, yWorld=Get_YMin(); y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++, yWorld+=Get_Cellsize())
		{
			for(x=0, xWorld=Get_XMin(); x<Get_NX(); x++, xWorld+=Get_Cellsize())
			{
				if( !Grid.Get_Value(xWorld, yWorld, Value, Interpolation, true) )
				{
					Set_NoData(x, y);
				}
				else switch( Operation )
				{
				case GRID_OPERATION_Addition:
					Add_Value(x, y,  Value);
					break;

				case GRID_OPERATION_Subtraction:
					Add_Value(x, y, -Value);
					break;

				case GRID_OPERATION_Multiplication:
					Mul_Value(x, y,  Value);
					break;

				case GRID_OPERATION_Division:
					if( Value != 0.0 )
					{
						Mul_Value(x, y, 1.0 / Value);
					}
					else
					{
						Set_NoData(x, y);
					}
					break;
				}
			}
		}

		SG_UI_Process_Set_Ready();

		//-------------------------------------------------
		switch( Operation )
		{
		case GRID_OPERATION_Addition:
			Get_History().Add_Entry(LNG("[HST] Grid addition")		, Grid.Get_Name());
			break;

		case GRID_OPERATION_Subtraction:
			Get_History().Add_Entry(LNG("[HST] Grid subtraction")	, Grid.Get_Name());
			break;

		case GRID_OPERATION_Multiplication:
			Get_History().Add_Entry(LNG("[HST] Grid multiplication"), Grid.Get_Name());
			break;

		case GRID_OPERATION_Division:
			Get_History().Add_Entry(LNG("[HST] Grid division")		, Grid.Get_Name());
			break;
		}

		Get_History().Assign(Grid.m_History, true);
	}

	return( *this );
}

//---------------------------------------------------------
CSG_Grid & CSG_Grid::_Operation_Arithmetic(double Value, TSG_Grid_Operation Operation)
{
	if( Operation == GRID_OPERATION_Division && Value == 0.0 )
	{
		return( *this );
	}

	//-----------------------------------------------------
	switch( Operation )
	{
	case GRID_OPERATION_Addition:
		Get_History().Add_Entry(LNG("[HST] Value addition")			, CSG_String::Format(SG_T("%f"), Value));
		break;

	case GRID_OPERATION_Subtraction:
		Get_History().Add_Entry(LNG("[HST] Value subtraction")		, CSG_String::Format(SG_T("%f"), Value));
		Value	= -Value;
		break;

	case GRID_OPERATION_Multiplication:
		Get_History().Add_Entry(LNG("[HST] Value multiplication")	, CSG_String::Format(SG_T("%f"), Value));
		break;

	case GRID_OPERATION_Division:
		Get_History().Add_Entry(LNG("[HST] Value division")			, CSG_String::Format(SG_T("%f"), Value));
		break;
	}

	//-----------------------------------------------------
	for(int y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
	{
		for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				switch( Operation )
				{
				case GRID_OPERATION_Addition:
				case GRID_OPERATION_Subtraction:
					Add_Value(x, y, Value);
					break;

				case GRID_OPERATION_Multiplication:
				case GRID_OPERATION_Division:
					Mul_Value(x, y, Value);
					break;
				}
			}
		}
	}

	SG_UI_Process_Set_Ready();

	return( *this );
}


///////////////////////////////////////////////////////////
//														 //
//					Grid-Operations - A					 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CSG_Grid::Invert(void)
{
	int		x, y;
	double	zMin, zMax;

	if( is_Valid() && Get_ZRange() > 0.0 )
	{
		zMin	= Get_ZMin();
		zMax	= Get_ZMax();

		for(y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
		{
			for(x=0; x<Get_NX(); x++)
			{
				if( !is_NoData(x, y) )
				{
					Set_Value(x, y, zMax - (asDouble(x, y) - zMin));
				}
			}
		}

		SG_UI_Process_Set_Ready();

		Get_History().Add_Entry(LNG("[HST] Grid operation"), LNG("[HST] Inversion"));
	}
}

//---------------------------------------------------------
void CSG_Grid::Flip(void)
{
	int		x, yA, yB;
	double	*Line, d;

	if( is_Valid() )
	{
		Line	= (double *)SG_Malloc(Get_NX() * sizeof(double));

		for(yA=0, yB=Get_NY()-1; yA<yB && SG_UI_Process_Set_Progress(2 * yA, Get_NY()); yA++, yB--)
		{
			for(x=0; x<Get_NX(); x++)
			{
				Line[x]	= asDouble(x, yA);
			}

			for(x=0; x<Get_NX(); x++)
			{
				d		= Line[x];
				Line[x]	= asDouble(x, yB);
				Set_Value(x, yB, d);
			}

			for(x=0; x<Get_NX(); x++)
			{
				Set_Value(x, yA, Line[x]);
			}
		}

		SG_UI_Process_Set_Ready();

		SG_Free(Line);

		Get_History().Add_Entry(LNG("[HST] Grid operation"), LNG("[HST] Vertically mirrored"));
	}
}

//---------------------------------------------------------
void CSG_Grid::Mirror(void)
{
	int		xA, xB, y;
	double	d;

	if( is_Valid() )
	{
		for(y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
		{
			for(xA=0, xB=Get_NX()-1; xA<xB; xA++, xB--)
			{
				d			=    asDouble(xA, y);
				Set_Value(xA, y, asDouble(xB, y));
				Set_Value(xB, y, d);
			}
		}

		SG_UI_Process_Set_Ready();

		Get_History().Add_Entry(LNG("[HST] Grid operation"), LNG("[HST] Horizontally mirrored"));
	}
}


///////////////////////////////////////////////////////////
//														 //
//					Grid-Operations - B					 //
//														 //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
void CSG_Grid::Normalise(void)
{
	int		x, y;
	double	d;

	if( is_Valid() )
	{
		Update_Statistics();

		if( m_Variance > 0.0 )
		{
			d	= sqrt(m_Variance);

			for(y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
			{
				for(x=0; x<Get_NX(); x++)
				{
					if( !is_NoData(x, y) )
					{
						Set_Value(x, y, (asDouble(x, y) - m_ArithMean) / d );
					}
				}
			}

			SG_UI_Process_Set_Ready();

			Get_History().Add_Entry(LNG("[HST] Grid normalisation"), CSG_String::Format(SG_T("%f / %f"), m_ArithMean, m_Variance));
		}
	}
}

//---------------------------------------------------------
void CSG_Grid::DeNormalise(double ArithMean, double Variance)
{
	int		x, y;

	if( is_Valid() )
	{
		Variance	= sqrt(Variance);

		for(y=0; y<Get_NY() && SG_UI_Process_Set_Progress(y, Get_NY()); y++)
		{
			for(x=0; x<Get_NX(); x++)
			{
				if( !is_NoData(x, y) )
				{
					Set_Value(x, y, Variance * asDouble(x, y) + ArithMean);
				}
			}
		}

		SG_UI_Process_Set_Ready();

		Get_History().Add_Entry(LNG("[HST] Grid denormalisation"), CSG_String::Format(SG_T("%f / %f"), ArithMean, Variance));
	}
}


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

//---------------------------------------------------------
int CSG_Grid::Get_Gradient_NeighborDir(int x, int y, bool bMustBeLower)	const
{
	int		i, ix, iy, Direction;
	double	z, dz, dzMax;

	Direction	= -1;

	if( is_InGrid(x, y) )
	{
		z		= asDouble(x, y);
		dzMax	= 0.0;

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

			if( !is_InGrid(ix, iy) )
			{
				return( i );
			}
			else
			{
				dz	= (z - asDouble(ix, iy)) / m_System.Get_Length(i);

				if( (bMustBeLower && dz > 0.0) || !bMustBeLower )
				{
					if( Direction < 0 || (dz > dzMax) )
					{
						Direction	= i;
						dzMax		= dz;
					}
				}
			}
		}
	}

	return( Direction );
}

//---------------------------------------------------------
bool CSG_Grid::Get_Gradient(int x, int y, double &Decline, double &Azimuth) const
{
	int		i, ix, iy, iDir;
	double	z, zm[4], G, H;

	if( is_InGrid(x, y) )
	{
		z		= asDouble(x, y);

		for(i=0, iDir=0; i<4; i++, iDir+=2)
		{
			ix		= m_System.Get_xTo(iDir, x);
			iy		= m_System.Get_yTo(iDir, y);

			if( is_InGrid(ix, iy) )
			{
				zm[i]	= asDouble(ix, iy) - z;
			}
			else
			{
				ix		= m_System.Get_xFrom(iDir, x);
				iy		= m_System.Get_yFrom(iDir, y);

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

		G		= (zm[0] - zm[2]) / (2.0 * Get_Cellsize());
        H		= (zm[1] - zm[3]) / (2.0 * Get_Cellsize());

		Decline	= atan(sqrt(G*G + H*H));

		if( G != 0.0 )
			Azimuth	= M_PI_180 + atan2(H, G);
		else
			Azimuth	= H > 0.0 ? M_PI_270 : (H < 0.0 ? M_PI_090 : -1.0);

		return( true );
	}

	Decline	= 0.0;
	Azimuth	= -1.0;

	return( false );
}


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

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