
///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//           Application Programming Interface           //
//                                                       //
//                  Library: SAGA_API                    //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//                  grid_operation.cpp                   //
//                                                       //
//          Copyright (C) 2005 by Olaf Conrad            //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'.                              //
//                                                       //
// This library is free software; you can redistribute   //
// it and/or modify it under the terms of the GNU Lesser //
// General Public License as published by the Free       //
// Software Foundation, either version 2.1 of the        //
// License, or (at your option) any later version.       //
//                                                       //
// This library is distributed in the hope that it will  //
// be useful, but WITHOUT ANY WARRANTY; without even the //
// implied warranty of MERCHANTABILITY or FITNESS FOR A  //
// PARTICULAR PURPOSE. See the GNU Lesser General Public //
// License for more details.                             //
//                                                       //
// You should have received a copy of the GNU Lesser     //
// General Public License along with this program; if    //
// not, see <http://www.gnu.org/licenses/>.              //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
//    e-mail:     oconrad@saga-gis.org                   //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include <memory.h>

#include "grid.h"


///////////////////////////////////////////////////////////
//                                                       //
//                     Assignments                       //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
bool CSG_Grid::Assign_NoData(void)
{
	double Value = Get_NoData_Value();

	#pragma omp parallel for
	for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
	{
		Set_Value(x, y, Value, false);
	}

	return( true );
}

//---------------------------------------------------------
bool CSG_Grid::Assign(double Value)
{
	if( !is_Valid() )
	{
		return( false );
	}

	if( Value == 0. && !is_Cached() )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++)
		{
			memset(m_Values[y], 0, Get_nLineBytes());
		}

		Set_Modified();
	}
	else
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			Set_Value(x, y, Value);
		}
	}

	//-----------------------------------------------------
	Get_History().Destroy(); m_Statistics.Invalidate();

	return( true );
}


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

//---------------------------------------------------------
bool CSG_Grid::Assign(CSG_Data_Object *pObject, bool bProgress)
{
	return( pObject && Assign(pObject->asGrid(), CSG_Grid_Resampling::Undefined) );
}

bool CSG_Grid::Assign(CSG_Grid *pGrid, CSG_Grid_Resampling Method, bool bProgress)
{
	if(	!is_Valid() || !pGrid || !pGrid->is_Valid() || is_Intersecting(pGrid->Get_Extent()) == INTERSECTION_None )
	{
		return( false );
	}

	bool bResult = false;

	//---------------------------------------------------------
	if( m_System == pGrid->m_System )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( pGrid->is_NoData(x, y) )
			{
				Set_NoData(x, y);
			}
			else
			{
				Set_Value(x, y, pGrid->asDouble(x, y));
			}
		}

		bResult = true; bProgress = false;
	}

	//---------------------------------------------------------
	else if(  Get_Cellsize() == pGrid->Get_Cellsize() // No-Scaling...
	&& fmod(Get_XMin() - pGrid->Get_XMin(), Get_Cellsize()) == 0.
	&& fmod(Get_YMin() - pGrid->Get_YMin(), Get_Cellsize()) == 0.	)
	{
		bResult = _Assign_Interpolated(pGrid, CSG_Grid_Resampling::NearestNeighbour, bProgress);
	}

	//---------------------------------------------------------
	else switch( Method )
	{
	case CSG_Grid_Resampling::NearestNeighbour:
	case CSG_Grid_Resampling::Bilinear        :
	case CSG_Grid_Resampling::Bicubic_1       :
	case CSG_Grid_Resampling::Bicubic_2       : bResult = _Assign_Interpolated(pGrid, Method                                   , bProgress); break;

	case CSG_Grid_Resampling::Mean_Nodes      :
	case CSG_Grid_Resampling::Mean_Cells      : bResult = _Assign_MeanValue   (pGrid, Method != CSG_Grid_Resampling::Mean_Nodes, bProgress); break;

	case CSG_Grid_Resampling::Minimum         :
	case CSG_Grid_Resampling::Maximum         : bResult = _Assign_ExtremeValue(pGrid, Method == CSG_Grid_Resampling::Maximum   , bProgress); break;

	case CSG_Grid_Resampling::Majority        : bResult = _Assign_Majority    (pGrid                                           , bProgress); break;

	default                                   : bResult = Get_Cellsize() < pGrid->Get_Cellsize()
		? _Assign_Interpolated(pGrid, CSG_Grid_Resampling::Bicubic_2, bProgress)  // Down-Scaling...
		: _Assign_MeanValue   (pGrid, true                          , bProgress); // Up-Scaling...
		break;
	}

	//---------------------------------------------------------
	if( bResult )
	{
		Set_Unit(pGrid->Get_Unit());

		if( pGrid->Get_Projection().is_Okay() )
		{
			Get_Projection() = pGrid->Get_Projection();
		}

		Get_History() = pGrid->Get_History();
	}

	//---------------------------------------------------------
	if( bProgress )
	{
		SG_UI_Process_Set_Ready();
	}

	return( bResult );
}

