
///////////////////////////////////////////////////////////
//                                                       //
//                         SAGA                          //
//                                                       //
//      System for Automated Geoscientific Analyses      //
//                                                       //
//                     Tool Library                      //
//                      Grid_Shapes                      //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//            Grid_Statistics_AddTo_Polygon.cpp          //
//                                                       //
//                 Copyright (C) 2003 by                 //
//                      Olaf Conrad                      //
//                                                       //
//                 quantile calculation:                 //
//                 Copyright (C) 2007 by                 //
//                   Johan Van de Wauw                   //
//                                                       //
//-------------------------------------------------------//
//                                                       //
// This file is part of 'SAGA - System for Automated     //
// Geoscientific Analyses'. SAGA is free software; you   //
// can redistribute it and/or modify it under the terms  //
// of the GNU General Public License as published by the //
// Free Software Foundation, either version 2 of the     //
// License, or (at your option) any later version.       //
//                                                       //
// SAGA 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 General Public        //
// License for more details.                             //
//                                                       //
// You should have received a copy of the GNU General    //
// Public License along with this program; if not, see   //
// <http://www.gnu.org/licenses/>.                       //
//                                                       //
//-------------------------------------------------------//
//                                                       //
//    e-mail:     oconrad@saga-gis.org                   //
//                                                       //
//    contact:    Olaf Conrad                            //
//                Institute of Geography                 //
//                University of Goettingen               //
//                Goldschmidtstr. 5                      //
//                37077 Goettingen                       //
//                Germany                                //
//                                                       //
///////////////////////////////////////////////////////////

//---------------------------------------------------------
#include "Grid_Statistics_AddTo_Polygon.h"


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

