/* 
* 
* 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 "FloatRaster.h"
#include "FlowDirGrid.h"
#include "Topology.h"
#include "Mesh.h"
#include "MeshTree.h"
#include "esriheader.h"
#include "Erreur.h"
#include <set>
#include <map>
#include <algorithm>

using namespace std;

int SegmentationTopo_ ( FlowDirGrid * F , FloatRaster * R , int * _root , double * _cumul , set<int> & Constraints , map<int,double> & Nodes , double * _downcontrib )
{
    int nup,i,nr,nc ;
    int * UPSTREAM ;
    set<int>::iterator isi;
    double downcontrib_i;

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

    (void) rec2RC ( &(R->Header) , _root , & nr, & nc );
    (*_downcontrib) = (double) (R->Data[nr][nc]) ;     // Lecture de la valeur locale de la grandeur a cumuler

    // Parcours des cellules amont

    for(i=0;i<nup;i++)
    {
        if ( SegmentationTopo_ ( F , R , UPSTREAM + i , _cumul , Constraints , Nodes , & downcontrib_i ) )
        {
            SetErreur("SegmentationTopo_ : recursive call failed\n");
            return(-1);
        }

        (*_downcontrib) += downcontrib_i;
    }

    // Si la contribution a la cellule (*_root) depasse "cumul" ou si cette cellule est un point de contrainte,
    // elle est rajoutee a la liste de noeuds et (*_downcontrib) est mis a zero

    isi = Constraints.find ( *_root );

    if ( (*_downcontrib) >= (*_cumul) | isi != Constraints.end() )
    {
         Nodes.insert ( make_pair ( *_root , * _downcontrib ));
         (*_downcontrib) = 0;
    }

    free(UPSTREAM);

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int SegmentationTopo ( int HandleFlowDir , int HandleRaster , int exu , double cumul , int * CONSTRAIN , int nconstrain, int * _nnodes , int ** _NODES , double ** _WEIGHTS )
{
    /***************************************************************************************************************

     Decoupe le domaine topographique a l'amont du pixel "exu" en sous-bassins contenant au moins la valeur "cumul"
     d'une grandeur locale cartographiee dans le raster d'identifiant "HandleRaster". Renvoie un tableau d'entiers
     contenant les identifiants des pixels exutoires de mailles.

     EXEMPLE : si la carte de grandeur locale ne contient que des 1, alors la fonction cree un maillage dans lequel
     chaque maille contient au moins "Cumul" pixels.

     ***************************************************************************************************************/

    FloatRaster * R ;
    FlowDirGrid * F ;
    double downcontrib ;
    int i;
    map<int,double> Nodes ;
    set<int> Constraints;
    map<int,double>::iterator imid;

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

    for (i=0;i<nconstrain;i++)
    {
        Constraints.insert(CONSTRAIN[i]);
    }

    Nodes.clear();

    if ( SegmentationTopo_ ( F , R , & exu , & cumul , Constraints , Nodes , & downcontrib ) )
    {
        SetErreur("SegmentationTopo : initial recursive call failed for cell %d\n",exu);
        return(-1);
    }

    Nodes.insert ( make_pair ( exu , downcontrib ) );  /* La cellule de l'appel primaire est necesssairement un noeud */

    (*_NODES) = (int *) malloc ( Nodes.size() * sizeof(int) ) ;
    (*_WEIGHTS) = (double *) malloc ( Nodes.size() * sizeof(double) ) ;
    (*_nnodes) = Nodes.size();
    
    i=0;

    for(imid=Nodes.begin();imid!=Nodes.end();imid++)
    {
        (*_NODES)[i] = (*imid).first;
        (*_WEIGHTS)[i] = (*imid).second;
        i++;
    }

    Nodes.clear();

    return(0);

}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MakeFlowLineThroughPoint ( int HandleFlowDir , int HandleRaster , int exu , int throughpoint , int ** _FlowLine , int * _length )
{
    FlowDirGrid * F ;
    FloatRaster * R ;
    int n,Amont,Aval ;

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

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

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

    // Extraction de la partie a l'amont de throughpoint

    (*_FlowLine) = NULL;
    n = 0;
    
    if ( MakeFlowLineThroughPoint_ ( F , R , & throughpoint , _FlowLine , & n ) )
    {
        SetErreur("MakeFlowLineThroughPoint : could not extract upstream part of the flowline\n");
        return(-1);
    }

    // Extraction de la partie a l'aval de throughpoint

    Amont = throughpoint;
    Aval = GetDownstream ( F , & throughpoint );

    while ( Aval>=0 & Amont!=exu & Aval!=Amont )
    {
        n++;
        (*_FlowLine) = (int *) realloc( (*_FlowLine) , n*sizeof(int) );
        (*_FlowLine)[n-1] = Aval;
        Amont = Aval;
        Aval = GetDownstream ( F , & Amont );
    }

    (*_length) = n;

    return(0);

}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MakeFlowLineThroughPoint_ ( FlowDirGrid * F , FloatRaster * R , int * _root , int ** _FlowLine , int * _length )
{
    int nup,i,nr,nc,iupmax ;
    float f,fupmax ;
    int * UPSTREAM ;

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

    if (nup>0)
    {
        fupmax = 0 ;
        iupmax = 0;
        
        for(i=0;i<nup;i++)
        {
            (void) rec2RC ( &(R->Header) , UPSTREAM + i , & nr, & nc );
            f = R->Data[nr][nc];

            if(f>fupmax)
            {
                fupmax = f;
                iupmax = i;
            }
        }

        if ( MakeFlowLineThroughPoint_ ( F , R , UPSTREAM + iupmax , _FlowLine , _length ) )
        {
            SetErreur("MakeFlowLineThroughPoint_ : recursive call failed\n");
            return(-1);
        }

        (*_length) += 1 ;
        (*_FlowLine) = (int *) realloc ( (*_FlowLine), (*_length) *sizeof(int) );
        (*_FlowLine)[(*_length)-1] = (*_root);
    }
    else
    {
        (*_length) = 1;
        (*_FlowLine) = (int *) malloc ( sizeof(int) );
        (*_FlowLine)[0] = (*_root);
    }

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetCatchmentSum_ ( FlowDirGrid * F , FloatRaster * R , int * _outlet, set<int> & Nodes, int ** _BB , float * _Sum , set<int> & UpNodes , int * _nr , int * _nc )
{
    int i , nup1 ;
    float Sum_i ;
    set<int>::iterator isi;
    int * UPSTREAM;

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

    (void) rec2RC ( &(F->Header) , _outlet , _nr , _nc );
    (*_Sum) = (R->Data)[(*_nr)][*_nc] ;

    // Mise a jour de la Bounding Box

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

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

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

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

    // Boucle sur les pixels amont

    for(i=0;i<nup1;i++)
    {

        isi = Nodes.find(UPSTREAM[i]);

        if ( isi == Nodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouve dans la liste des noeuds d'arret
        {
            if ( GetCatchmentSum_ ( F , R , UPSTREAM + i , Nodes , _BB , & Sum_i , UpNodes , _nr , _nc ) )     // Appel recursif
            {
                SetErreur("GetCatchmentSum_ : recursive call failed\n");
                return(-1);
            }

            (*_Sum) += Sum_i ;
        }
        else
        {
            UpNodes.insert ( *isi );
        }
    }

    free(UPSTREAM);

    return (0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetCatchmentSum ( int HandleFlowDir , int HandleRaster , int outlet, int nnodes, int * NODES , int ** _BB , float * _Sum , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    FloatRaster * R;
    set<int> Nodes,UpNodes;
    map<int,int> MapNodes;
    set<int>::iterator isi;
    map<int,int>::iterator imii;
    int nr;
    int nc;

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

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

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

    // Recopie des elements du tableau 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;

    if ( GetCatchmentSum_ ( F , R , & outlet , Nodes, _BB , _Sum , UpNodes , & nr , & nc ) )
    {
        SetErreur("GetCatchmentSum : recursive call failed\n");
	return(-1);
    }

    // Recopie des elements du set UpNodes 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);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetCatchmentSum1_ ( FlowDirGrid * F , int * _outlet, set<int> & Nodes, int ** _BB , float * _Sum , set<int> & UpNodes , int * _nr , int * _nc )
{
    int i , nup1 ;
    float Sum_i ;
    set<int>::iterator isi;
    int * UPSTREAM;

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

    (void) rec2RC ( &(F->Header) , _outlet , _nr , _nc );
    (*_Sum) = 1 ;

    // Mise a jour de la Bounding Box

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

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

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

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

    // Boucle sur les pixels amont

    for(i=0;i<nup1;i++)
    {

        isi = Nodes.find(UPSTREAM[i]);

        if ( isi == Nodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouve dans la liste des noeuds d'arret
        {
            if ( GetCatchmentSum1_ ( F , UPSTREAM + i , Nodes , _BB , & Sum_i , UpNodes , _nr , _nc ) )     // Appel recursif
            {
                SetErreur("GetCatchmentSum1_ : recursive call failed\n");
                return(-1);
            }

            (*_Sum) += Sum_i ;
        }
        else
        {
            UpNodes.insert ( *isi );
        }
    }

    free(UPSTREAM);

    return (0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetCatchmentSum1 ( int HandleFlowDir , int outlet, int nnodes, int * NODES , int ** _BB , float * _Sum , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    FloatRaster * R;
    set<int> Nodes,UpNodes;
    map<int,int> MapNodes;
    set<int>::iterator isi;
    map<int,int>::iterator imii;
    int nr;
    int nc;

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

    *_Sum = 0;

    // Recopie des elements du tableau UPNODES dans le set UpNodes

    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;

    if ( GetCatchmentSum1_ ( F , & outlet , Nodes, _BB , _Sum , UpNodes , & nr , & nc ) )
    {
        SetErreur("GetCatchmentSum1 : recursive call failed\n");
	return(-1);
    }

    // Recopie des elements 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);
}
//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentMean_ ( FlowDirGrid * F , FloatRaster * Loc , FloatRaster * Acc , int * _outlet , set<int> & Inlets , float * _Sum , int * _Sum1 , int * _nr_tmp , int * _nc_tmp )
{
    int i , nup ;
    float Sum_i,f ;
    int Sum1_i ;
    int * UPSTREAM;
    set<int>::iterator isi;

    (*_Sum) = 0;
    (*_Sum1) = 1;

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

    for(i=0;i<nup;i++)
    {
	isi = Inlets.find(UPSTREAM[i]);

	if (isi==Inlets.end())
	{		
	        if ( MapCatchmentMean_ ( F , Loc , Acc , UPSTREAM + i , Inlets , & Sum_i , & Sum1_i , _nr_tmp , _nc_tmp ) )
        	{
	            SetErreur("MapCatchmentMean_ : recursive call failed\n");
        	    return(-1);
	        }

        	(*_Sum) += Sum_i ;
	        (*_Sum1) += Sum1_i ;
	}
    }

    free(UPSTREAM);

    (void) rec2RC ( &(F->Header) , _outlet , _nr_tmp , _nc_tmp );
    f = Loc->Data[*_nr_tmp][*_nc_tmp];
    (*_Sum) += f;

    Acc->Data[*_nr_tmp][*_nc_tmp] = (*_Sum)/((float) (*_Sum1));

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentMean ( int HandleFlowDir, int HandleLocalRaster , int exu , int HandleCumulRaster , int ninlets , int * INLETS )
{
    FlowDirGrid * F;
    FloatRaster * Loc , * Acc;
    int nr,nc,i;
    float Sum;
    int Sum1;
    set<int> Inlets;

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

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

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

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

    if ( CheckSameHeader( &(Loc->Header) , &(Acc->Header) ) )
    {
        SetErreur("MapCatchmentMean : raster # %d and # %d do not have identical referencing headers\n",HandleLocalRaster,HandleCumulRaster);
        return(-1);
    }

    Inlets.clear();

    for(i=0;i<ninlets;i++)
	    Inlets.insert(INLETS[i]);


    if ( MapCatchmentMean_ ( F , Loc , Acc , & exu , Inlets , & Sum , & Sum1 , & nr , & nc ) )
    {
        SetErreur("MapCatchmentMean : initial recursive call failed at cell %d\n",exu);
        return(-1);
    }

    Inlets.clear();

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentSum_ ( FlowDirGrid * F , FloatRaster * Loc , FloatRaster * Acc , int * _outlet , set<int> & Inlets , float * _Sum , int * _nr_tmp , int * _nc_tmp )
{
    int i , nup ;
    float Sum_i,f ;
    int * UPSTREAM;
    set<int>::iterator isi;

    (*_Sum) = 0;

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

    for(i=0;i<nup;i++)
    {
	isi = Inlets.find(UPSTREAM[i]);

	if (isi==Inlets.end())
	{		
	        if ( MapCatchmentSum_ ( F , Loc , Acc , UPSTREAM + i , Inlets , & Sum_i , _nr_tmp , _nc_tmp ) )
        	{
	            SetErreur("MapCatchmentSum_ : recursive call failed\n");
        	    return(-1);
	        }

        	(*_Sum) += Sum_i ;
	}
    }

    free(UPSTREAM);

    (void) rec2RC ( &(F->Header) , _outlet , _nr_tmp , _nc_tmp );
    f = Loc->Data[*_nr_tmp][*_nc_tmp];
    (*_Sum) += f;

    Acc->Data[*_nr_tmp][*_nc_tmp] = (*_Sum);

    return(0);
}
//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentSum ( int HandleFlowDir, int HandleLocalRaster , int exu , int HandleCumulRaster , int ninlets , int * INLETS )
{
    FlowDirGrid * F;
    FloatRaster * Loc , * Acc;
    int nr,nc,i;
    float Sum;
    set<int> Inlets;

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

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

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

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

    if ( CheckSameHeader( &(Loc->Header) , &(Acc->Header) ) )
    {
        SetErreur("MapCatchmentSum : raster # %d and # %d do not have identical referencing headers\n",HandleLocalRaster,HandleCumulRaster);
        return(-1);
    }

    Inlets.clear();

    for(i=0;i<ninlets;i++)
	    Inlets.insert(INLETS[i]);


    if ( MapCatchmentSum_ ( F , Loc , Acc , & exu , Inlets , & Sum , & nr , & nc ) )
    {
        SetErreur("MapCatchmentSum : initial recursive call failed at cell %d\n",exu);
        return(-1);
    }

    Inlets.clear();

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentSum1_ ( FlowDirGrid * F , FloatRaster * Acc , int * _outlet , set<int> & Inlets , int * _Sum , int * _nr_tmp , int * _nc_tmp )
{
    int i , nup ;
    int Sum_i,f ;
    int * UPSTREAM;
    set<int>::iterator isi;

    (*_Sum) = 1;

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

    for(i=0;i<nup;i++)
    {
	isi = Inlets.find(UPSTREAM[i]);

	if (isi==Inlets.end())
	{		
	        if (MapCatchmentSum1_ ( F , Acc , UPSTREAM + i , Inlets , & Sum_i , _nr_tmp , _nc_tmp ) )
        	{
	            SetErreur("MapCatchmentSum1_ : recursive call failed\n");
        	    return(-1);
	        }

        	(*_Sum) += Sum_i ;
	}
    }

    free(UPSTREAM);

    (void) rec2RC ( &(F->Header) , _outlet , _nr_tmp , _nc_tmp );
    Acc->Data[*_nr_tmp][*_nc_tmp] = (float)( (*_Sum) );

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int MapCatchmentSum1 ( int HandleFlowDir, int exu , int HandleCumulRaster , int ninlets , int * INLETS )
{
    FlowDirGrid * F;
    FloatRaster * Acc;
    int nr,nc,i;
    int Sum;
    set<int> Inlets;

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

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

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

    Inlets.clear();

    for(i=0;i<ninlets;i++)
	    Inlets.insert(INLETS[i]);


    if ( MapCatchmentSum1_ ( F , Acc , & exu , Inlets , & Sum , & nr , & nc ) )
    {
        SetErreur("MapCatchmentSum1 : initial recursive call failed at cell %d\n",exu);
        return(-1);
    }

    Inlets.clear();

    return(0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetUpstreamNodes_ ( FlowDirGrid * F , int * _outlet, set<int> & Nodes, set<int> & UpNodes )
{
    int i , nup1 ;
    set<int>::iterator isi;
    int * UPSTREAM;

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

    // Boucle sur les pixels amont

    for(i=0;i<nup1;i++)
    {

        isi = Nodes.find(UPSTREAM[i]);

        if ( isi == Nodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouve dans la liste des noeuds d'arret
        {
            if ( GetUpstreamNodes_ ( F , UPSTREAM + i , Nodes , UpNodes ) )     // Appel recursif
            {
                SetErreur("GetUpstreamNodes_ : recursive call failed\n");
                return(-1);
            }
        }
        else
        {
            UpNodes.insert ( *isi );
        }
    }

    free(UPSTREAM);

    return (0);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int GetUpstreamNodes ( int HandleFlowDir , int outlet, int nnodes, int * NODES , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    set<int> Nodes,UpNodes;
    map<int,int> MapNodes;
    set<int>::iterator isi;
    map<int,int>::iterator imii;

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

    // Recopie des elements du tableau UPNODES dans le set UpNodes

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

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

    if (  GetUpstreamNodes_ ( F , & outlet , Nodes, UpNodes ) )
    {
        SetErreur("GetUpstreamNodes : recursive call failed\n");
	return(-1);
    }

    // Recopie des elements 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);
}

//*********************************************************************************************************************************************************
//*********************************************************************************************************************************************************

int OptimMeshVariance ( int HandleFlowDir , int HandleCumulRaster , int HandlePkRaster , int nnodes , int * NODES , int * CNT , double damper , int ** _NEWNODES , double ** _WEIGHTS )
{
    int i,j;
    double meanweight ;
    MeshTree * MTree ;
    FloatRaster * CW , * PK ;
    FlowDirGrid * F ;

    (*_NEWNODES) = NULL ;
    (*_WEIGHTS) = NULL ;

    ClearErreur();

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

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

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


    if ( MakeMeshTree ( F , CW , PK , nnodes , NODES , CNT , & MTree , & meanweight ) )
    {
        SetErreur("OptimMeshVariance : could not create mesh tree\n");
        return(-1);
    }

    for(i=0;i<nnodes;i++)
    {
        if ( MTree[i].cnt | !(MTree[i].Parent) )
        {
            meanweight = (MTree[i].TreeWeight) / ((double)(1+MTree[i].ndescendants));

            if ( OptimMeshVariance_ ( MTree + i , & meanweight , & meanweight , & damper ) )
            {
                SetErreur("OptimMeshVariance : optimization failed upstream node %d\n", MTree[i].RootCell );
                return(-1);
            }
        }
    }

    (*_NEWNODES) = (int *) malloc ( nnodes * sizeof (int) );
    (*_WEIGHTS) = (double *) malloc ( nnodes * sizeof (double) );

    for(i=0;i<nnodes;i++)
    {
        (*_NEWNODES)[i] = MTree[i].RootCell ;
        (*_WEIGHTS)[i] = MTree[i].Weight ;

/*        SetErreur("\n\n Noeud %d : RootCell  = %d\n",i,MTree[i].RootCell);
        SetErreur("            nchildren = %d\n",MTree[i].nchildren);
        SetErreur("            weight    = %f\n",MTree[i].Weight);
        SetErreur("            position = %d / %d\n",MTree[i].CurrentPosition + 1,MTree[i].RangeLength);
        for(j=0;j<MTree[i].nchildren;j++)
            SetErreur("            Child[%d] = %d\n",j,MTree[i].Child[j]->RootCell);
*/
        free(MTree[i].Child);
        free(MTree[i].CWeight);
        free(MTree[i].PositionRange);
    }

    free(MTree);

    return(0);
}