//---------------------------------------------------------
bool CSG_Grid::_Assign_Interpolated(CSG_Grid *pGrid, CSG_Grid_Resampling Interpolation, bool bProgress)
{
	for(int y=0; y<Get_NY() && (!bProgress || SG_UI_Process_Set_Progress(y, Get_NY())); y++)
	{
		double py = Get_YMin() + y * Get_Cellsize();

		#pragma omp parallel for
		for(int x=0; x<Get_NX(); x++)
		{
			double px = Get_XMin() + x * Get_Cellsize(), z;

			if( pGrid->Get_Value(px, py, z, Interpolation) )
			{
				Set_Value(x, y, z);
			}
			else
			{
				Set_NoData(x, y);
			}
		}
	}

	return( true );
}

//---------------------------------------------------------
bool CSG_Grid::_Assign_ExtremeValue(CSG_Grid *pGrid, bool bMaximum, bool bProgress)
{
	if( Get_Cellsize() < pGrid->Get_Cellsize() )
	{
		return( false );
	}

	//-----------------------------------------------------
	Assign_NoData();

	double d = pGrid->Get_Cellsize() / Get_Cellsize();

	double ax = 0.5 + (pGrid->Get_XMin() - Get_XMin()) / Get_Cellsize();
	double py = 0.5 + (pGrid->Get_YMin() - Get_YMin()) / Get_Cellsize();

	for(int y=0; y<pGrid->Get_NY() && (!bProgress || SG_UI_Process_Set_Progress(y, pGrid->Get_NY())); y++, py+=d)
	{
		int iy = (int)floor(py);

		if( iy >= 0 && iy < Get_NY() )
		{
			#pragma omp parallel for
			for(int x=0; x<pGrid->Get_NX(); x++)
			{
				if( !pGrid->is_NoData(x, y) )
				{
					int ix = (int)floor(ax + x * d);
					
					if( ix >= 0 && ix < Get_NX() )
					{
						double z = pGrid->asDouble(x, y);

						if( is_NoData(ix, iy)
						|| (bMaximum ==  true && z > asDouble(ix, iy))
						|| (bMaximum == false && z < asDouble(ix, iy)) )
						{
							Set_Value(ix, iy, z);
						}
					}
				}
			}
		}
	}

	//-----------------------------------------------------
	return( true );
}

