/* 
* 
* Copyright (C) 2013 UMPC
* 
* This file is part of CLaRiNet.
*  
* CLaRiNet 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 3 of the License, or
* (at your option) any later version.
* 
* CLaRiNet 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 CLaRiNet.  If not, see <http://www.gnu.org/licenses/>.
*/
#include <stdlib.h>
#include <stddef.h>
//#include <sys/resource.h>
#include "GeoMorph.h"
#include "RasterCalc.h"
#include "Erreur.h"
#include "esriheader.h"
#include "FlowDirGrid.h"
#include "FloatRaster.h"
#include "Topology.h"
#include "CellTree.h"

#include<map>
#include<set>

using namespace std;

int GetRasterValue ( const int HandleRaster , int * CELLS , const int ncells , float ** _VALUES )
{
    int i,nr,nc;
    FloatRaster * R;

    if ( !(R = GetFloatRaster(HandleRaster)))
    {
        SetErreur("GetRasterValue : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    for (i=0;i<ncells;i++)
    {
        (void) rec2RC ( &(R->Header) , CELLS + i , & nr, & nc );
        (*_VALUES)[i] = R->Data[nr][nc];
    }
    
    return(0);
}

int RasterFind ( const int HandleRaster , const float * _a, const float * _b, int * _nfound , int ** _CELLS )
{
     /* Trouve les cellules irec pour lesquelles la condition
                   a*Raster[irec]+b >0
     est vÃƒÂ©rifiÃƒÂ©e
     */
    
    int nr,nc,nrows,ncols,nfound;
    FloatRaster * R;
    float f,a,b;

    if ( !(R = GetFloatRaster(HandleRaster)))
    {
        SetErreur("RasterFind : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    ncols = R->Header.ncols;
    nrows = R->Header.nrows;

    /* Initialisation des valeurs de retour */

    nfound = 0;
    *_CELLS = NULL;
    a = *_a;
    b = *_b;

    for (nr=0;nr<nrows;nr++)
    {
        for (nc=0;nc<ncols;nc++)
        {
            f = R->Data[nr][nc];

            if ( (a*f+b)>0 )
            {
                if ( !(*_CELLS = (int *) realloc(*_CELLS , (nfound+1)*sizeof(int)) ) )
                {
                    SetErreur("RasterFind : could not reallocate integer pointer for output\n");
                    return(-1);
                }

                (*_CELLS)[nfound] = nr*ncols+nc;
                nfound++;
            }
        }
    }

    *_nfound = nfound;

    return(0);
}


int RasterFindTransition_ ( FloatRaster * R , FlowDirGrid * F , const float * _a , const float * _b , int * _root , float * _evalexpr , int * _nit , set<int> & Cells )
{
    int i,nup1 , nit_i ;
    int nr,nc ;
    int * UPSTREAM ;
    float evalexpr_i;
        
    if ( (*_nit)<=0 )
        return(0);


    /* RÃƒÂ©cupÃƒÂ©ration des indices des cellules amont */

    if ( GetUpstream ( F , _root , & nup1 , & UPSTREAM ) )
    {
            SetErreur("RasterFindTransition_ : could not get upstream cells\n");
            return(-1);
    }

    if(!nup1)
        return(0);

    /* Parcours des cellules amont */

    for(i=0;i<nup1;i++)
    {
        /* Calcul de la valeur de a*Raster[UPSTREAM[i]]+b */

        (void) rec2RC ( &(R->Header) , UPSTREAM + i , & nr, & nc );
        
	evalexpr_i = (*_a) * (R->Data[nr][nc]) + (*_b);	


	if ( (! ((*_evalexpr)>0) ) & ( evalexpr_i>0 ) )
        {
            nit_i=(*_nit)-1;
            Cells.insert(UPSTREAM[i]);
        }
        else
        {
            nit_i=(*_nit);
        }

        /* Appel récursif */

        if ( RasterFindTransition_ ( R , F , _a , _b , UPSTREAM + i , & evalexpr_i , & nit_i , Cells ) )
        {
            SetErreur("RasterFindTransition_ : recursive call failed\n");
            return(-1);
        }

     }

    free(UPSTREAM);

    return (0);
}


int RasterFindTransition (  const int HandleRaster ,
                            const int HandleFlowDir ,
                            const float * _a,
                            const float * _b,
                            int * _root,
                            int * _nitmax,
                            int * _nfound,
                            int ** _CELLS )
{
    /* Trouve les cellules irec pour lesquelles la condition
                   a*Raster[irec]+b >0
     est verifiee alors qu'elle ne l'est par pour la cellule aval.
     L'analyse est effectuée de façon récursive à partir de la cellule "root"
     */

    FloatRaster * R;
    FlowDirGrid * F;

    float f;
    int nr,nc,nit,nfound,i;
    float evalexpr;

    set<int>Cells ;
    set<int>::iterator isi;

//// Trying to fix recursion depth issue ... /////////////////////

/*    const rlim_t kStackSize = 512L * 1024L * 1024L;   // min stack size = 512 Mb
    struct rlimit rl;
    int result;

    result = getrlimit(RLIMIT_STACK, &rl);
    if (result == 0)
    {
        if (rl.rlim_cur < kStackSize)
        {
            rl.rlim_cur = kStackSize;
            result = setrlimit(RLIMIT_STACK, &rl);
            if (result != 0)
            {
                SetErreur("setrlimit returned result = %d\n", result);
		return(-1);
            }
        }
    } 
    else
    {
        SetErreur("getrlimit returned result = %d\n", result);
	return(-1);
    }
*/
///////////////////////////////////////////////////////////////////////////////////


    if ( !(R = GetFloatRaster(HandleRaster)) )
    {
        SetErreur("RasterFindTransition : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    if ( !(F = GetFlowDirGrid(HandleFlowDir)) )
    {
        SetErreur("RasterFindTransition : flow direction grid # %d does not exist\n",HandleFlowDir);
        return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("RasterFindTransition : raster # %d and flow direction grid # %d do not have identical referencing headers\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    /* Initialisation des valeurs de retour */
    
    (void) rec2RC ( &(R->Header) , _root , & nr, & nc );
    f = R->Data[nr][nc];
    evalexpr = (*_a)*f+(*_b);
    Cells.clear();

    if ( RasterFindTransition_ ( R , F , _a , _b , _root , & evalexpr , _nitmax , Cells ) )
    {
        SetErreur("RasterFindTransition : recursive call failed\n");
	return(-1);
    }

    (*_nfound) = Cells.size();
    if ( !((*_CELLS) = (int *) malloc ( (*_nfound) * sizeof(int) ) ) )
    {
     	SetErreur("RasterFindTransition : could not allocate integer pointer for output\n");
	return(-1);
    }
    
    i=0;
    for (isi=Cells.begin();isi!=Cells.end();isi++)
    {
	(*_CELLS)[i] = (*isi);
	i++;
    }

    Cells.clear();
    return(0);
}

typedef struct DistanceMapToken {

	int id ;
	int parentoffset ;
	int nr ;
	int nc ;
	int n1 ;
	int nsqrt2 ;

} DistanceMapToken ;


int MakeDistanceMap_iter ( int HandleFlowDir, int HandleRaster, int outlet )
{
	FlowDirGrid * F;
	FloatRaster * D;
	int i,nup,stacksize,parentoffset,currentoffset,prevoffset,nchildren ;
	int push ;
	int nr,nc,nr_i,nc_i;
	int * UPSTREAM ;
	DistanceMapToken * Stack , * Stack_tmp ;

	/* Initialisation : erreur et valeurs de retour en cas d'erreur */

	ClearErreur();

	if ( !(F = GetFlowDirGrid(HandleFlowDir)) )
	{
		SetErreur("MakeDistanceMap_iter : flow direction grid %d does not exist\n",HandleFlowDir);
		return(-1);
	}

	if ( !(D = GetFloatRaster(HandleRaster)) )
	{
		SetErreur("MakeDistanceMap_iter : float raster %d does not exist\n",HandleRaster);
		return(-1);
	}

	// Initialize stack

	if ( !(Stack = (DistanceMapToken *) malloc ( sizeof(DistanceMapToken) )) )
	{
		SetErreur("MakeDistanceMap_iter - could not allocate stack\n");
		return(-1);
	}
	
	stacksize = 1;
	Stack[0].id = outlet ;
	Stack[0].parentoffset = 0 ;
	(void) rec2RC ( &(F->Header) , & outlet , & nr , & nc );
	Stack[0].nr = nr ;
	Stack[0].nc = nc ;
	Stack[0].n1 = 0 ;
	Stack[0].nsqrt2 = 0 ;

	currentoffset = 0 ;
	push = 1 ;
	UPSTREAM = NULL ;

	while ( stacksize>0 )
	{
		while ( push ) 
		{		
			GetUpstream ( F , & (Stack[currentoffset].id) , & nup , & UPSTREAM ) ;
			nr = Stack[currentoffset].nr ;
			nc = Stack[currentoffset].nc ;

			if ( nup>0) 
			{
				// Push children on stack

				Stack_tmp = (DistanceMapToken *) realloc ( Stack , (stacksize+nup)*sizeof(DistanceMapToken) );
				if (Stack_tmp==NULL)
				{
					SetErreur("MakeDistanceMap_iter - realloc failed, could not push elements on stack\n");
					return(-1);
				}
				else Stack = Stack_tmp ;

				for(i=0;i<nup;i++)
				{
					Stack[stacksize+i].id = UPSTREAM[i];
					(void) rec2RC ( &(F->Header) , UPSTREAM+i , & nr_i , & nc_i );
					Stack[stacksize+i].nr = nr_i ;
					Stack[stacksize+i].nc = nc_i ;
					if ( nr==nr_i | nc==nc_i)
					{
					    Stack[stacksize+i].n1 = Stack[currentoffset].n1 + 1;
					    Stack[stacksize+i].nsqrt2 = Stack[currentoffset].nsqrt2 ;
					}   
					else
					{
					    Stack[stacksize+i].n1 = Stack[currentoffset].n1;
					    Stack[stacksize+i].nsqrt2 = Stack[currentoffset].nsqrt2 +1;
					}

					if (i>0)
						Stack[stacksize+i].parentoffset = -1 ;
					else
						Stack[stacksize+i].parentoffset = currentoffset ;				
				}	
				stacksize+=nup;
				currentoffset = stacksize-1;
				push = 1 ;
				free(UPSTREAM);
			}
			else
			{
				push = 0 ;
			}
		} // end of while ( push ) sequence


		// ******  Compute result here (token is leaf) ***********
		// Result = write distance value to raster

		D->Data[nr][nc] = ((float)(Stack[currentoffset].n1)) + sqrt(2.0)*((float)(Stack[currentoffset].nsqrt2));		

		if ( Stack[currentoffset].parentoffset < 0 )
		{
			currentoffset-- ;
			push = 1 ;
		}

		// ******  Begin pop sequence **********

		while ( (!push) & (stacksize>1) )
		{
			parentoffset = Stack[currentoffset].parentoffset ;
			nchildren = stacksize-currentoffset ;
				
			// Compute the result for the parent, children are at offsets currentoffset:(currentoffset+nchildren-1)

			nr = Stack[parentoffset].nr ;
			nc = Stack[parentoffset].nc ;
			D->Data[nr][nc] = ((float)(Stack[parentoffset].n1)) + sqrt(2.0)*((float)(Stack[parentoffset].nsqrt2));

			// Remove all tokens with offset >= currentoffset from the stack ;
					
			Stack_tmp = (DistanceMapToken *) realloc ( Stack , currentoffset * sizeof(DistanceMapToken) );
			if (Stack_tmp==NULL)
			{
				SetErreur("MakeDistanceMap_iter - realloc failed, could not pop elements from stack\n");
				return(-1);
			}
			else Stack = Stack_tmp ;

			stacksize = currentoffset ;
			currentoffset = parentoffset ;			

			if ( Stack[currentoffset].parentoffset < 0 )
			{
				currentoffset-- ;
				push = 1 ;
			}
		}

		if(stacksize==1)
		{
			stacksize = 0;
		}
	
	}
	
	free(Stack);
	free(UPSTREAM);
	return(0);
}


int MakeDistanceMap_ ( FlowDirGrid * F , FloatRaster * Dst , int * _outlet , int * _n1 , int * _nsqrt2 )
{
    int outlet_i , i , nup1 , nr , nc , nr_i , nc_i ;
    int n1_i, nsqrt2_i;
    float D2;
    int * UPSTREAM;

    (void) rec2RC ( &(F->Header) , _outlet , & nr , & nc );
    Dst->Data[nr][nc] = ((float)(*_n1)) + sqrt(2.0)*((float)(*_nsqrt2));

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("MakeDistanceMap_ : could not get upstream cells\n");
            return(-1);
    }

    for(i=0;i<nup1;i++)
    {
        outlet_i = UPSTREAM[i];
        (void) rec2RC ( &(F->Header) , & outlet_i , & nr_i , & nc_i );

        if ( nr==nr_i | nc==nc_i)
        {
            n1_i=(*_n1)+1;
            nsqrt2_i = (*_nsqrt2);
        }   
        else
        {
            n1_i=(*_n1);
            nsqrt2_i = (*_nsqrt2)+1;            
        }

        if ( MakeDistanceMap_ ( F , Dst , & outlet_i , & n1_i , & nsqrt2_i ) )
        {
            SetErreur("MakeDistanceMap_ : recursive call failed\n");
            return(-1);
        }
    }
	
    free(UPSTREAM);

    return(0);
}


int MakeDistanceMap ( int HandleFlowDir, int HandleRaster, int outlet )
{
    FlowDirGrid * F;
    FloatRaster * Dst;
    int root ;

    int n1,nsqrt2;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("MakeDistanceMap : flow direction grid %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    if( !(Dst = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("MakeDistanceMap : raster %d does not exist\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(Dst->Header) ) )
    {
        SetErreur("MakeDistanceMap : raster # %d and flow direction grid # %d do not have identical referencing headers\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    n1 = 0;
    nsqrt2 = 0;
    root = outlet ;

    if ( MakeDistanceMap_ ( F , Dst , & root , & n1 , & nsqrt2 ) )
    {
        SetErreur("MakeDistanceMap : recursive call failed\n");
	return(-1);
    }

    return(0);
}


int PutRasterValue ( const int HandleRaster , const int ncells , int * CELLS , const float * VALUES )
{
    int i,nr,nc;
    FloatRaster * R;

    if ( !(R = GetFloatRaster(HandleRaster)))
    {
        SetErreur("PutRasterValue : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    for (i=0;i<ncells;i++)
    {
        (void) rec2RC ( &(R->Header) , CELLS + i , & nr, & nc );
        R->Data[nr][nc] = VALUES[i];
    }

    return(0);
}

void PutRasterValue1 ( FloatRaster * R , int * _rec , float * _f )
{
    int nr,nc;
    (void) rec2RC ( &(R->Header) ,_rec , & nr, & nc );
    R->Data[nr][nc] = *_f ;
 }


int MakeAccumulationMap ( const int HandleFlowDir, const int HandleRaster )
{
    FlowDirGrid * F;
    FloatRaster * Acc;
    int ncells,nr,nc,nup,outlet;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("MakeAccumulationMap : flow direction grid # %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    if( !(Acc = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("MakeAccumulationMap : raster # %d does not exist\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(Acc->Header) ) )
    {
        SetErreur("MakeAccumulationMap : raster # %d and flow direction grid # %d do not have identical referencing headers\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    ncells = (Acc->Header.ncols)*(Acc->Header.nrows);

    if ( SetRasterToNODATA ( HandleRaster ) )
    {
        SetErreur("MakeAccumulationMap : could not initialize raster # %d to NODATA\n",HandleRaster);
        return(-1);
    }

    for (outlet=0;outlet<ncells;outlet++)
    {
        (void) rec2RC ( &(F->Header) , & outlet , & nr , & nc );
        
        if ( ( Acc->Data[nr][nc] ) == ( Acc->Header.NODATA_value ) )   // Si on n'a pas encore calculé la surface amont
        {
            if ( MakeAccumulationMap_ ( F , Acc , & outlet , & nup , & nr , & nc ) )
            {
                SetErreur("MakeAccumulationMap : recursive call failed in cell %d\n",outlet);
                return(-1);
            }
        }
    }
	
    return(0);
}

int MakeAccumulationMap_ ( FlowDirGrid * F , FloatRaster * Acc, int * _outlet , int * _nup , int * _nr_tmp , int * _nc_tmp )
{
    int nup_i , i , nup1 ;
    int * UPSTREAM;

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("MakeAccumulationMap_ : could not get upstream cells for cell %d\n",*_outlet);
            return(-1);
    }

	*_nup = 1;
	
    for(i=0;i<nup1;i++)
    {
        if ( MakeAccumulationMap_ ( F , Acc , UPSTREAM + i , & nup_i , _nr_tmp , _nc_tmp ) )
        {
            SetErreur("MakeAccumulationMap_ : recursive call failed\n");
            return(-1);
        }

		(*_nup) += nup_i ;
    }

	free(UPSTREAM);
    
	(void) rec2RC ( &(F->Header) , _outlet , _nr_tmp , _nc_tmp );
	Acc->Data[*_nr_tmp][*_nc_tmp] = (float) (*_nup) ;
	
    return(0);
}

int RasterCDF ( const int HandleRaster , const int HandleMask , const int nf , const float * FREQ , float ** _FVALUES , int * _N )
{
    float nodata_R , nodata_M , r , m ,w0,w1,r0,r1,ix_f;
    int nr,nc,nrows,ncols,i,N,ix0;
    FloatRaster * R , * M ;
    multiset<float> Sample;
    multiset<float>::iterator imsf;

    if( !(R = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("RasterCDF : raster # %d does not exist\n",HandleRaster);
	return(-1);
    }

    if( !(M = GetFloatRaster(HandleMask) ) )
    {
        SetErreur("RasterCDF : raster # %d does not exist\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(R->Header) , &(M->Header) ) )
    {
        SetErreur("RasterCDF : rasters # %d and # %d do not have identical referencing headers\n",HandleRaster, HandleMask);
        return(-1);
    }

    nodata_R = (R->Header).NODATA_value;
    nodata_M = (M->Header).NODATA_value;
    nrows = (R->Header).nrows;
    ncols = (R->Header).ncols;

    // Récupération de l'ensemble des valeurs

    Sample.clear();

    for (nr=0;nr<nrows;nr++)
    {
        for (nc=0;nc<ncols;nc++)
        {
            r = R->Data[nr][nc];
            m = M->Data[nr][nc];

            if ( (r!=nodata_R) & (m!=nodata_M) )
                Sample.insert(r);
        }
    }

    // Calcul des quantiles

    N = Sample.size();
    (*_N) = N;
    
	if (N==0)
	{
		for(i=0;i<nf;i++)
			(*_FVALUES)[i] = nodata_R;
		
		return(0);
	}
	
    imsf = Sample.end();imsf--;
    r = (*imsf);
    Sample.insert(r); // On répète le max pour pouvoir garder la C.L. q = w0*r0 + w1*r1 même pour f=1
	
    for(i=0;i<nf;i++)
    {
        ix_f = ((float)(N-1))*FREQ[i];
        ix0 = (int) floor(ix_f);
        w1 = ix_f - floor(ix_f);
        w0 = 1.0 - w1;

        imsf = Sample.begin();
        advance ( imsf , ix0);
        r0 = (*imsf);
        imsf++;
        r1 = (*imsf);

        (*_FVALUES)[i] = w0*r0 + w1*r1;
    }

    Sample.clear();

    return(0);

}

int RasterMSE ( const int HandleRaster1 , const int HandleRaster2 , const int HandleMask , float * _MSE )
{
    // Calcul de l'erreur quadratique moyenne entre deux rasters (avec possibilité d'appliquer un masque)

    float nodata_R1 , nodata_R2 , nodata_M , r1 , r2 , m , SSE ;
    int nr , nc , nrows , ncols , N ;
    FloatRaster * R1 , * R2 , * M ;

    if( !(R1 = GetFloatRaster(HandleRaster1) ) )
    {
        SetErreur("RasterCDF : raster # %d does not exist\n",HandleRaster1);
	return(-1);
    }

    if( !(R2 = GetFloatRaster(HandleRaster2) ) )
    {
        SetErreur("RasterCDF : raster # %d does not exist\n",HandleRaster2);
	return(-1);
    }


    if( !(M = GetFloatRaster(HandleMask) ) )
    {
        SetErreur("RasterCDF : raster # %d does not exist\n",HandleMask);
	return(-1);
    }

    if ( CheckSameHeader( &(R1->Header) , &(R2->Header) ) )
    {
        SetErreur("RasterCDF : rasters # %d et # %d do not have identical referencing headers\n",HandleRaster1, HandleRaster2);
        return(-1);
    }

    if ( CheckSameHeader( &(R1->Header) , &(M->Header) ) )
    {
        SetErreur("RasterCDF : rasters # %d et # %d do not have identical referencing headers\n",HandleRaster1, HandleMask);
        return(-1);
    }

    nodata_R1 = (R1->Header).NODATA_value;
    nodata_R2 = (R2->Header).NODATA_value;
    nodata_M = (M->Header).NODATA_value;
    nrows = (R1->Header).nrows;
    ncols = (R1->Header).ncols;

    // Récupération de l'ensemble des valeurs

    SSE = 0;
    N = 0;

    for (nr=0;nr<nrows;nr++)
    {
        for (nc=0;nc<ncols;nc++)
        {
            r1 = R1->Data[nr][nc];
            r2 = R2->Data[nr][nc];
            m = M->Data[nr][nc];

            if ( (r1!=nodata_R1) & (r2!=nodata_R2) & (m!=nodata_M) )
            {
                SSE+=(r1-r2)*(r1-r2);
                N++;
            }
        }
    }

    (*_MSE) = SSE / ((float) N);

    return(0);

}

int GetRasterRowAddr ( const int HandleRaster , const int row , int * _ncols , float ** _ptrRow )
{
    FloatRaster * R;

    if ( !(R = GetFloatRaster(HandleRaster)))
    {
        SetErreur("PutRasterRow : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    if ( row<0 | row>=(R->Header.nrows) )
    {
        SetErreur("PutRasterRow : invalid row index (%d)\n",row);
        return(-1);
    }

    (*_ptrRow) = R->Data[row] ;
    (*_ncols) = R->Header.ncols;

    return(0);
}


int PutRasterRow ( const int HandleRaster , const int row , const int nv , const float * VALUES )
{
    int i;
    FloatRaster * R;

    if ( !(R = GetFloatRaster(HandleRaster)))
    {
        SetErreur("PutRasterRow : raster # %d does not exist\n",HandleRaster);
        return(-1);
    }

    if ( row<0 | row>=(R->Header.nrows) )
    {
        SetErreur("PutRasterRow : invalid row index (%d)\n",row);
        return(-1);
    }

    if ( nv!=(R->Header.ncols) )
    {
        SetErreur("PutRasterRow : the number of values (%d) is different from the number of columns in the raster (%d)\n",nv,R->Header.ncols);
        return(-1);
    }

    for (i=0;i<nv;i++)
        R->Data[row][i] = VALUES[i];

    return(0);
}

int FlagCatchmentCells_ ( FloatRaster * R , FlowDirGrid * F , int * _root , float * _flag , set<int> & Nodes , int * _nr_tmp , int * _nc_tmp , set<int> & UpNodes ,  int * _nup , int ** _BB )
{
    int i,nup,nup_i ;
    int * UPSTREAM;

    if ( GetUpstream ( F , _root , & nup , & UPSTREAM ) )
    {
            SetErreur("FlagCatchmentCells_ : could not get upstream cells for cell %d\n",*_root);
            return(-1);
    }

    *_nup = 1;

    for(i=0;i<nup;i++)
    {
        if ( Nodes.find(UPSTREAM[i]) == Nodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouvé dans la liste des noeuds d'arrêt :
        {
            if ( FlagCatchmentCells_ ( R , F , UPSTREAM + i , _flag , Nodes , _nr_tmp , _nc_tmp , UpNodes , & nup_i , _BB ) )         // Appel récursif
            {
                SetErreur("FlagCatchmentCells_ : recursive call failed\n");
                return(-1);
            }

            (*_nup) += nup_i ;
        }
        else
        {
            UpNodes.insert(UPSTREAM[i]);
        }
    }

    free(UPSTREAM);

    (void) rec2RC ( &(F->Header) , _root , _nr_tmp , _nc_tmp );
    R->Data[*_nr_tmp][*_nc_tmp] = (*_flag) ;

    /* Mise à jour de la Bounding Box */

    if ( (*_nc_tmp) < (*_BB)[0] )
        (*_BB)[0] = (*_nc_tmp);     // col min

    if ( (*_nr_tmp) < (*_BB)[1] )
        (*_BB)[1] = (*_nr_tmp);     // row min

    if ( (*_nc_tmp) > (*_BB)[2] )
        (*_BB)[2] = (*_nc_tmp);     // col max

    if ( (*_nr_tmp) > (*_BB)[3] )
        (*_BB)[3] = (*_nr_tmp);     // row max

    return(0);
}

int FlagCatchmentCells ( int HandleRaster , int HandleFlowDir , int root , float flag , int nnodes, int * NODES , int ** _BB , int * _nup , int ** _IXUPNODES , int * _nixupnodes )
{
    FlowDirGrid * F;
    FloatRaster * R;
    int nr,nc,i;
    set<int> Nodes,UpNodes;
    map<int,int> MapNodes;
    set<int>::iterator isi;
    map<int,int>::iterator imii;
    
    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("FlagCatchmentCells : grid # %d does not exist\n",HandleFlowDir);
	return(-1);
    }

    if( !(R = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("FlagCatchmentCells : raster # %d does not exist\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("FlagCatchmentCells : raster # %d and flow direction grid # %d do not have identical referencing headers\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    /* Recopie du vecteur NODES dans le set Nodes */

    Nodes.clear();
    UpNodes.clear();
    MapNodes.clear();

    for (i=0;i<nnodes;i++)
    {
        Nodes.insert( NODES[i] );
        MapNodes.insert ( make_pair ( NODES[i] , i ));
    }

    (*_BB)[0] = F->Header.ncols - 1;
    (*_BB)[1] = F->Header.nrows - 1;
    (*_BB)[2] = 0;
    (*_BB)[3] = 0;

    /* Appel racine à FlagCatchmentCells_ */

    if ( FlagCatchmentCells_ ( R , F , & root , & flag , Nodes , & nr , & nc , UpNodes , _nup ,  _BB ) )
    {
        SetErreur("FlagCatchmentCells : initial recursvice call failed\n");
        return(-1);
    }

    // Recopie des éléments du set ixUpNodes dans le tableau IXUPNODES

    (*_IXUPNODES) = (int *) malloc ( UpNodes.size()*sizeof(int) );
    (*_nixupnodes) = UpNodes.size();

    i=0;
    for (isi=UpNodes.begin();isi!=UpNodes.end();isi++)
    {
	imii = MapNodes.find(*isi);
        (*_IXUPNODES)[i] = (*imii).second;
        i++;
    }

    // Vidage des containers

    Nodes.clear();
    UpNodes.clear();
    MapNodes.clear();

    return(0);
}

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





