#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 , float * _cumul , set<int> & Constraints , set<int> & Nodes , float * _downcontrib )
{
    int nup,i,nr,nc ;
    int * UPSTREAM ;
    set<int>::iterator isi;
    float downcontrib_i;

    if ( GetUpstream ( F , _root , & nup , & UPSTREAM ) )
    {
            SetErreur("SegmentationTopo_ : echec de recuperation des cellules amont\n");
            return(-1);
    }

    (void) rec2RC ( &(R->Header) , _root , & nr, & nc );
    (*_downcontrib) = 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_ : echec de l'appel recursif\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 ( *_root );
         (*_downcontrib) = 0;
    }

    return(0);
}

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

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

     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 ;
    float downcontrib ;
    int i;
    set<int> Nodes ;
    set<int> Constraints;
    set<int>::iterator isi ;

    if ( !(R = GetFloatRaster(HandleRaster)) )
    {
        SetErreur("RasterFindTransition : le raster # %d n'existe pas\n",HandleRaster);
        return(-1);
    }

    if ( !(F = GetFlowDirGrid(HandleFlowDir)) )
    {
        SetErreur("RasterFindTransition : la grille d'ecoulement # %d n'existe pas\n",HandleFlowDir);
        return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("RasterFindTransition : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

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

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

    if ( SegmentationTopo_ ( F , R , & exu , & cumul , Constraints , Nodes , & downcontrib ) )
    {
        SetErreur("SegmentationTopo : echec de l'appel initial au pixel %d\n",exu);
        return(-1);
    }
 

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

    for(isi=Nodes.begin();isi!=Nodes.end();isi++)
    {
        (*_NODES)[i] = (*isi);
        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 : la grille %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

    if( !(R = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("MakeFlowLineThroughPoint : le raster %d n'existe pas\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("MakeFlowLineThroughPoint : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\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 : echec d'extraction de la partie amont de la ligne d'ÃƒÆ’Ã†â€™Ãƒâ€ Ã¢â‚¬â„¢ÃƒÆ’Ã¢â‚¬Å¡Ãƒâ€šÃ‚Â©coulement\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_ : echec de recuperation des cellules amont\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_ : echec de l'appel recursif\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> & UpNodes, int ** _BB , float * _Sum , set<int> & ixUpNodes , int * _nr , int * _nc )
{
    int i , nup1 ;
    float Sum_i ;
    set<int>::iterator isi;
    int * UPSTREAM;

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("GetCatchmentSum_ : echec de recuperation des cellules amont\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 = UpNodes.find(UPSTREAM[i]);

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

            (*_Sum) += Sum_i ;
        }
        else
        {
            ixUpNodes.insert ( (int) distance ( UpNodes.begin() , isi ) );
        }
    }

    free(UPSTREAM);

    return (0);
}

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

int GetCatchmentSum ( int HandleFlowDir , int HandleRaster , int outlet, int nupnodes, int * UPNODES , int ** _BB , float * _Sum , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    FloatRaster * R;
    set<int> UpNodes,ixUpNodes;
    set<int>::iterator isi;
    int nr;
    int nc;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetCatchmentSize : la grille %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

        if( !(R = GetFloatRaster(HandleRaster) ) )
    {
        SetErreur("MakeFlowLineThroughPoint : le raster %d n'existe pas\n",HandleRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(R->Header) ) )
    {
        SetErreur("MakeFlowLineThroughPoint : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleRaster, HandleFlowDir);
        return(-1);
    }

    *_Sum = 0;

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

    UpNodes.clear();
    ixUpNodes.clear();

    for (i=0;i<nupnodes;i++)
    {
        UpNodes.insert( UPNODES[i] );
    }

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

    if ( GetCatchmentSum_ ( F , R , & outlet , UpNodes, _BB , _Sum , ixUpNodes , & nr , & nc ) )
    {
        SetErreur("GetCatchmentSum : echec de l'extraction recursive\n");
	return(-1);
    }

    // Recopie des elements du set ixUpNodes dans le tableau IXUPNODES

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

    i=0;
    for (isi=ixUpNodes.begin();isi!=ixUpNodes.end();isi++)
    {
        (*_IXUPNODES)[i] = (*isi);
        i++;
    }

    // Vidage des containers

    ixUpNodes.clear();
    UpNodes.clear();

    return(0);
}

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

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

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("GetCatchmentSum_ : echec de recuperation des cellules amont\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 = UpNodes.find(UPSTREAM[i]);

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

            (*_Sum) += Sum_i ;
        }
        else
        {
            ixUpNodes.insert ( (int) distance ( UpNodes.begin() , isi ) );
        }
    }

    free(UPSTREAM);

    return (0);
}

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

int GetCatchmentSum1 ( int HandleFlowDir , int outlet, int nupnodes, int * UPNODES , int ** _BB , float * _Sum , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    FloatRaster * R;
    set<int> UpNodes,ixUpNodes;
    set<int>::iterator isi;
    int nr;
    int nc;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetCatchmentSum1 : la grille %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

    *_Sum = 0;

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

    UpNodes.clear();
    ixUpNodes.clear();

    for (i=0;i<nupnodes;i++)
    {
        UpNodes.insert( UPNODES[i] );
    }

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

    if ( GetCatchmentSum1_ ( F , & outlet , UpNodes, _BB , _Sum , ixUpNodes , & nr , & nc ) )
    {
        SetErreur("GetCatchmentSum : echec de l'extraction recursive\n");
	return(-1);
    }

    // Recopie des elements du set ixUpNodes dans le tableau IXUPNODES

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

    i=0;
    for (isi=ixUpNodes.begin();isi!=ixUpNodes.end();isi++)
    {
        (*_IXUPNODES)[i] = (*isi);
        i++;
    }

    // Vidage des containers

    ixUpNodes.clear();
    UpNodes.clear();

    return(0);
}

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

int Accumulate ( int HandleFlowDir, int HandleLocalRaster , int exu , int HandleCumulRaster )
{
    FlowDirGrid * F;
    FloatRaster * Loc , * Acc;
    int nr,nc;
    float Sum;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("Accumulate : la grille # %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

    if( !(Loc = GetFloatRaster(HandleLocalRaster) ) )
    {
        SetErreur("Accumulate : le raster # %d n'existe pas\n",HandleLocalRaster);
	return(-1);
    }

    if( !(Acc = GetFloatRaster(HandleCumulRaster) ) )
    {
        SetErreur("Accumulate : le raster # %d n'existe pas\n",HandleCumulRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(Loc->Header) ) )
    {
        SetErreur("Accumulate : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleLocalRaster, HandleFlowDir);
        return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(Acc->Header) ) )
    {
        SetErreur("Accumulate : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleCumulRaster, HandleFlowDir);
        return(-1);
    }

    if ( Accumulate_ ( F , Loc , Acc , & exu , & Sum , & nr , & nc ) )
    {
        SetErreur("Accumulate : echec de la construction recursive, cellule %d\n",exu);
        return(-1);
    }

    return(0);
}

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

int Accumulate_ ( FlowDirGrid * F , FloatRaster * Loc , FloatRaster * Acc , int * _outlet , float * _Sum , int * _nr_tmp , int * _nc_tmp )
{
    int i , nup ;
    float Sum_i,f ;
    int * UPSTREAM;

    (*_Sum) = 0;

    if ( GetUpstream ( F , _outlet , & nup , & UPSTREAM ) )
    {
            SetErreur("Accumulate_ : echec de recuperation des cellules amont, cellule %d\n",*_outlet);
            return(-1);
    }

    for(i=0;i<nup;i++)
    {
        if ( Accumulate_ ( F , Loc , Acc , UPSTREAM + i , & Sum_i , _nr_tmp , _nc_tmp ) )
        {
            SetErreur("Accumulate_ : echec de l'appel recursif\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 Accumulate1 ( int HandleFlowDir, int exu , int HandleCumulRaster )
{
    FlowDirGrid * F;
    FloatRaster * Acc;
    int nr,nc;
    float Sum;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("Accumulate : la grille # %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

    if( !(Acc = GetFloatRaster(HandleCumulRaster) ) )
    {
        SetErreur("Accumulate : le raster # %d n'existe pas\n",HandleCumulRaster);
	return(-1);
    }

    if ( CheckSameHeader( &(F->Header) , &(Acc->Header) ) )
    {
        SetErreur("Accumulate : le raster # %d et la grille d'ecoulement # %d n'ont pas la meme reference spatiale\n",HandleCumulRaster, HandleFlowDir);
        return(-1);
    }

    if ( Accumulate1_ ( F , Acc , & exu , & Sum , & nr , & nc ) )
    {
        SetErreur("Accumulate : echec de la construction recursive, cellule %d\n",exu);
        return(-1);
    }

    return(0);
}

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

int Accumulate1_ ( FlowDirGrid * F , FloatRaster * Acc , int * _outlet , float * _Sum , int * _nr_tmp , int * _nc_tmp )
{
    int i , nup ;
    float Sum_i,f ;
    int * UPSTREAM;

    (*_Sum) = 1;

    if ( GetUpstream ( F , _outlet , & nup , & UPSTREAM ) )
    {
            SetErreur("Accumulate_ : echec de recuperation des cellules amont, cellule %d\n",*_outlet);
            return(-1);
    }

    for(i=0;i<nup;i++)
    {
        if ( Accumulate1_ ( F , Acc , UPSTREAM + i , & Sum_i , _nr_tmp , _nc_tmp ) )
        {
            SetErreur("Accumulate_ : echec de l'appel recursif\n");
            return(-1);
        }

        (*_Sum) += Sum_i ;
    }

    free(UPSTREAM);

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

    return(0);
}

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

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

    if ( GetUpstream ( F , _outlet , & nup1 , & UPSTREAM ) )
    {
            SetErreur("GetUpstreamNodes_ : echec de recuperation des cellules amont\n");
            return(-1);
    }

    // Boucle sur les pixels amont

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

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

        if ( isi == UpNodes.end() )       // Si le pixel UPSTREAM[i] n'est pas trouve dans la liste des noeuds d'arret
        {
            if ( GetUpstreamNodes_ ( F , UPSTREAM + i , UpNodes , ixUpNodes ) )     // Appel recursif
            {
                SetErreur("GetUpstreamNodes_ : echec de l'appel recursif\n");
                return(-1);
            }
        }
        else
        {
            ixUpNodes.insert ( (int) distance ( UpNodes.begin() , isi ) );
        }
    }

    free(UPSTREAM);

    return (0);
}

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

int GetUpstreamNodes ( int HandleFlowDir , int outlet, int nupnodes, int * UPNODES , int ** _IXUPNODES , int * _nixupnodes )
{
    int i;
    FlowDirGrid * F;
    set<int> UpNodes,ixUpNodes;
    set<int>::iterator isi;

    if( !(F = GetFlowDirGrid(HandleFlowDir) ) )
    {
        SetErreur("GetCatchmentSum1 : la grille %d n'existe pas\n",HandleFlowDir);
	return(-1);
    }

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

    UpNodes.clear();
    ixUpNodes.clear();

    for (i=0;i<nupnodes;i++)
    {
        UpNodes.insert( UPNODES[i] );
    }

    if (  GetUpstreamNodes_ ( F , & outlet , UpNodes, ixUpNodes ) )
    {
        SetErreur("GetUpstreamNodes : echec de l'extraction recursive\n");
	return(-1);
    }

    // Recopie des elements du set ixUpNodes dans le tableau IXUPNODES

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

    i=0;
    for (isi=ixUpNodes.begin();isi!=ixUpNodes.end();isi++)
    {
        (*_IXUPNODES)[i] = (*isi);
        i++;
    }

    // Vidage des containers

    ixUpNodes.clear();
    UpNodes.clear();

    return(0);
}

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

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

    (*_NEWNODES) = NULL ;

    if ( !(CW = GetFloatRaster(HandleCumulRaster)) )
    {
        SetErreur("OptimMeshVariance : le raster # %d n'existe pas\n",HandleCumulRaster);
        return(-1);
    }

    if ( !(PK = GetFloatRaster(HandlePkRaster)) )
    {
        SetErreur("OptimMeshVariance : le raster # %d n'existe pas\n",HandlePkRaster);
        return(-1);
    }

    if ( !(F = GetFlowDirGrid(HandleFlowDir)) )
    {
        SetErreur("OptimMeshVariance : la grille d'ecoulement # %d n'existe pas\n",HandleFlowDir);
        return(-1);
    }


    if ( MakeMeshTree ( F , CW , PK , nnodes , NODES , CNT , & MTree , & meanweight ) )
    {
        SetErreur("OptimMeshVariance : echec de création de l'arbre de mailles\n");
        return(-1);
    }

    for(i=0;i<nnodes;i++)
    {
        if ( MTree[i].cnt | !(MTree[i].Parent) )
        {
            if ( OptimMeshVariance_ ( MTree + i , & meanweight ) )
            {
                SetErreur("OptimMeshVariance : echec de l'optimisation à l'amont du noeud %d\n", MTree[i].RootCell );
                return(-1);
            }
        }
    }

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

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

    return(0);
}