//---------------------------------------------------------
bool CSG_Grid::_Assign_MeanValue(CSG_Grid *pGrid, bool bAreaProportional, bool bProgress)
{
	if( Get_Cellsize() < pGrid->Get_Cellsize() )
	{
		return( false );
	}

	//-----------------------------------------------------
	double d = Get_Cellsize() / pGrid->Get_Cellsize();

	double ox = (Get_XMin(true) - pGrid->Get_XMin()) / pGrid->Get_Cellsize();
	double py = (Get_YMin(true) - pGrid->Get_YMin()) / pGrid->Get_Cellsize();

	for(int y=0; y<Get_NY() && (!bProgress || SG_UI_Process_Set_Progress(y, Get_NY())); y++, py+=d)
	{
		int ay = (int)(bAreaProportional ? floor(py    ) : ceil (py    ));
		int by = (int)(bAreaProportional ? ceil (py + d) : floor(py + d));

		//-------------------------------------------------
		#pragma omp parallel for
		for(int x=0; x<Get_NX(); x++)
		{
			double px = ox + x * d;

			int ax = (int)(bAreaProportional ? floor(px    ) : ceil (px    ));
			int bx = (int)(bAreaProportional ? ceil (px + d) : floor(px + d));

			CSG_Rect rMean(px, py, px + d, py + d); CSG_Simple_Statistics s;

			for(int iy=ay; iy<=by; iy++)
			{
				if( iy >= 0 && iy < pGrid->Get_NY() )
				{
					for(int ix=ax; ix<=bx; ix++)
					{
						if( ix >= 0 && ix < pGrid->Get_NX() && !pGrid->is_NoData(ix, iy) )
						{
							if( bAreaProportional )
							{
								CSG_Rect r(ix - 0.5, iy - 0.5, ix + 0.5, iy + 0.5);

								if( r.Intersect(rMean) )
								{
									s.Add_Value(pGrid->asDouble(ix, iy), r.Get_Area());
								}
							}
							else
							{
								s.Add_Value(pGrid->asDouble(ix, iy));
							}
						}
					}
				}
			}

			//---------------------------------------------
			if( s.Get_Count() > 0 )
			{
				Set_Value(x, y, s.Get_Mean());
			}
			else
			{
				Set_NoData(x, y);
			}
		}
	}

	//-----------------------------------------------------
	return( true );
}


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

//---------------------------------------------------------
bool CSG_Grid::_Assign_Majority(CSG_Grid *pGrid, bool bProgress)
{
	if( Get_Cellsize() < pGrid->Get_Cellsize() )
	{
		return( false );
	}

	//-----------------------------------------------------
	Set_NoData_Value(pGrid->Get_NoData_Value());

	Assign_NoData();

	//-----------------------------------------------------
	int ay, by = (int)(1. + (((0 - 0.5) * Get_Cellsize() + Get_YMin()) - pGrid->Get_YMin()) / pGrid->Get_Cellsize());

	for(int y=0; y<Get_NY() && (!bProgress || SG_UI_Process_Set_Progress(y, Get_NY())); y++)
	{
		ay = by; by = (int)(1. + (((y + 0.5) * Get_Cellsize() + Get_YMin()) - pGrid->Get_YMin()) / pGrid->Get_Cellsize());

		if( ay < pGrid->Get_NY() && by > 0 )
		{
			if( ay < 0 ) { ay = 0; } if( by > pGrid->Get_NY() ) { by = pGrid->Get_NY(); }

			int ax, bx = (int)(1. + (((0 - 0.5) * Get_Cellsize() + Get_XMin()) - pGrid->Get_XMin()) / pGrid->Get_Cellsize());

			for(int x=0; x<Get_NX(); x++)
			{
				ax = bx; bx = (int)(1. + (((x + 0.5) * Get_Cellsize() + Get_XMin()) - pGrid->Get_XMin()) / pGrid->Get_Cellsize());

				if( ax < pGrid->Get_NX() && bx > 0 )
				{
					CSG_Unique_Number_Statistics s;

					if( ax < 0 ) { ax = 0; } if( bx > pGrid->Get_NX() ) { bx = pGrid->Get_NX(); }

					for(int iy=ay; iy<by; iy++) for(int ix=ax; ix<bx; ix++)
					{
						if( !pGrid->is_NoData(ix, iy) )
						{
							s += pGrid->asDouble(ix, iy);
						}
					}

					int n; double z;

					if( s.Get_Majority(z, n) )//&& n > 1 )
					{
						Set_Value(x, y, z);
					}
				}
			}
		}
	}

	//-----------------------------------------------------
	return( true );
}


