/* 
* 
* 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 "GeoMorph.h"
#include "CellTree.h"
#include "Erreur.h"
#include "FlowDirGrid.h"
#include "Topology.h"
#include "esriheader.h"
#include "FloatRaster.h"
#include <map>

using namespace std;

extern CellTree ** ptr_CellTree = NULL;
extern int MaxTree = 0;

typedef struct BuildCellTreeToken {

	int id ;
	int parentoffset ;
	CellTree * Tree ;

} BuildCellTreeToken ;



int CreateCellTree_iter ( int HandleFlowDir, int outlet, int HandleWeight )
{
	FlowDirGrid * F;
	FloatRaster * W;
	int HandleTree ;
	int i,nup,stacksize,parentoffset,currentoffset,prevoffset,nchildren ;
	int push ;
	int nr,nc;
	int * UPSTREAM ;
	CellTree ** ptr_CellTree_tmp;
	CellTree * TmpTree ;
	multimap<float,int> Children;
	multimap<float,int>::reverse_iterator rimmfi ;
	BuildCellTreeToken * Stack , * Stack_tmp ;

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

	if (HandleWeight>=0)
	{
		if ( !(W = GetFloatRaster(HandleWeight)) )
		{
			SetErreur("CreateCellTree_iter : float raster %d does not exist\n",HandleWeight);
			return(-1);
		}
	}
	else
		W = NULL;

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

	ClearErreur();

	/* Recherche d'un Handle libre ou crÃ©ation d'un nouveau */

	HandleTree = MaxTree - 1;

	while (HandleTree >= 0)
	{
		if(ptr_CellTree[HandleTree] == NULL) break;	/* Handle disponible */
		HandleTree--;
	}

	if(HandleTree < 0)							/* Pas de handle disponible */
	{
		if (MaxTree == 0)
		{
			ptr_CellTree = (CellTree **) malloc((MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree == NULL)
			{
				SetErreur("CreateCellTree_iter - could not allocate handle\n");
				return(-1);
			}
		}
		else
		{
			ptr_CellTree_tmp = (CellTree **) realloc(ptr_CellTree, (MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree_tmp == NULL)
			{
				SetErreur("CreateCellTree_iter - could not allocate handle\n");
				return(-1);
			}
                        else ptr_CellTree = ptr_CellTree_tmp;
		}
		HandleTree = MaxTree;
		MaxTree = MaxTree + 1;
	}

	// Initialize stack

	if ( !(Stack = (BuildCellTreeToken *) malloc ( sizeof(BuildCellTreeToken) )) )
	{
		SetErreur("CreateCellTree_iter - could not allocate stack\n");
		return(-1);
	}
	
	stacksize = 1;
	Stack[0].id = outlet ;
	Stack[0].parentoffset = 0 ;

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

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

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

				Stack_tmp = (BuildCellTreeToken *) realloc ( Stack , (stacksize+nup)*sizeof(BuildCellTreeToken) );
				if (Stack_tmp==NULL)
				{
					SetErreur("CreateCellTree_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];
					
					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 = 1 allocation

		TmpTree = (CellTree *) malloc ( sizeof(CellTree) ) ;
		if ( ! TmpTree )
		{
			SetErreur("CreateCellTree_iter: could not allocate token");
			return(-1);
		}
		TmpTree->RootCell = Stack[currentoffset].id ;
		TmpTree->FirstChild = NULL ;
		if (W)
		{
			(void) rec2RC ( &(W->Header) , & (Stack[currentoffset].id) , & nr , & nc );
			TmpTree->Size = W->Data[nr][nc];
		}
		else
		{
			TmpTree->Size = 1.0 ;
		}
		Stack[currentoffset].Tree = TmpTree ;

		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)

			TmpTree = (CellTree *) malloc ( sizeof(CellTree) ) ;
			if ( ! TmpTree )
			{
				SetErreur("CreateCellTree_iter: could not allocate token");
				return(-1);
			}
			TmpTree->RootCell = Stack[parentoffset].id ;
			if (W)
			{
				(void) rec2RC ( &(W->Header) , & (Stack[parentoffset].id) , & nr , & nc );
				TmpTree->Size = W->Data[nr][nc];
			}
			else
			{
				TmpTree->Size = 1.0 ;
			}

			Children.clear();
			for (i=0;i<nchildren;i++)
			{
				TmpTree->Size += (Stack[currentoffset+i].Tree)->Size ;
				Children.insert(make_pair( (Stack[currentoffset+i].Tree)->Size ,currentoffset+i));
			}

			rimmfi = Children.rbegin();
			TmpTree->FirstChild = Stack[(*rimmfi).second].Tree ;
			TmpTree->NextSibling = NULL;
			Stack[parentoffset].Tree = TmpTree ;

			prevoffset = (*rimmfi).second;
			rimmfi++;

			while(rimmfi!=Children.rend())
			{
				TmpTree = Stack[prevoffset].Tree ;
				TmpTree->NextSibling = Stack[(*rimmfi).second].Tree ;
				prevoffset = (*rimmfi).second;
				rimmfi++;
			}
			TmpTree = Stack[prevoffset].Tree ;
			TmpTree->NextSibling = NULL ;

			// Remove all tokens with offset >= currentoffset from the stack ;
					
			Stack_tmp = (BuildCellTreeToken *) realloc ( Stack , currentoffset * sizeof(BuildCellTreeToken) );
			if (Stack_tmp==NULL)
			{
				SetErreur("CreateCellTree_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)
		{
			ptr_CellTree[HandleTree] = Stack[0].Tree ;
			ptr_CellTree[HandleTree]->NextSibling = NULL ;
			stacksize = 0;
		}
	
	}
	
	free(Stack);
	free(UPSTREAM);
	return(HandleTree);
}



int BuildCellTree_ ( FlowDirGrid * F , int * _RootCell , CellTree * _Tree , multimap<float,int> & MapUpstream , multimap<float,int>::reverse_iterator & rimmfi )
{
    int i,nup ;
    int * UPSTREAM ;
    CellTree ** _UpstreamTrees;

    _Tree->RootCell = (*_RootCell);
    _Tree->Size = 1.0;
    _Tree->NextSibling = NULL;
    _Tree->FirstChild = NULL;

    // Récupération des indices des cellules amont

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

    if (!nup)
    {
        free(UPSTREAM);
        return(0);
    }

    if ( !( _UpstreamTrees = (CellTree **) malloc( nup * sizeof(CellTree *)) ) )
    {
            SetErreur("BuildCellTree_ : could not allocate upstream trees\n");
            return(-1);
    }

    for (i=0;i<nup;i++)
    {
        // Allocation d'un nouvel arbre

        if ( !( _UpstreamTrees[i] = (CellTree *) malloc( sizeof(CellTree)) ) )
        {
            SetErreur("BuildCellTree_ : could not allocate new tree\n");
            return(-1);
        }

        (_UpstreamTrees[i])->FirstChild = NULL;
        (_UpstreamTrees[i])->NextSibling = NULL;

        if ( BuildCellTree_ ( F , UPSTREAM + i , _UpstreamTrees[i] , MapUpstream , rimmfi ) )
        {
            SetErreur("BuildCellTree_ : recursive construction failed\n");
            FreeCellTree_ ( _UpstreamTrees[i] );
            return(-1);
        }

        // Affectation du parent :

//        (_UpstreamTrees[i])->Parent = _Tree;
    }

    MapUpstream.clear();

    for (i=0;i<nup;i++)     // Classement des arbres fils par taille croissante
    {
        MapUpstream.insert( make_pair ( (_UpstreamTrees[i])->Size , i ) );
        (_Tree->Size) += (_UpstreamTrees[i])->Size ;
    }

    i=0;
    for (rimmfi=MapUpstream.rbegin();rimmfi!=MapUpstream.rend();rimmfi++)
    {
        UPSTREAM[i] = (*rimmfi).second; // UPSTREAM[i] contient le numéro du (i+1)ème plus grand arbre à l'amont de RootCell
        i++;
    }

    // Affectation du premier fils : bassin amont de plus grande taille

    _Tree->FirstChild = _UpstreamTrees[UPSTREAM[0]];

    // Affectation des frères aux autres bassins amont

    for (i=0;i<(nup-1);i++)
        (_UpstreamTrees[UPSTREAM[i]])->NextSibling = _UpstreamTrees[UPSTREAM[i+1]] ;

    (_UpstreamTrees[UPSTREAM[nup-1]])->NextSibling = NULL ;

    free(_UpstreamTrees);
    free(UPSTREAM);

    return(0);

}

int CreateCellTree ( const int HandleFlowDir , int RootCell )
{
    FlowDirGrid * F;
    int HandleTree = MaxTree;
    int * _RootCell;
    CellTree ** ptr_CellTree_tmp;
    CellTree * T;
    multimap<float,int> MapUpstream;
    multimap<float,int>::reverse_iterator rimmfi;


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

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

	ClearErreur();

	/* Recherche d'un Handle libre ou crÃ©ation d'un nouveau */

	HandleTree = MaxTree - 1;

	while (HandleTree >= 0)
	{
		if(ptr_CellTree[HandleTree] == NULL) break;	/* Handle disponible */
		HandleTree--;
	}

	if(HandleTree < 0)							/* Pas de handle disponible */
	{
		if (MaxTree == 0)
		{
			ptr_CellTree = (CellTree **) malloc((MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree == NULL)
			{
				SetErreur("CreateCellTree - could not allocate handle\n");
				return(-1);
			}
		}
		else
		{
			ptr_CellTree_tmp = (CellTree **) realloc(ptr_CellTree, (MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree_tmp == NULL)
			{
				SetErreur("CreateCellTree - could not allocate handle\n");
				return(-1);
			}
                        else ptr_CellTree = ptr_CellTree_tmp;
		}
		HandleTree = MaxTree;
		MaxTree = MaxTree + 1;
	}

	/* Creation d'un nouvel arbre rattaché au handle */

	if (!(ptr_CellTree[HandleTree] = (CellTree *) malloc(sizeof(CellTree))))
	{
		SetErreur("CreateCellTree - could not allocate tree structure\n");
		return(-1);
	}

	T = ptr_CellTree[HandleTree];
//        T->Parent = NULL;

        // Construction récursive de l'arbre

//        (*_RootCell) = RootCell ;

        if ( BuildCellTree_ ( F , & RootCell ,  T , MapUpstream , rimmfi ) )
        {
             SetErreur("CreateCellTree - recursive construction failed\n");
             return(-1);
        }

        return(HandleTree);
}


int FreeCellTree_ ( CellTree * T )
{
    // Désallocation récursive d'un arbre

    if ( T->FirstChild )
    {
       if( FreeCellTree_ (T->FirstChild) )
       {
            SetErreur("DestroyCellTree_ : recursive deallocation failed\n");
            return(-1);
       }
    }

    if ( T->NextSibling )
    {
       if( FreeCellTree_ (T->NextSibling ) )
       {
            SetErreur("DestroyCellTree_ : recursive deallocation failed\n");
            return(-1);
       }
    }

    free(T);
    return(0);

}


int DestroyCellTree ( const int HandleTree )
{
    CellTree * T;

    if ( !(T = GetCellTree(HandleTree)) )
    {
        SetErreur("DestroyCellTree : tree %d does not exist\n",HandleTree);
        return(-1);
    }

    if ( FreeCellTree_ (ptr_CellTree[HandleTree]) )
    {
        SetErreur("DestroyCellTree (%d) : deallocation failed\n",HandleTree);
        return(-1);
    }

    ptr_CellTree[HandleTree] = NULL;

    return(0);
}

CellTree * GetCellTree (const int HandleTree)
{
		if ((HandleTree < 0) | (HandleTree >= MaxTree)) return(NULL);
                else return(ptr_CellTree[HandleTree]);
}

int MakeStreamLines ( const int HandleTree , float * _threshold , int * _todivide , int * _nstreams , int ** _StreamLength , int *** _StreamConnect , int *** _StreamCells )
{
    CellTree * T;

    if ( !(T = GetCellTree( HandleTree )) )
    {
        SetErreur("MakeStreamLines : tree %d does not exist\n",HandleTree);
        return(-1);
    }

    if ( T->Size < (*_threshold) )
    {
        SetErreur("MakeStreamLines : total tree size (%d) is smaller than the specified extraction threshold (%d)\n",T->Size,(*_threshold));
        return(-1);
    }

    (*_nstreams) = 1;

    (*_StreamLength) =  (int * ) malloc ( sizeof(int) );
    (*_StreamLength)[(*_nstreams)-1] = 0;

    (*_StreamConnect)[0] = (int * ) malloc ( sizeof(int) );
    (*_StreamConnect)[1] = (int * ) malloc ( sizeof(int) );
    (*_StreamConnect)[0][(*_nstreams)-1] = 1;
    (*_StreamConnect)[1][(*_nstreams)-1] = T->RootCell;

    (*_StreamCells) = (int **) malloc ( sizeof (int *) );

    if ( MakeStreamLines_ ( T , _threshold , _todivide , T , _nstreams , _StreamLength , _StreamConnect , _StreamCells ) )
    {
        return(-1);
    }

    return(0);

}



int MakeStreamLines_ ( CellTree * T , float * _threshold , int * _todivide , CellTree * Ancestor , int * _nstreams , int ** _StreamLength ,  int *** _StreamConnect , int *** _StreamCells )
{
    int len,n_tmp;

    /* On récupère la longueur de la ligne courante */

    len = (*_StreamLength)[(*_nstreams)-1] ;

    if (len==0)
    {
        (*_StreamCells)[(*_nstreams)-1] = (int * ) malloc ( sizeof(int) );
    }
    else
    {
        (*_StreamCells)[(*_nstreams)-1] = (int * ) realloc ( (*_StreamCells)[(*_nstreams)-1] , (len+1)*sizeof(int) );
    }

    (*_StreamCells)[(*_nstreams)-1][len] = T->RootCell ;
    (*_StreamLength)[(*_nstreams)-1] += 1;
    n_tmp = (*_nstreams);

    if ( T->FirstChild )    // Si la cellule a un BV amont, on continue à construire la ligne courante
    {
        if ( 	(*_todivide) 	|	 ( (T->FirstChild)->Size > (*_threshold) )	 )  // Si le BV amont est supérieur au seuil ou qu'on extrait jusqu'à la ligne de crête
        {
            if ( MakeStreamLines_ ( T->FirstChild , _threshold , _todivide , T , _nstreams , _StreamLength , _StreamConnect , _StreamCells ) )
            {
                return(-1);
            }
        }
    }

    if ( T->NextSibling )
    {
       if ( (T->NextSibling)->Size > (*_threshold) ) // Si le BV est supérieur au seuil
        {
            (*_nstreams)+=1;

            (*_StreamLength) =  (int * ) realloc ( (*_StreamLength) , (*_nstreams)*sizeof(int) );
            (*_StreamLength)[(*_nstreams)-1] = 0;

            (*_StreamConnect)[0] = (int * ) realloc ( (*_StreamConnect)[0] , (*_nstreams)*sizeof(int) );
            (*_StreamConnect)[1] = (int * ) realloc ( (*_StreamConnect)[1] , (*_nstreams)*sizeof(int) );
            (*_StreamConnect)[0][(*_nstreams)-1] = n_tmp;
            (*_StreamConnect)[1][(*_nstreams)-1] = Ancestor->RootCell ; // T->RootCell; BUG ?

            (*_StreamCells) = (int ** ) realloc ( (*_StreamCells) , (*_nstreams)*sizeof(int *) );

            if ( MakeStreamLines_ ( T->NextSibling , _threshold , _todivide , Ancestor , _nstreams , _StreamLength , _StreamConnect , _StreamCells ) )
            {
                return(-1);
            }
       }
    }

    return(0);

}

int SampleMainStream_ ( CellTree * T , int * _ntribut , int ** _TributOutlet , int ** _TributSize )
{

    CellTree ** _Ts;

    // On collecte toute la fratrie de l'arbre T (i.e., les arbres confluant avec le drain principal)

    _Ts = & T ;

    while ( (*_Ts)->NextSibling )
    {
        (*_ntribut)+=1;

        (*_TributOutlet) =  (int * )  realloc ( (*_TributOutlet) , (*_ntribut)*sizeof(int) );
        (*_TributOutlet)[(*_ntribut)-1] = ((*_Ts)->NextSibling)->RootCell ;

        (*_TributSize) =  (int * )  realloc ( (*_TributSize) , (*_ntribut)*sizeof(int) );
        (*_TributSize)[(*_ntribut)-1] = ((*_Ts)->NextSibling)->Size ;

        _Ts = & ((*_Ts)->NextSibling);
    }

    // Si l'arbre T a un fils, on continue récursivement sur le drain principal

    if ( T->FirstChild )
        SampleMainStream_ ( T->FirstChild , _ntribut , _TributOutlet , _TributSize );

    return(0);
}

int SampleMainStream ( const int HandleTree , int * _ntribut , int ** _TributOutlet , int ** _TributSize )
{

    CellTree * T;

    if ( !(T = GetCellTree( HandleTree )) )
    {
        SetErreur("SampleMainStream : tree %d does not exist\n",HandleTree);
        return(-1);
    }

    (*_ntribut) = 0;
    (*_TributOutlet) =  NULL ;
    (*_TributSize) =  NULL ;

    if ( T->FirstChild )
        SampleMainStream_ ( T->FirstChild , _ntribut , _TributOutlet , _TributSize );

    return(0);
}

int MainBifurcationNode ( int HandleTree , int * _A , int * _B , int * _C , int * _D )
{
    /* Identifie le point de bifurcation principal (Moussa, 2008) D
     * du bassin dÃ©fini par l'exutoire O */

    CellTree * T;
    float SB;

    if ( !(T = GetCellTree( HandleTree )) )
    {
        SetErreur("MainBifurcationNode : tree %d does not exist\n",HandleTree);
		return(-1);
    }

    SB = 0.0;
    *_A = -1;
    *_B = -1;
    *_C = -1;
    *_D = -1;

    if ( MainBifurcationNode_ ( T , _A , _B , _C , _D , & SB ) )
    {
        SetErreur("MainBifurcationNode : recursive call failed\n");
        return(-1);
    }

    return(0);

}

int MainBifurcationNode_ ( CellTree * T , int * _A , int * _B , int * _C , int * _D , float * _SB )
{
	CellTree * Tc;
	
	Tc = T->FirstChild;
	
	if ( ! Tc )
		return(0);

    if ( Tc->NextSibling )
	{
		if ( (Tc->NextSibling)->Size > (*_SB) )
		{
			(*_SB) = (Tc->NextSibling)->Size ;
			(*_B) = (Tc->NextSibling)->RootCell ;
			(*_C) = Tc->RootCell ;
			(*_D) = T->RootCell ;
		}
	}
		
    if ( Tc->Size > (*_SB) )
    {
        if ( MainBifurcationNode_ ( T->FirstChild , _A , _B , _C , _D , _SB ) )
		{
			SetErreur("MainBifurcationNode_ : recursive call failed\n");
			return(-1);
		}
    }
    else
    {
        (*_A) = T->RootCell ;
    }

    return(0);
}

int UpstreamAverage_ ( CellTree * T , FloatRaster * R , FloatRaster * Rm, float * _sum )
{
    int nr,nc;
    float sumc;
    CellTree * Tc ;
    CellTree * Tc2 ;

    (void) rec2RC ( &(R->Header) , &(T->RootCell)  , & nr , & nc );
    (*_sum) = R->Data[nr][nc];

    if ( T->FirstChild )
    {
        if ( UpstreamAverage_ ( T->FirstChild , R , Rm , & sumc ) )
        {
            SetErreur("UpstreamAverage_ : recursive call failed for first child tree\n");
            return(-1);
        }

        (*_sum) = (*_sum) + sumc ;

        Tc = T->FirstChild;

        while ( Tc->NextSibling )   // si le fils de T a un petit frère
        {
            Tc2 = Tc->NextSibling;   // le frère courant est le suivant

            if ( UpstreamAverage_ ( Tc2 , R , Rm , & sumc ) )
            {
                SetErreur("UpstreamAverage_ : recursive call failed for next sibling tree\n");
                return(-1);
            }

            (*_sum) = (*_sum) + sumc ;
            Tc = Tc2;
        }

    }

    Rm->Data[nr][nc] = (*_sum) / ((float)(T->Size)) ;    // moyenne = somme / taille de l'arbre

    return(0);
}


int UpstreamAverage ( int HandleTree , int HandleRasterLocal , int HandleRasterAveraged )
{
    CellTree * T;
    FloatRaster * R;
    FloatRaster * Rm;
    float sum;

    if( !(T = GetCellTree(HandleTree) ) )
    {
        SetErreur("UpstreamAverage : tree # %d does not exist\n",HandleTree);
	return(-1);
    }

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

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

    if ( UpstreamAverage_ ( T , R , Rm , & sum ) )
    {
        SetErreur("UpstreamAverage : recursive call failed at root cell\n");
        return(-1);
    }

    return(0);
}


int BuildWeightedCellTree_ ( FlowDirGrid * F , int * _RootCell , FloatRaster * W , CellTree * _Tree , multimap<float,int> & MapUpstream , multimap<float,int>::reverse_iterator & rimmfi , int * _nr_tmp , int * _nc_tmp )
{
    int i,nup ;
    int * UPSTREAM ;
    CellTree ** _UpstreamTrees;

    _Tree->RootCell = (*_RootCell);
	
	// Get rootcell weight
	(void) rec2RC ( &(W->Header) , _RootCell , _nr_tmp , _nc_tmp );
    _Tree->Size = W->Data[*_nr_tmp][*_nc_tmp];
	
    _Tree->NextSibling = NULL;
    _Tree->FirstChild = NULL;

    // Récupération des indices des cellules amont

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

    if (!nup)
    {
        free(UPSTREAM);
        return(0);
    }

    if ( !( _UpstreamTrees = (CellTree **) malloc( nup * sizeof(CellTree *)) ) )
    {
            SetErreur("BuildWeightedCellTree_ : could not allocate upstream trees\n");
            return(-1);
    }

    for (i=0;i<nup;i++)
    {
        // Allocation d'un nouvel arbre

        if ( !( _UpstreamTrees[i] = (CellTree *) malloc( sizeof(CellTree)) ) )
        {
            SetErreur("BuildWeightedCellTree_ : could not allocate new tree\n");
            return(-1);
        }

        (_UpstreamTrees[i])->FirstChild = NULL;
        (_UpstreamTrees[i])->NextSibling = NULL;

        if ( BuildWeightedCellTree_ ( F , UPSTREAM + i , W , _UpstreamTrees[i] , MapUpstream , rimmfi , _nr_tmp , _nc_tmp ) )
        {
            SetErreur("BuildWeightedCellTree_ : recursive construction failed\n");
            FreeCellTree_ ( _UpstreamTrees[i] );
            return(-1);
        }

        // Affectation du parent :

//        (_UpstreamTrees[i])->Parent = _Tree;
    }

    MapUpstream.clear();

    for (i=0;i<nup;i++)     // Classement des arbres fils par taille croissante
    {
        MapUpstream.insert( make_pair ( (_UpstreamTrees[i])->Size , i ) );
        (_Tree->Size) += (_UpstreamTrees[i])->Size ;
    }

    i=0;
    for (rimmfi=MapUpstream.rbegin();rimmfi!=MapUpstream.rend();rimmfi++)
    {
        UPSTREAM[i] = (*rimmfi).second; // UPSTREAM[i] contient le numéro du (i+1)ème plus grand arbre à l'amont de RootCell
        i++;
    }

    // Affectation du premier fils : bassin amont de plus grande taille

    _Tree->FirstChild = _UpstreamTrees[UPSTREAM[0]];

    // Affectation des frères aux autres bassins amont

    for (i=0;i<(nup-1);i++)
        (_UpstreamTrees[UPSTREAM[i]])->NextSibling = _UpstreamTrees[UPSTREAM[i+1]] ;

    (_UpstreamTrees[UPSTREAM[nup-1]])->NextSibling = NULL ;

    free(_UpstreamTrees);
    free(UPSTREAM);

    return(0);

}

int CreateWeightedCellTree ( const int HandleFlowDir , int RootCell , int HandleWeight )
{
    FlowDirGrid * F;
	FloatRaster * W;
    int HandleTree = MaxTree;
    int * _RootCell;
    CellTree ** ptr_CellTree_tmp;
    CellTree * T;
    multimap<float,int> MapUpstream;
    multimap<float,int>::reverse_iterator rimmfi;
	int nr,nc;

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

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

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

	ClearErreur();

	/* Recherche d'un Handle libre ou crÃ©ation d'un nouveau */

	HandleTree = MaxTree - 1;

	while (HandleTree >= 0)
	{
		if(ptr_CellTree[HandleTree] == NULL) break;	/* Handle disponible */
		HandleTree--;
	}

	if(HandleTree < 0)							/* Pas de handle disponible */
	{
		if (MaxTree == 0)
		{
			ptr_CellTree = (CellTree **) malloc((MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree == NULL)
			{
				SetErreur("CreateCellTree - could not allocate handle\n");
				return(-1);
			}
		}
		else
		{
			ptr_CellTree_tmp = (CellTree **) realloc(ptr_CellTree, (MaxTree + 1) * sizeof(CellTree *));
			if (ptr_CellTree_tmp == NULL)
			{
				SetErreur("CreateCellTree - could not allocate handle\n");
				return(-1);
			}
                        else ptr_CellTree = ptr_CellTree_tmp;
		}
		HandleTree = MaxTree;
		MaxTree = MaxTree + 1;
	}

	/* Creation d'un nouvel arbre rattaché au handle */

	if (!(ptr_CellTree[HandleTree] = (CellTree *) malloc(sizeof(CellTree))))
	{
		SetErreur("CreateCellTree - could not allocate tree structure\n");
		return(-1);
	}

	T = ptr_CellTree[HandleTree];
//        T->Parent = NULL;

        // Construction récursive de l'arbre

//        (*_RootCell) = RootCell ;

        if ( BuildWeightedCellTree_ ( F , & RootCell , W , T , MapUpstream , rimmfi , & nr , & nc ) )
        {
             SetErreur("CreateCellTree - recursive construction failed\n");
             return(-1);
        }

        return(HandleTree);
}