//---------------------------------------------------------
CGrid_Statistics_AddTo_Polygon::CGrid_Statistics_AddTo_Polygon(void)
{
	Set_Name		(_TL("Zonal Grid Statistics for Polygons"));

	Set_Author		("O.Conrad (c) 2003, Quantile Calculation (c) 2007 by Johan Van de Wauw");

	Set_Description	(_TW(
		"Zonal grid statistics. For each polygon statistics based on all covered grid cells will be calculated. "
		"Statistics can also be aggregated by a category provided by an attribute. "
	));

	//-----------------------------------------------------
	Parameters.Add_Grid_List("",
		"GRIDS"       , _TL("Grids"),
		_TL(""),
		PARAMETER_INPUT
	);

	Parameters.Add_Shapes("",
		"POLYGONS"    , _TL("Polygons"),
		_TL(""),
		PARAMETER_INPUT, SHAPE_TYPE_Polygon
	);

	Parameters.Add_Shapes("",
		"RESULT"      , _TL("Polygon Statistics"),
		_TL("If not set polygon statistics will be appended to the input polygons layer."),
		PARAMETER_OUTPUT_OPTIONAL, SHAPE_TYPE_Polygon
	);

	Parameters.Add_Bool("RESULT",
		"ATTRIBUTES"  , _TL("Copy Attributes"),
		_TL("Copy attributes from input polygons layer."),
		false
	);

	Parameters.Add_Table_Field("POLYGONS",
		"CATEGORY"    , _TL("Category"),
		_TL(""),
		false
	);

	Parameters.Add_Table("",
		"CATEGORIES"  , _TL("Category Statistics"),
		_TL(""),
		PARAMETER_OUTPUT
	);

	Parameters.Add_Choice("",
		"OUTPUT"      , _TL("Output"),
		_TL(""),
		CSG_String::Format("%s|%s|%s",
			_TL("polygon statistics"),
			_TL("polygon and category statistics"),
			_TL("category statistics")
		), 0
	);

	Parameters.Add_Choice("",
		"NAMING"      , _TL("Field Naming"),
		_TL(""),
		CSG_String::Format("%s|%s",
			_TL("grid number"),
			_TL("grid name")
		), 1
	);

	Parameters.Add_Choice("",
		"METHOD"      , _TL("Method"),
		_TL(""),
		CSG_String::Format("%s|%s|%s|%s",
			_TL("simple and fast"),
			_TL("polygon wise (cell centers)"),
			_TL("polygon wise (cell area)"),
			_TL("polygon wise (cell area weighted)")
		), 0
	);

	Parameters.Add_Bool("METHOD",
		"PARALLELIZED", _TL("Use Multiple Cores"),
		_TL(""),
		true
	);

	//-----------------------------------------------------
	Parameters.Add_Node("",
		"STATISTICS"  , _TL("Statistics"),
		_TL("")
	);

	Parameters.Add_Bool  ("STATISTICS", "COUNT"    , _TL("Number of Cells"   ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "MIN"      , _TL("Minimum"           ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "MAX"      , _TL("Maximum"           ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "RANGE"    , _TL("Range"             ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "SUM"      , _TL("Sum"               ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "MEAN"     , _TL("Mean"              ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "VAR"      , _TL("Variance"          ), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "STDDEV"   , _TL("Standard Deviation"), _TL(""),  true);
	Parameters.Add_Bool  ("STATISTICS", "GINI"     , _TL("Gini"              ), _TL(""), false);
	Parameters.Add_String("STATISTICS", "QUANTILES", _TL("Percentiles"       ),
		_TL("Separate the desired percentiles by semicolon, e.g. \"5; 25; 50; 75; 95\""), ""
	);
}


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

//---------------------------------------------------------
int CGrid_Statistics_AddTo_Polygon::On_Parameters_Enable(CSG_Parameters *pParameters, CSG_Parameter *pParameter)
{
	if( pParameter->Cmp_Identifier("RESULT") )
	{
		pParameters->Set_Enabled("ATTRIBUTES"  , pParameter->asPointer());
	}

	if( pParameter->Cmp_Identifier("OUTPUT") )
	{
		pParameters->Set_Enabled("RESULT"      , pParameter->asInt() <= 1);
		pParameters->Set_Enabled("CATEGORY"    , pParameter->asInt() >= 1);
		pParameters->Set_Enabled("CATEGORIES"  , pParameter->asInt() >= 1);
	}

	if( pParameter->Cmp_Identifier("METHOD") )
	{
		pParameters->Set_Enabled("PARALLELIZED", pParameter->asInt() != 0 && SG_OMP_Get_Max_Num_Threads() > 1);
	}

	return( CSG_Tool_Grid::On_Parameters_Enable(pParameters, pParameter) );
}


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

//---------------------------------------------------------
bool CGrid_Statistics_AddTo_Polygon::On_Execute(void)
{
	CSG_Shapes *pPolygons = Parameters["POLYGONS"].asShapes();

	if( pPolygons->Get_Count() < 1 )
	{
		Error_Set(_TL("no polygons in input layer"));

		return( false );
	}

	//-----------------------------------------------------
	CSG_Parameter_Grid_List *pGrids = Parameters["GRIDS"].asGridList();

	if( pGrids->Get_Grid_Count() < 1 )
	{
		Error_Set(_TL("no grids in selection"));

		return( false );
	}

	if( !Get_System().Get_Extent().Intersects(pPolygons->Get_Extent()) )
	{
		Error_Set(_TL("no spatial intersection between grid(s) and polygon layer"));

		return( false );
	}

	//-----------------------------------------------------
	int	nVars   = 0;

	int	fCOUNT  = Parameters["COUNT" ].asBool() ? nVars++ : -1;
	int	fMIN    = Parameters["MIN"   ].asBool() ? nVars++ : -1;
	int	fMAX    = Parameters["MAX"   ].asBool() ? nVars++ : -1;
	int	fRANGE  = Parameters["RANGE" ].asBool() ? nVars++ : -1;
	int	fSUM    = Parameters["SUM"   ].asBool() ? nVars++ : -1;
	int	fMEAN   = Parameters["MEAN"  ].asBool() ? nVars++ : -1;
	int	fVAR    = Parameters["VAR"   ].asBool() ? nVars++ : -1;
	int	fSTDDEV = Parameters["STDDEV"].asBool() ? nVars++ : -1;
	int	fGINI   = Parameters["GINI"  ].asBool() ? nVars++ : -1;

	CSG_Vector Percentiles;

	{
		CSG_Strings Values = SG_String_Tokenize(Parameters["QUANTILES"].asString(), ";");

		for(int i=0; i<Values.Get_Count(); i++)
		{
			double Value;

			if( Values[i].asDouble(Value) && Value >= 0. && Value <= 100. )
			{
				Percentiles.Add_Row(Value);
			}
		}
	}

	int fQUANTILE = Percentiles.Get_N() > 0 ? nVars : -1;

	nVars += Percentiles.Get_N();

	if( nVars == 0 )
	{
		Error_Set(_TL("no statistical variable in selection"));

		return( false );
	}

	//-----------------------------------------------------
	bool bParallelized = Parameters["PARALLELIZED"].is_Enabled() && Parameters["PARALLELIZED"].asBool();

	int Naming = Parameters["NAMING"].asInt();
	int Method = Parameters["METHOD"].asInt();

	CSG_Grid Index;

	if( Method == 0 && !Get_Simple_Index(pPolygons, Index) )
	{
		Error_Set(_TL("no grids in selection"));

		return( false );
	}

	//-----------------------------------------------------
	CSG_Shapes *pOutput = Parameters["OUTPUT"].asInt() <= 1 ? Parameters["RESULT"].asShapes() : NULL;

	if( pOutput && pOutput != pPolygons )
	{
		if( Parameters["ATTRIBUTES"].asBool() ) // copy construct
		{
			pOutput->Create(*pPolygons);
		}
		else // copy polygons
		{
			pOutput->Create(SHAPE_TYPE_Polygon, NULL, NULL, pPolygons->Get_Vertex_Type());

			for(sLong i=0; i<pPolygons->Get_Count(); i++)
			{
				pOutput->Add_Shape(pPolygons->Get_Shape(i));
			}
		}

		pOutput->Fmt_Name("%s [%s]", Parameters["POLYGONS"].asShapes()->Get_Name(), _TL("Grid Statistics"));
	}
	else if( Parameters["OUTPUT"].asInt() <= 1 ) // append statistics to input polygons
	{
		pOutput = pPolygons;
	}

	//-----------------------------------------------------
	CSG_Table *pCategories = NULL; int Category = Parameters["CATEGORY"].asInt();

	if( Parameters["OUTPUT"].asInt() >= 1 && Category >= 0 && Category < pPolygons->Get_Field_Count() )
	{
		CSG_Unique_String_Statistics Uniques;

		for(sLong i=0; i<pPolygons->Get_Count(); i++)
		{
			Uniques += pPolygons->Get_Shape(i)->asString(Category);
		}

		if( Uniques.Get_Count() > 0 )
		{
			pCategories = Parameters["CATEGORIES"].asTable();
			pCategories->Create();
			pCategories->Fmt_Name("%s.%s [%s]", pPolygons->Get_Name(), pPolygons->Get_Field_Name(Category), _TL("Grid Statistics"));
			pCategories->Add_Field(pPolygons->Get_Field_Name(Category), SG_DATATYPE_String);

			for(int i=0; i<Uniques.Get_Count(); i++)
			{
				pCategories->Add_Record()->Set_Value(0, Uniques.Get_Value(i));
			}
		}
	}

	//-----------------------------------------------------
	for(int iGrid=0; iGrid<pGrids->Get_Grid_Count() && Process_Get_Okay(); iGrid++)
	{
		Process_Set_Text("[%d/%d] %s", 1 + iGrid, pGrids->Get_Grid_Count(), pGrids->Get_Grid(iGrid)->Get_Name());

		CSG_Simple_Statistics *Statistics = new CSG_Simple_Statistics[pPolygons->Get_Count()], *Categories = pCategories ? new CSG_Simple_Statistics[pCategories->Get_Count()] : NULL;

		if( (Method == 0 && Get_Simple (pGrids->Get_Grid(iGrid), pPolygons, Statistics, Percentiles.Get_N() > 0 || fGINI > 0, Index        ))
		||  (Method != 0 && Get_Precise(pGrids->Get_Grid(iGrid), pPolygons, Statistics, Percentiles.Get_N() > 0 || fGINI > 0, bParallelized)) )
		{
			int nOutFields = pOutput ? pOutput->Get_Field_Count() : 0, nCatFields = pCategories ? pCategories->Get_Field_Count() : 0;

			#define ADD_FIELD(var, type) { CSG_String name(var);\
				if( !Naming ) { name.Printf("G%02d_%s", iGrid + 1, name.c_str()); }\
				else          { name.Printf("%s (%s)", pGrids->Get_Grid(iGrid)->Get_Name(), name.c_str()); }\
				if( pOutput ) { pOutput->Add_Field(name, type); } if( pCategories ) { pCategories->Add_Field(name, type); }\
			}

			if( fCOUNT    >= 0 ) { ADD_FIELD("CELLS"   , SG_DATATYPE_Int   ); }
			if( fMIN      >= 0 ) { ADD_FIELD("MIN"     , SG_DATATYPE_Double); }
			if( fMAX      >= 0 ) { ADD_FIELD("MAX"     , SG_DATATYPE_Double); }
			if( fRANGE    >= 0 ) { ADD_FIELD("RANGE"   , SG_DATATYPE_Double); }
			if( fSUM      >= 0 ) { ADD_FIELD("SUM"     , SG_DATATYPE_Double); }
			if( fMEAN     >= 0 ) { ADD_FIELD("MEAN"    , SG_DATATYPE_Double); }
			if( fVAR      >= 0 ) { ADD_FIELD("VARIANCE", SG_DATATYPE_Double); }
			if( fSTDDEV   >= 0 ) { ADD_FIELD("STDDEV"  , SG_DATATYPE_Double); }
			if (fGINI     >= 0 ) { ADD_FIELD("GINI"    , SG_DATATYPE_Double); }
			if( fQUANTILE >= 0 )
			{
				for(int iPercentile=0; iPercentile<Percentiles.Get_N(); iPercentile++)
				{
					ADD_FIELD(CSG_String("PCTL" + SG_Get_String(Percentiles[iPercentile], -2)).c_str(), SG_DATATYPE_Double);
				}
			}

			//---------------------------------------------
			for(sLong i=0; i<pPolygons->Get_Count() && Set_Progress(i, pPolygons->Get_Count()); i++)
			{
				if( pOutput )
				{
					CSG_Shape *pPolygon = pOutput->Get_Shape(i);

					if( Statistics[i].Get_Count() == 0 )
					{
						if( fCOUNT    >= 0 ) { pPolygon->Set_NoData(nOutFields + fCOUNT ); }
						if( fMIN      >= 0 ) { pPolygon->Set_NoData(nOutFields + fMIN   ); }
						if( fMAX      >= 0 ) { pPolygon->Set_NoData(nOutFields + fMAX   ); }
						if( fRANGE    >= 0 ) { pPolygon->Set_NoData(nOutFields + fRANGE ); }
						if( fSUM      >= 0 ) { pPolygon->Set_NoData(nOutFields + fSUM   ); }
						if( fMEAN     >= 0 ) { pPolygon->Set_NoData(nOutFields + fMEAN  ); }
						if( fVAR      >= 0 ) { pPolygon->Set_NoData(nOutFields + fVAR   ); }
						if( fSTDDEV   >= 0 ) { pPolygon->Set_NoData(nOutFields + fSTDDEV); }
						if( fGINI     >= 0 ) { pPolygon->Set_NoData(nOutFields + fGINI  ); }
						if( fQUANTILE >= 0 )
						{
							for(int iPercentile=0, iField=nOutFields + fQUANTILE; iPercentile<Percentiles.Get_N(); iPercentile++, iField++)
							{
								pPolygon->Set_NoData(iField);
							}
						}
					}
					else
					{
						if( fCOUNT    >= 0 ) { pPolygon->Set_Value(nOutFields + fCOUNT , Statistics[i].Get_Count   ()); }
						if( fMIN      >= 0 ) { pPolygon->Set_Value(nOutFields + fMIN   , Statistics[i].Get_Minimum ()); }
						if( fMAX      >= 0 ) { pPolygon->Set_Value(nOutFields + fMAX   , Statistics[i].Get_Maximum ()); }
						if( fRANGE    >= 0 ) { pPolygon->Set_Value(nOutFields + fRANGE , Statistics[i].Get_Range   ()); }
						if( fSUM      >= 0 ) { pPolygon->Set_Value(nOutFields + fSUM   , Statistics[i].Get_Sum     ()); }
						if( fMEAN     >= 0 ) { pPolygon->Set_Value(nOutFields + fMEAN  , Statistics[i].Get_Mean    ()); }
						if( fVAR      >= 0 ) { pPolygon->Set_Value(nOutFields + fVAR   , Statistics[i].Get_Variance()); }
						if( fSTDDEV   >= 0 ) { pPolygon->Set_Value(nOutFields + fSTDDEV, Statistics[i].Get_StdDev  ()); }
						if( fGINI     >= 0 ) { pPolygon->Set_Value(nOutFields + fGINI  , Statistics[i].Get_Gini    ()); }
						if( fQUANTILE >= 0 )
						{
							for(int iPercentile=0, iField=nOutFields + fQUANTILE; iPercentile<Percentiles.Get_N(); iPercentile++, iField++)
							{
								pPolygon->Set_Value(iField, Statistics[i].Get_Percentile(Percentiles[iPercentile]));
							}
						}
					}
				}

				if( pCategories && Statistics[i].Get_Count() > 0 )
				{
					sLong iCategory = -1; CSG_String sCategory(pPolygons->Get_Shape(i)->asString(Category));

					for(sLong j=0; iCategory<0 && j<pCategories->Get_Count(); j++)
					{
						if( sCategory.Cmp(pCategories->Get_Record(j)->asString(0)) == 0 )
						{
							iCategory = j;
						}
					}

					if( iCategory >= 0 )
					{
						Categories[iCategory] += Statistics[i];
					}
				}
			}

			//---------------------------------------------
			if( pCategories )
			{
				for(sLong i=0; i<pCategories->Get_Count(); i++)
				{
					CSG_Table_Record *pCategory = pCategories->Get_Record(i);

					if( fCOUNT    >= 0 ) { pCategory->Set_Value(nCatFields + fCOUNT , Categories[i].Get_Count   ()); }
					if( fMIN      >= 0 ) { pCategory->Set_Value(nCatFields + fMIN   , Categories[i].Get_Minimum ()); }
					if( fMAX      >= 0 ) { pCategory->Set_Value(nCatFields + fMAX   , Categories[i].Get_Maximum ()); }
					if( fRANGE    >= 0 ) { pCategory->Set_Value(nCatFields + fRANGE , Categories[i].Get_Range   ()); }
					if( fSUM      >= 0 ) { pCategory->Set_Value(nCatFields + fSUM   , Categories[i].Get_Sum     ()); }
					if( fMEAN     >= 0 ) { pCategory->Set_Value(nCatFields + fMEAN  , Categories[i].Get_Mean    ()); }
					if( fVAR      >= 0 ) { pCategory->Set_Value(nCatFields + fVAR   , Categories[i].Get_Variance()); }
					if( fSTDDEV   >= 0 ) { pCategory->Set_Value(nCatFields + fSTDDEV, Categories[i].Get_StdDev  ()); }
					if( fGINI     >= 0 ) { pCategory->Set_Value(nCatFields + fGINI  , Categories[i].Get_Gini    ()); }
					if( fQUANTILE >= 0 )
					{
						for(int iPercentile=0, iField=nCatFields + fQUANTILE; iPercentile<Percentiles.Get_N(); iPercentile++, iField++)
						{
							pCategory->Set_Value(iField, Categories[i].Get_Percentile(Percentiles[iPercentile]));
						}
					}
				}
			}
		}

		//-------------------------------------------------
		delete[](Statistics);

		if( Categories )
		{
			delete[](Categories);
		}
	}

	//-----------------------------------------------------
	DataObject_Update(pPolygons);

	return( true );
}


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

//---------------------------------------------------------
bool CGrid_Statistics_AddTo_Polygon::Get_Precise(CSG_Grid *pGrid, CSG_Shapes *pPolygons, CSG_Simple_Statistics *Statistics, bool bHoldValues, bool bParallelized)
{
	int Method = Parameters["METHOD"].asInt();

	if( bParallelized )
	{
		#pragma omp parallel for
		for(sLong i=0; i<pPolygons->Get_Count(); i++)
		{
			Get_Precise(pGrid, (CSG_Shape_Polygon *)pPolygons->Get_Shape(i), Statistics[i], bHoldValues, Method);
		}
	}
	else
	{
		for(sLong i=0; i<pPolygons->Get_Count() && Set_Progress(i, pPolygons->Get_Count()); i++)
		{
			Get_Precise(pGrid, (CSG_Shape_Polygon *)pPolygons->Get_Shape(i), Statistics[i], bHoldValues, Method);
		}
	}

	return( true );
}

//---------------------------------------------------------
bool CGrid_Statistics_AddTo_Polygon::Get_Precise(CSG_Grid *pGrid, CSG_Shape_Polygon *pPolygon, CSG_Simple_Statistics &Statistics, bool bHoldValues, int Method)
{
	CSG_Shapes Intersect(SHAPE_TYPE_Polygon);

	CSG_Shape_Polygon *pCell = Method != 3 ? NULL : Intersect.Add_Shape()->asPolygon();
	CSG_Shape_Polygon *pArea = Method != 3 ? NULL : Intersect.Add_Shape()->asPolygon();

	//-----------------------------------------------------
	Statistics.Create(bHoldValues);

	int ax = Get_System().Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMin()) - 1; if( ax < 0         ) { ax = 0           ; }
	int bx = Get_System().Get_xWorld_to_Grid(pPolygon->Get_Extent().Get_XMax()) + 1; if( bx >= Get_NX() ) { bx = Get_NX() - 1; }
	int ay = Get_System().Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMin()) - 1; if( ay < 0         ) { ay = 0           ; }
	int by = Get_System().Get_yWorld_to_Grid(pPolygon->Get_Extent().Get_YMax()) + 1; if( by >= Get_NY() ) { by = Get_NY() - 1; }

	TSG_Point Center; TSG_Rect Cell;

	//-----------------------------------------------------
	Center.y  = Get_System().Get_yGrid_to_World(ay);
	Cell.yMin = Center.y - 0.5 * Get_Cellsize();
	Cell.yMax = Cell.yMin + Get_Cellsize();
		
	for(int y=ay; y<=by && Process_Get_Okay(); y++, Center.y+=Get_Cellsize(), Cell.yMin+=Get_Cellsize(), Cell.yMax+=Get_Cellsize())
	{
		Center.x  = Get_System().Get_xGrid_to_World(ax);
		Cell.xMin = Center.x - 0.5 * Get_Cellsize();
		Cell.xMax = Cell.xMin + Get_Cellsize();

		for(int x=ax; x<=bx; x++, Center.x+=Get_Cellsize(), Cell.xMin+=Get_Cellsize(), Cell.xMax+=Get_Cellsize())
		{
			if( !pGrid->is_NoData(x, y) )
			{
				switch( Method )
				{
				//-----------------------------------------
				default:	// polygon wise (cell centers)
					if( pPolygon->Contains(Center) )
					{
						Statistics += pGrid->asDouble(x, y);
					}
					break;

				//-----------------------------------------
				case  2:	// polygon wise (cell area)
					if( pPolygon->Intersects(Cell) )
					{
						Statistics += pGrid->asDouble(x, y);
					}
					break;

				//-----------------------------------------
				case  3:	// polygon wise (cell area weighted)
					switch( pPolygon->Intersects(Cell) )
					{
					case INTERSECTION_None     : break;
					case INTERSECTION_Identical:
					case INTERSECTION_Contains : Statistics.Add_Value(pGrid->asDouble(x, y), 1.);	break;
					case INTERSECTION_Contained: Statistics.Add_Value(pGrid->asDouble(x, y), pPolygon->Get_Area() / Get_Cellarea());	break;
					case INTERSECTION_Overlaps :
						pCell->Del_Parts();

						pCell->Add_Point(Cell.xMin, Cell.yMin);	pCell->Add_Point(Cell.xMin, Cell.yMax);
						pCell->Add_Point(Cell.xMax, Cell.yMax);	pCell->Add_Point(Cell.xMax, Cell.yMin);

						if( SG_Shape_Get_Intersection(pPolygon, pCell, pArea) )
						{
							Statistics.Add_Value(pGrid->asDouble(x, y), pArea->Get_Area() / Get_Cellarea());
						}
						break;
					}
					break;
				}
			}
		}
	}

	return( true );
}


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

//---------------------------------------------------------
bool CGrid_Statistics_AddTo_Polygon::Get_Simple(CSG_Grid *pGrid, CSG_Shapes *pPolygons, CSG_Simple_Statistics *Statistics, bool bHoldValues, CSG_Grid &Index)
{
	for(sLong i=0; i<pPolygons->Get_Count(); i++)
	{
		Statistics[i].Create(bHoldValues);
	}

	for(int y=0; y<Get_NY() && Set_Progress_Rows(y); y++)
	{
		for(int x=0, i; x<Get_NX(); x++)
		{
			if( !pGrid->is_NoData(x, y) && (i = Index.asInt(x, y)) >= 0 && i < pPolygons->Get_Count() )
			{
				Statistics[i] += pGrid->asDouble(x, y);
			}
		}
	}

	return( true );
}

//---------------------------------------------------------
bool CGrid_Statistics_AddTo_Polygon::Get_Simple_Index(CSG_Shapes *pPolygons, CSG_Grid &Index)
{
	Index.Create(Get_System(), pPolygons->Get_Count() < 32767 ? SG_DATATYPE_Short : SG_DATATYPE_Int);
	Index.Assign(-1.);

	bool *bCrossing = (bool *)SG_Malloc(Get_NX() * sizeof(bool));

	//-----------------------------------------------------
	for(sLong iShape=0; iShape<pPolygons->Get_Count() && Set_Progress(iShape, pPolygons->Get_Count()); iShape++)
	{
		CSG_Shape *pPolygon = pPolygons->Get_Shape(iShape); TSG_Rect Extent = pPolygon->Get_Extent();

		int xStart = Get_System().Get_xWorld_to_Grid(Extent.xMin) - 1;	if( xStart < 0 )		xStart	= 0;
		int xStop  = Get_System().Get_xWorld_to_Grid(Extent.xMax) + 1;	if( xStop >= Get_NX() )	xStop	= Get_NX() - 1;

		CSG_Point Left (Get_XMin() - 1., 0.);
		CSG_Point Right(Get_XMax() + 1., 0.);

		//-------------------------------------------------
		for(int y=0; y<Get_NY(); y++)
		{
			double yPos = Get_YMin() + y * Get_Cellsize();

			if( yPos >= Extent.yMin && yPos <= Extent.yMax )
			{
				memset(bCrossing, 0, Get_NX() * sizeof(bool));

				Left.y = Right.y = yPos;

				//-----------------------------------------
				for(int iPart=0; iPart<pPolygon->Get_Part_Count(); iPart++)
				{
					CSG_Point B = pPolygon->Get_Point(pPolygon->Get_Point_Count(iPart) - 1, iPart);

					for(int iPoint=0; iPoint<pPolygon->Get_Point_Count(iPart); iPoint++)
					{
						CSG_Point A = B; B = pPolygon->Get_Point(iPoint, iPart);

						if( ((A.y <= yPos && yPos <  B.y)
						||   (A.y >  yPos && yPos >= B.y)) )
						{
							TSG_Point C; SG_Get_Crossing(C, A, B, Left, Right, false);

							int ix = (int)((C.x - Get_XMin()) / Get_Cellsize() + 1.);

							if( ix < 0 )
							{
								ix = 0;
							}
							else if( ix >= Get_NX() )
							{
								continue;
							}

							bCrossing[ix] = !bCrossing[ix];
						}
					}
				}

				//-----------------------------------------
				for(int x=xStart, bFill=0; x<=xStop; x++)
				{
					if( bCrossing[x] )
					{
						bFill = bFill ? 0 : 1;
					}

					if( bFill )
					{
						Index.Set_Value(x, y, (double)iShape);
					}
				}
			}
		}
	}

	//-----------------------------------------------------
	SG_Free(bCrossing);

	return( true );
}


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

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