///////////////////////////////////////////////////////////
//                                                       //
//                      Operators                        //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSG_Grid & CSG_Grid::operator =  (const CSG_Grid &Grid)
{
	Assign(Grid.asGrid(), CSG_Grid_Resampling::Undefined, false);

	return( *this );
}

CSG_Grid & CSG_Grid::operator =  (double Value)
{
	Assign(Value);

	return( *this );
}

//---------------------------------------------------------
CSG_Grid CSG_Grid::operator   +  (const CSG_Grid &Grid) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
}

CSG_Grid CSG_Grid::operator   +  (double Value) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
}

CSG_Grid & CSG_Grid::operator += (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
}

CSG_Grid & CSG_Grid::operator += (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
}

CSG_Grid & CSG_Grid::Add         (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Addition) );
}

CSG_Grid & CSG_Grid::Add         (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Addition) );
}

//---------------------------------------------------------
CSG_Grid CSG_Grid::operator   -  (const CSG_Grid &Grid) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
}

CSG_Grid CSG_Grid::operator   -  (double Value) const
{
	CSG_Grid	g(*this);

	return( g._Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
}

CSG_Grid & CSG_Grid::operator -= (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
}

CSG_Grid & CSG_Grid::operator -= (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
}

CSG_Grid & CSG_Grid::Subtract    (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Subtraction) );
}

CSG_Grid & CSG_Grid::Subtract    (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Subtraction) );
}

//---------------------------------------------------------
CSG_Grid CSG_Grid::operator   *  (const CSG_Grid &Grid) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
}

CSG_Grid CSG_Grid::operator   *  (double Value) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
}

CSG_Grid & CSG_Grid::operator *= (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
}

CSG_Grid & CSG_Grid::operator *= (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
}

CSG_Grid & CSG_Grid::Multiply    (const CSG_Grid &Grid)
{
	return( _Operation_Arithmetic(Grid, GRID_OPERATION_Multiplication) );
}

CSG_Grid & CSG_Grid::Multiply    (double Value)
{
	return( _Operation_Arithmetic(Value, GRID_OPERATION_Multiplication) );
}

//---------------------------------------------------------
CSG_Grid CSG_Grid::operator   /  (const CSG_Grid &Grid) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Grid, GRID_OPERATION_Division) );
}

CSG_Grid CSG_Grid::operator  /   (double Value) const
{
	CSG_Grid g(*this);

	return( g._Operation_Arithmetic(Value, GRID_OPERATION_Division) );
}

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) );
}


///////////////////////////////////////////////////////////
//                                                       //
//                      Operations                       //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
CSG_Grid & CSG_Grid::_Operation_Arithmetic(const CSG_Grid &Grid, TSG_Grid_Operation Operation)
{
	if( is_Intersecting(Grid.Get_Extent()) )
	{
		CSG_Grid_Resampling	Interpolation	=
			Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_XMin() - Grid.Get_XMin(), Get_Cellsize()) == 0.
		&&	Get_Cellsize() == Grid.Get_Cellsize() && fmod(Get_YMin() - Grid.Get_YMin(), Get_Cellsize()) == 0.
		?	CSG_Grid_Resampling::NearestNeighbour
		:	CSG_Grid_Resampling::Bicubic_2;

		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++)
		{
			double yWorld = Get_YMin() + y * Get_Cellsize();

			for(int x=0; x<Get_NX(); x++)
			{
				if( !is_NoData(x, y) )
				{
					double xWorld = Get_XMin() + x * Get_Cellsize(), Value;

					if( Grid.Get_Value(xWorld, yWorld, Value, Interpolation) )
					{
						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. )
							{
								Mul_Value(x, y, 1. / Value);
							}
							else
							{
								Set_NoData(x, y);
							}
							break;
						}
					}
				}
			}
		}
	}

	return( *this );
}

//---------------------------------------------------------
CSG_Grid & CSG_Grid::_Operation_Arithmetic(double Value, TSG_Grid_Operation Operation)
{
	switch( Operation )
	{
	case GRID_OPERATION_Addition      : if( Value == 0. ) { return( *this ); }                     break;
	case GRID_OPERATION_Subtraction   : if( Value == 0. ) { return( *this ); } Value =     -Value; break;
	case GRID_OPERATION_Multiplication: if( Value == 1. ) { return( *this ); }                     break;
	case GRID_OPERATION_Division      : if( Value == 0. ) { return( *this ); } Value = 1. / Value; break;
	}

	//-----------------------------------------------------
	#pragma omp parallel for
	for(int y=0; 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;
			}
		}
	}

	//-----------------------------------------------------
	return( *this );
}


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

//---------------------------------------------------------
void CSG_Grid::Invert(void)
{
	if( is_Valid() && Get_Range() > 0. )
	{
		double Min = Get_Min(), Max = Get_Max();

		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				Set_Value(x, y, Max - (asDouble(x, y) - Min));
			}
		}
	}
}

//---------------------------------------------------------
void CSG_Grid::Flip(void)
{
	if( is_Valid() )
	{
		#pragma omp parallel for
		for(int x=0; x<Get_NX(); x++) for(int yA=0, yB=Get_NY()-1; yA<yB; yA++, yB--)
		{
			double v = asDouble(x, yA); Set_Value(x, yA, asDouble(x, yB)); Set_Value(x, yB, v);
		}
	}
}

//---------------------------------------------------------
void CSG_Grid::Mirror(void)
{
	if( is_Valid() )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int xA=0, xB=Get_NX()-1; xA<xB; xA++, xB--)
		{
			double v = asDouble(xA, y); Set_Value(xA, y, asDouble(xB, y)); Set_Value(xB, y, v);
		}
	}
}


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

//---------------------------------------------------------
bool CSG_Grid::Normalise(void)
{
	if( is_Valid() && Get_Range() > 0. )
	{
		double Min = Get_Min(), Range = Get_Range();

		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				Set_Value(x, y, (asDouble(x, y) - Min) / Range);
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CSG_Grid::DeNormalise(double Minimum, double Maximum)
{
	if( is_Valid() && Minimum < Maximum )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				Set_Value(x, y, Minimum + asDouble(x, y) * (Maximum - Minimum));
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CSG_Grid::Standardise(void)
{
	if( is_Valid() && Get_StdDev() > 0. )
	{
		double Mean = Get_Mean(), StdDev = Get_StdDev();

		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				Set_Value(x, y, (asDouble(x, y) - Mean) / StdDev);
			}
		}

		return( true );
	}

	return( false );
}

//---------------------------------------------------------
bool CSG_Grid::DeStandardise(double Mean, double StdDev)
{
	if( is_Valid() && StdDev > 0. )
	{
		#pragma omp parallel for
		for(int y=0; y<Get_NY(); y++) for(int x=0; x<Get_NX(); x++)
		{
			if( !is_NoData(x, y) )
			{
				Set_Value(x, y, Mean + asDouble(x, y) * StdDev);
			}
		}

		return( true );
	}

	return( false );
}


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

//---------------------------------------------------------
/**
  * Returns the direction to which the downward gradient is
  * steepest. This implements the Deterministic 8 (D8) algorithm
  * for identifying a single flow direction based on elevation
  * data. Direction is numbered clock-wise beginning with 0 for
  * the North. If no direction can be identified result will be a -1.
  * If 'bDown' is not true the cell that the direction is pointing
  * to might have a higher value than the center cell. If 'bNoEdges'
  * is true a -1 will be returned for all cells that are at the
  * edge of the data area.
*/
//---------------------------------------------------------
int CSG_Grid::Get_Gradient_NeighborDir(int x, int y, bool bDown, bool bNoEdges)	const
{
	int Direction = -1;

	if( is_InGrid(x, y) )
	{
		double z = asDouble(x, y), dzMax = 0.;

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

			if( is_InGrid(ix, iy) )
			{
				double dz = (z - asDouble(ix, iy)) / m_System.Get_Length(i);

				if( (!bDown || dz > 0.) && (Direction < 0 || dzMax < dz) )
				{
					Direction = i; dzMax = dz;
				}
			}
			else if( bNoEdges )
			{
				return( -1 );
			}
		}
	}

	return( Direction );
}

//---------------------------------------------------------
/**
  * Calculates the gradient of a cell interpreting all grid cell values
  * as elevation surface. Calculation uses the formulas proposed
  * by Zevenbergen & Thorne (1986). Slope and aspect values are calculated
  * in radians. Aspect is zero for the North direction and increases
  * clockwise.
*/
//---------------------------------------------------------
bool CSG_Grid::Get_Gradient(int x, int y, double &Slope, double &Aspect) const
{
	if( is_InGrid(x, y) )
	{
		double z = asDouble(x, y), dz[4];

		for(int i=0, iDir=0, ix, iy; i<4; i++, iDir+=2)
		{
			if( is_InGrid(
				ix = m_System.Get_xTo  (iDir, x),
				iy = m_System.Get_yTo  (iDir, y)) )
			{
				dz[i] = asDouble(ix, iy) - z;
			}
			else if( is_InGrid(
				ix = m_System.Get_xFrom(iDir, x),
				iy = m_System.Get_yFrom(iDir, y)) )
			{
				dz[i] = z - asDouble(ix, iy);
			}
			else
			{
				dz[i] = 0.;
			}
		}

		double G = (dz[0] - dz[2]) / (2. * Get_Cellsize());
        double H = (dz[1] - dz[3]) / (2. * Get_Cellsize());

		Slope  = atan(sqrt(G*G + H*H));
		Aspect = G != 0. ? M_PI_180 + atan2(H, G) : H > 0. ? M_PI_270 : H < 0. ? M_PI_090 : -1.;

		return( true );
	}

	Slope  =  0.;
	Aspect = -1.;

	return( false );
}

//---------------------------------------------------------
/**
  * Calculates the gradient for the given world coordinate.
  * Calculation uses the formulas proposed by Zevenbergen & Thorne (1986).
  * Slope and aspect values are calculated in radians.
  * Aspect is zero for the North direction and increases clockwise.
*/
//---------------------------------------------------------
bool CSG_Grid::Get_Gradient(double x, double y, double &Slope, double &Aspect, CSG_Grid_Resampling Interpolation) const
{
	double z, iz, dz[4];

	if( Get_Value(x, y, z, Interpolation) )
	{
		for(int i=0, iDir=0; i<4; i++, iDir+=2)
		{
			if( Get_Value(
				x + Get_Cellsize() * m_System.Get_xTo  (iDir),
				y + Get_Cellsize() * m_System.Get_yTo  (iDir), iz, Interpolation) )
			{
				dz[i] = iz - z;
			}
			else if( Get_Value(
				x + Get_Cellsize() * m_System.Get_xFrom(iDir),
				y + Get_Cellsize() * m_System.Get_yFrom(iDir), iz, Interpolation) )
			{
				dz[i] = z - iz;
			}
			else
			{
				dz[i] = 0.;
			}
		}

		double G = (dz[0] - dz[2]) / (2. * Get_Cellsize());
        double H = (dz[1] - dz[3]) / (2. * Get_Cellsize());

		Slope  = atan(sqrt(G*G + H*H));
		Aspect = G != 0. ? M_PI_180 + atan2(H, G) : H > 0. ? M_PI_270 : H < 0. ? M_PI_090 : -1.;

		return( true );
	}

	Slope  =  0.;
	Aspect = -1.;

	return( false );
}

//---------------------------------------------------------
/**
  * Calculates the gradient for the given world coordinate.
  * Calculation uses the formulas proposed by Zevenbergen & Thorne (1986).
  * Slope and aspect values are calculated in radians.
  * Aspect is zero for the North direction and increases clockwise.
*/
//---------------------------------------------------------
bool CSG_Grid::Get_Gradient(const TSG_Point &p, double &Incline, double &Azimuth, CSG_Grid_Resampling Interpolation) const
{
	return( Get_Gradient(p.x, p.y, Incline, Azimuth, Interpolation) );
}


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

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